Category: Our tool & Methodology

R – Python

Posted on by sojung lucia KIM in R

한동안 연구소 일 때문에 Python 를 쓰다가 다시 논문으로 R 코드를 일년만에 들여다 보니, R-Studio 환경과 깔끔한 명령어 타이핑이 주는 이 친숙함으로인해 다시 python 환경으로 돌아가고 싶지 않아졌다.. 나 같은 사용자를 위해

현재 쓰고 있는 우분투 환경에서 최대한 R 과  Python 의 장점을 취할 수 있는 방법을 정리한다.

 

첫째, 라이브러리간의 dependency  를 유지하면서 프로젝트별로 안정하게 개발하려면 가상환경을 적극 활용한다.

  • 몇년 사이에 아나콘다는 정말 많이 무거워졌다. 그래도  dependency 문제로 python 이라면 아나콘다 의 conda 를 통해 라이브러리를 관리하는 것이 이후에 발생가능한 골치아픈 문제를 많이 해결해준다.
  • 아나콘다에서 R 의 주요 라이브러리까지 포함하는 R-essential  라이브러리를 설치한다.
  • anaconda-navigator 를 이용해서 관리하고 설치하고 런칭할 수도 있다.
  • 예를 들어 R 가상 환경을 만들고 ,  conda activate R_env
  • 일단 R 가상환경에서는 Rstudio 에서 패키지 설치시 , install.packages(“라이브러리명”)을 실행하더라도 anaconda conda 프레임에서 동작한다.  문제는 클라우드에서 제공되지 않는 R 라이브러리가 많다는 것이다. 이경우에는 git-hub의 라이브러리를 다운받아야 하는데 devtools::install_gitbuhb(“깃허브명/라이브러리명”)으로 다운을 받는다.  간혹 unzip 시에 tar 디렉토리를 못찾는 경우가 있다. 이 경우는 다음과 같이 리눅스의 디렉토리를 연결시켜준다.   : options(unzip = "internal")   sudo ln -s /bin/tar /bin/gtar   아나콘다 클라우드에 존재하는 라이브러리인 경우 터미널에서 다음과 같이 설치할 수도 있다. conda install -c r r-devtools

둘째, Python 사용이 개발용이 아니라 연구용이라면, 변수에 대한 glimpse,  디렉토리 구조까지 포함하는 Rstudio와 가장 비슷한 환경의 Spyder 의 사용을 고려한다.

  • spyder 의 환경은 R-Studio 와 비슷하지만, (1) 로딩시간이 많이 든다 (2) plot 과 관련된 인터페이스는 spyder4 에서 다소 개선되었다고 해도 앞으로 많이 개선되었으면 하는 부분이다  (3) 사용되는 메모리를 보여주는 기능은 훌륭하다. 그런데 메모리가 부족한 환경에서는   plot 을 보여주지도 못한다.

셋째,  우분투환경에서 R에서 한글을 처리한다면 Rstudio 가 아니라 Notebook 을 고려한다.

  • Python3 부터는 unicode 가 디폴트로, 이것만으로도 한글처리에 있어서 많은 문제가 해결되었다(euc-KR 비호환 등). R 은 이 부분에서는 대응이 거의 없는 편인듯 하다. 단일 윈도우즈 컴퓨터 환경 사용자라면 이 부분이 얼마나 심각한지 모를 수 있다. 대부분의 처리를 윈도우즈에서 해주기 때문. 문제는 리눅스에서 R 을 사용하는 경우이다. 다행히도 리눅스의 한글처리를 위한 많은 (그리고 정말 고마운) IME가 있다. 덕분에 R 콘솔에서 한글 입력은 문제가 없지만, 우분투에서 실행된 R-studio 에서는 한글입력이 불가능하다. 윈도우즈에서 작업한 RData 역시 우분투 R stdio 에서는 제대로 보이지 않는 경우가 많다.
  • 여기에 대한 하나의 솔루션은 Jupyter Notebook 을 사용하여 한글입출력을 브라우저에서 처리할 수 있도록 하는 것..

 

우분투서버에서 실행해본 결과 R의 실행시간과 IDE 로딩 속도가 python 대비 느리지는 않다 (오히려 빠른 편. ).  tensorflow 와 같은 프레임워크가 제공하는 gpu를 통한 메모리 최적화 기능이 가능한지는 확인해봐야 한다.

딥러닝을 포함 엄청난 속도로 기여하는 오픈소스 커뮤니티도 없다. 반면 주기적으로 발표되는 통계 논문과 그에 따른 라이브러리가 쏟아져 나온다.  아나콘다 파운데이션처럼 라이브러리 의존성을 관리해줄 주체도 없다. 이런 오만함에, 소위 비-ascii 처리는 별로 지원해줄 것 같다. R 사용자는 ‘최소한의 도덕성이 없다’고 모 구글 팀장이 말할 정도로, 우리가 프로그래밍에서 기대하는 철학이 별로 안 보인다.  이런 단점에도 불구하고,  해들리윅햄의 독선적이고도 우아한 ggplot 문법, data.frame, tibble,  등 데이터 프레임에 잘 맞는 직관적인 사용법을 가진 라이브러리는 놓치기 싫은 자산이다.

 

 

 

 

 

 

 

R 메모리 관리

Posted on by sojung lucia KIM in R, Uncategorized

R은 – C로 구현된 다양한 라이브러리를 제공한다. 빅데이터 처리에 처리 속도와 메모리 부족 문제를 자주 경험하고  있는 사람은 다음을 참고하자.

해들리 윅햄

http://adv-r.had.co.nz/memory.html

 

벡터는 4바이트 메모리를 차지한다. (C integer 는 4바이트 된다).

따라서, 32비트 운영체제라면  24 × 8 − 1 (231, about two billion) elements  을 저장할 수 있다 (C도 비슷하다). R 3.0  이상은 252 까지 지원된다.

The length of the vector (4 bytes). By using only 4 bytes, you might expect that R could only support vectors up to 24 × 8 − 1 (231, about two billion) elements. But in R 3.0.0 and later, you can actually have vectors up to 252elements. Read R-internals to see how support for long vectors was added without having to change the size of this field.

만약 벡터가 128 바이트 이상이라면 8 바이트씩 allocate 한다.

 

 

메모리에 여러개의 오브젝트를 동시에 올려 놓고 작업하다보면 내가 사용하고 있는 메모리가 궁금해진다.   pryr  라이브러리의 다음 함수가 유용하다.

각각의 오브젝트 사이즈는

object_size

R이 사용하는  전체 메모리는 mem_used()

메모리 변화를 보려면 mem_change

object_size (x <- 1:1e6) )
x <- 1:1e6  ### 1billlion
object_size (x)
mem_used()
mem_change(x)

dynamic programming 의 유래

Posted on by sojung lucia KIM in gossip

연구하기 좋지 않은 시절 과학자가 살아 남는 법

 

“전쟁 직후인 1950년대, 만연한 병리학적 두려움과 혐오로 한가로이 수학을 연구하기 좋지 않은 시절, 국방부는 ‘연구’라는 말을  꺼려했고 의회는 이런 혐오스런 단어에 대한 예산 승인을 거부했다.

그래서 나는 “프로그래밍”이라는 단어를 사용하기로 결정했다.. 그리고 아이디어간의 across 가 있으면서, 다이내믹 한것, 즉   multi-stage . 그런데 혐오감 주지 않은 단어..   그래 dynamic  programing 의 조합이 다른 단어랑 결합되서 혐오감을 주지 않을 것 같다. 의회도 거절하지 못하겠다 싶었다. 그래서 나는 dynamic programming 이라는 표현을 쓰기로 결정했다. ”

  • Richard Bellman (허리케인의 눈 : 자서전)

 

…The 1950s were not good years for mathematical research. [the] Secretary of
Defense …had a pathological fear and hatred of the word, research…
I decided therefore to use the word, “programming”.
I wanted to get across the idea that this was dynamic, this was multistage… I thought,
let’s … take a word that has an absolutely precise meaning, namely dynamic… it’s
impossible to use the word, dynamic, in a pejorative sense. Try thinking of some
combination that will possibly give it a pejorative meaning. It’s impossible.
Thus, I thought dynamic programming was a good name. It was something not even a
Congressman could object to.”
Richard Bellman, “Eye of the Hurricane: an autobiography” 1984.

source : https://web.stanford.edu/class/cs124/lec/med.pdf

common methods variance bias (Lindell & Whitney 2001)

Posted on by sojung lucia KIM in recipe, Uncategorized

아주 오래 전 일이다. 내 논문에서 응답자가 단일인 경우 CMV 처리에 대한 부분이 지적되었다.  당시에는 한개의 요인의 누적 분산량을 리포팅하는 harman’s one factor analysis로 처리를 했고 넘어갔지만 뒤에 marker variable (이론적으로 독립/종속 변수와 아무 상관이 없는 변수) 을 모델에 포함하는 것이 더 일반적이라는 것을 알게 되었다 (lindell & whitney 2001)  몇년뒤 작성한 후속 논문에선 수퍼바이저(NUS 교수)의 충고에 따라 marker variable을 사용했다.  그리고 놀랍게도 path 에 변화가 생기고,  고심 끝에 이 부분은 보고하지 않은채  출간이 되었다.

10년이 지난 지금, 내가 새롭게 수집한 자료에 아주아주 심한 CMV 가 의심되지만 처리하지 않고 넘어간 덕에 ISR 리뷰어에게 다시 지적 받았다.    CMV 처리를 하는 것 자체는 어렵지 않은데, 도대체 어떤 응답자가 어떤 이유에서 더 많은 bias를 보여줄까?

posakoff 는 CMV 에 대해  아주 집요하게 연구한 학자이며 일반적으로 CMV에 대한 처리는 그의 논문을 인용한다. 본 포스팅에선 좀더 쉽게 정리한 페이퍼를 소개한다. 여기서  marker  variable 을 활용한 방법을 참조하자.

http://www.galileoco.com/comSciJliterature/vishwanathEgnotoOrtega12.pdf

  • one factor analysis : unrotated single  factor의 variance 가 majority 인지 여부를 확인하는 것
  • marker variable : rM 을 활용하여 바람직한  r A 을 찾는 방법

r M : marker  와 다른 변수들간의 correlation

r U : un-restricted model  correlation

cmv

Lindell & Whitney (2001)  : 이런 종류의 top-tier 논문이 그러하듯  Journal of applied psychology 에 출간된 논문이다.  Lindell 는 당시 소속이  Texas A&M Univ. Hazard reduction & recovery center.  Whiteny 는 심리학과 교수이다.

모델은 X1, X2 가 Y 에 영향을 주는 경로를 표현한다.  r* 이 latent 변수들간의 correlation 이라면 r 은 observed 변수들간의 correlation 이다.

cmv_lindell_model

Lindell  & Whitney (2001) 의 equation 4번의 좌측항이 rA (Marker 가 통제된 상황에서의 Xi – Y 의 correlation)

cmv_lindell

 

reverse causality

Posted on by sojung lucia KIM in recipe, Uncategorized

OLS 의 기본 가정은 독립변수 X 와 오차항 E 사이에 상관관계가 없어야 한다는 것이 있다.

OLS 로 테스트함에 있어 종속변수가 다시 독립변수에 영향을 줄 수 있다면 reverse causality 의 문제가 생긴다.

 

다음 예를 보자

Y= XB + E        // 이러한 식이 있다.

X = YC + W     // 그런데 X 는 Y 에 의해 영향을 받는다.

X = (XB + E) C + W  // X 와 E 는 correlated 되어 있고, 이는 다시 OLS 가정에 위배된다.

만약 이런 문제 (endogenity, reverse causality)가 존재한다면 가설 검정에도 문제가 있다.

 

 

<1> ad hoc solution

만약 종속변수가 잠재적으로 endogenous 하다면 그런 문제가 생기지 않을 proxy를 찾아야 한다.

If a dependent variable is potentially endogenous, it is intuitively appealing to look for a proxy that does not suffer from the same problem.

일반적인 방법은, 문제가 될 여지가 있는 변수를 찾아, 한 두 period 지연된  변수를 사용하는 법

 

<2> IV (instrumental variables) technique  (도구 변수 방법)

위의 X를 둘로 나눈다. E 와 관련된 것, E 와 관련되지 않은 것

E 와 관련되지 않은 것을 사용하면 올바른 파라메터 추정을 할 수 있다 (위의 예시 – 계수 B (beta) 를 찾는다)

가장 일반적인 IV 추정방법은  Two Stage Least Squares (TSLS) 이다.

Intuitively, IV estimation works as follows:

잠재적으로 내생성을 가진 regressor 와 강력한 상관관계를 가진 변수를 찾은 후 (X 와 밀접한 관계를 가진 IV 라 하자,  IV 는 오차항과는 상관관계가 없어야 한다)

Find a genuinely exogenous variable (instrument) that is strongly correlated with the potentially endogenous regressor.

그 변수가 잠재적으로 내생성을 가진 독립변수를 통해 종속변수에 영향을 준다는 것을 보여준다.

IV ->X ->  Y

Ensure that the instrument only influences the dependant variable through the potentially endogenous independent variable.

위에서 언급한 바와 같이 2단계 방법은 다음과 같다.

First Stage :  X 를 Z 에 대해 회귀분석하여  X 의 추정치를 찾는다 (이를 X_hat 이라 한다. X_hat 역시 n개의 관측치를 가진 vector 이다)

X = b0 + b1 * Z   // 이 회귀방정식을 통해 X_hat을 구한다.

Second Stage : X_hat 을 원래 변수로 대치하여 추정치를 분석하여 beta 를 찾는다.

Y = B0 + B1 * X_hat  //  이 방정식에서 B1이 내생성이 통제된 베타값이 된다.

 

 

R에서 도구변수를 위해 많이 사용되는 패키지는 AER(Applied Economics with R) 이다.

예를 들어 ,  ITG 데이터로 아래와 같은 regression 을 만들었다.

fit = lm (performance ~ governance, data = ITG)

거버넌스가 잘되면 성과가 높아진다. 반대로 성과는 다시 거버넌스의 효과를 높인다.  이 경우 오직 거버넌스에만 강력한 영향을 주는 변수를 찾는다.  공식적 미팅은 거버넌스에만 영향을 줄 것이다.

fit1 = lm (governance ~ communication, data = ITG)

여기에서 나온 fit1$fitted.value 가 바로 X_hat 이다.

fit2 = lm (governance ~ hat_governance, data = cbind(ITG, data.frame(hat_governance =  fit1$fitted.value ))

이 복잡한(?) 과정을 거치지 않고 AER 의 ivreg 를 사용한다.

fit2= ivreg(performace ~ governance | communication, data = ITG)

 

 

 

 

 

 

 

 

 

 

 

Heteroscedasticity

Posted on by sojung lucia KIM in terminology

it refers the circumstannce in which the variablity of a variable (y) is unequal across the rage of values of a second variable (x) that predict it (y).

scatterplot을 그려보면 cone-like 모양을 그린다.  예를 들어 나이에 따른 소득은, 어릴때는 나이가 어릴 수록 소득이 높은 구조이다. 하지만 30대 40대 가 되어가면 전혀 다른 carrier path를 가지게 되고 변동성이 더 커진다.  크게 보면 나이가 늘 수록 소득이 높아지지만, 그 편차가 커진다.

 

regression 에 heterosedasticity 여부를 확인하는 방법으로 white test 가 있다.

bstats::white.test  . (cran에서 더이상 제공되지 않는다)

het.test::white.test (VAR object를 인수로 사용한다. VAR를 이용하기 위해서는 vars 패키지가 우선 설치되어 있어야 한다.

lmtest::bptest

Breusch-Pagan test ( bptest)를 수행해서 p value 가 낮다면 H1 : hetero 하다..다음을 확인해보자.

If a model is estimated using the following code:

 lm(y~x1+x2)->p

1. bptest(p) does the Breuch Pagan test to formally check presence of heteroscedasticity. To use bptest, you will have to call lmtest library.

2. If the test is positive (low p value), you should see if any transformation of the dependent variable helps you eliminate heteroscedasticity. Also check if the right hand side of the model is okay.

3. If 2 does not work, you can use the white’s heteroscedasticity-corrected covariance matrices to make inference. Package car has a function hccm that gives you the heteroscedasticity-corrected covariance matrix (there is a similar function in package sandwich also). coeftest(p,vcov=hccm(p)) will give you the results of the tests using this matrix. Use these results instead of summary(p).

library(lmtest)
library(car)
bptest(p)
coeftest(p,vcov=hccm(p))

 

Simpson’s paradox

Posted on by sojung lucia KIM in terminology

https://en.wikipedia.org/wiki/Simpson%27s_paradox

통계분석할때 흔히하는 오류의 하나이다.  그룹별로 데이터 분석을 하면 드러나지만 그룹이 함께 있으면 사라지거나 반대가 되는 현상

에드워드 심슨이 이러한 현상을 발표했지만 이전에 Pearson 과 Yule 가 이미 비슷한 효과를 주장한적이 있다.

1973년 UC버클리대학이 지원한 대학원생을 뽑을때  남성과 여성의 차별을 하고 있다는 소송을 당했는데, 모든 학과별로 종합된 데이타에 의하면 남성의 admitted 비율이 여성보다 높다.

Applicants Admitted
Men 8442 44%
Women 4321 35%

 

하지만 학과별로 살펴보면 통계적으로 유의하게 여성의 진학률이 더 높은  더 선호하는 학과가 더 많았다.

원인은 여성은 경쟁이 몰려 admitted 비율이 낮은 학과에 지원하는 경향이 높은 반면 (영어 등, 남성 300명 여성 500명 지원, admitted 비율이 30%),

남성은  admitted 비율이 높은 학과에 지원하는 경향이 높았다 (공학 등 남성 800명, 여성 100명 지원 비율은 admitted 비율이 80%)

 

 

Department Men Women
Applicants Admitted Applicants Admitted
A 825 62% 108 82%
B 560 63%  25 68%
C 325 37% 593 34%
D 417 33% 375 35%
E 191 28% 393 24%
F 373  6% 341  7%