윈도우에서 R 업데이트 하기 

아래 방법으로 최신 버전의 R을 설치할 수 있습니다. 

install.packages("installr")

library(installr)

updateR()

이전 버전의 패키지 가져오기

1. 이전 R 버전 폴더의 모든 라이브러리들을 새로운 R 버전 폴더로 copy and paste 합니다.

저는 R 3.4.* 버전에서 R 3.5.3. 버전으로 업데이트하였는데, 경로는 아래와 같습니다. 

C:\Users\Username\Documents\R\win-library\3.4

C:\Users\Username\Documents\R\win-library\3.5

3.4 폴더에 있는 모든 라이브러리 폴더들을 3.5 폴더로 복사 붙여넣기하면 됩니다. 하지만 이렇게만 하면 패키지 버전이 호환되지 않기 때문에 제대로 로드가 안될 수도 있습니다. 따라서 3.5 폴더에 있는 라이브러리들을 업데이트 시켜줘야합니다.

2. 패키지 업데이트

update.packages(checkBuilt=TRUE)

를 입력하면 팝업 박스가 나오면서 yes no 를 선택할 수 있는데 yes 를 계속해서 클릭해줍니다. 

패키지 ‘ucminf’를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다

패키지 ‘triebeard’를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다

패키지 ‘ordinal’를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다

패키지 ‘magic’를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다

패키지 ‘lpSolve’를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다

패키지 ‘RcppProgress’를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다

패키지 ‘desc’를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다

패키지 ‘xopen’를 성공적으로 압축해제하였고 MD5 sums 이 확인되었습니다

...... 

그러면 이런식으로 패키지를 최신 버전으로 업데이트를 알아서 해줍니다. 

3. 패키지가 모두 제대로 설치가 되었는지 확인합니다. 

version

packageStatus()

이제 최신 버전의 R을 이용할 수 있습니다. 

R - 다중회귀분석 (Multiple Linear Regression)

R 을 통해 다중회귀분석을 수행하는 기본적인 절차를 포스팅합니다. 

라이브러리 임포트

datarium 의 marketing 이라고 하는 데이터셋을 이용하여 다중회귀분석을 수행한다.

library(devtools)
# install_github("kassambara/datarium")

library(tidyverse)
library(datarium)

data("marketing", package = "datarium")
head(marketing, 4)

탐색적 자료분석 - 데이터 분포 살펴보기

가장 먼저 할 일은 데이터의 분포를 살펴보는 것이다. 

ggplot(data = marketing, aes(x = newspaper, y = sales)) +
  geom_point(size=5)

위와 같이 하나하나 산점도를 볼 수 있지만, 아래처럼 모든 변수들의 분포 (Variation)와 두 변수의 관계 (Covariation) 한 번에 볼 수 있도록 하는 라이브러리들이 많이 존재한다. 예를 들어, PerformanceAnalytics 라이브러리를 이용하면 아래와같이 한 번에 분포와 상관계수까지 구할 수 있다. 

library("PerformanceAnalytics")
chart.Correlation(marketing, histogram=TRUE, pch=19)

모델 만들기

sales = b0 + b1youtube + b2facebook + b3*newspaper 

분포를 봤을 때, youtube, facebook, newspaper 가 sales와 선형관계를 갖고 있다고 가정하고 위와 같이 다중 회귀 모델을 만들 수 있으며, 아래 R 코드를 통해 모델을 적합시킬 수 있다. 

model <- lm(sales ~ youtube + facebook + newspaper, data = marketing)
summary(model)

Call:
lm(formula = sales ~ youtube + facebook + newspaper, data = marketing)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.5932  -1.0690   0.2902   1.4272   3.3951 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.526667   0.374290   9.422   <2e-16 ***
youtube      0.045765   0.001395  32.809   <2e-16 ***
facebook     0.188530   0.008611  21.893   <2e-16 ***
newspaper   -0.001037   0.005871  -0.177     0.86    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.023 on 196 degrees of freedom
Multiple R-squared:  0.8972,	Adjusted R-squared:  0.8956 
F-statistic: 570.3 on 3 and 196 DF,  p-value: < 2.2e-16

결과 해석

Regression 결과를 해석하는데, 가장 먼저할 것은 F-statistics 를 확인하는 것이다. 가장 밑에 F-statistics 를 보면, p-value가 매우 작아, 최소 한 변수가 유의하다는 것을 알 수 있다. 어떤 변수가 유의한지를 확인하기 위해, 아래 coefficient table 을 보면 된다. 이 때, t-value와 p-value 는 귀무가설 평균=0 로 놓고 구한 값으로, 귀무가설로부터 얼마나 벗어난지에 관한 지표이다. 

X의 estimate 값은 다른 변수가 고정되어있을 때, X가 한 단위 변화했을 때, Y의 평균적인 변화량이다. 예를 들어, youtube의 estimate 는 0.04 인데, 이는 youtube 가 1 변화했을 때, sales의 평균적인 변화량이다. newpaper 의 경우, p-value > 0.05 이므로 이러한 estimate 의 값이 유의하지 않다는 것을 알 수 있다. 

t value 와 p-value 는 무엇인가?

$$ H_0 : \beta = 0 $$

$$ t = \frac{\hat{\beta}}{se(\hat{\beta})} $$ 

이 때 위의 t 값은 자유도가 n-2 인 t 분포를 따른다. (이렇게 되는 이유)

p-value는 이렇게 계수 추정치가 t 분포를 따른다는 사실을 이용해서 구한 값이다. 

결정계수와 F-stat 

다음과 같이 정의해보자.

$$ SSE = \sum^{n}_{i=1}{(y_i - \hat{y})^2} $$

$$ SSR = \sum^{n}_{i=1}{(\bar{y} - \hat{y_i})^2} $$

$$ SST = \sum^{n}_{i=1}{(y_i - \bar{y})^2} $$

이 때, 각각의 통계량에 대한 자유도는 아래와 같다. 

df(SSE) : n-p

df(SSR) : p-1

df(SST) : n-1

이 때, p는 intercept를 포함한 모델의 총 파라미터의 갯수이며, n은 샘플 수이다. 

$$ MSE = \sum^{n}_{i=1}{(y_i - \hat{y})^2} / (n-p) $$

$$ MSR = \sum^{n}_{i=1}{(\bar{y} - \hat{y_i})^2} / (p-1) $$

$$ MST = \sum^{n}_{i=1}{(y_i - \bar{y})^2} / (n-1) $$

귀무가설을 아래와 같이 설정하자

$$ H_0 : \beta_1 = \beta_2 = \beta_3 = ... = \beta_{p-1} $$ 

이 때, 

$$ \frac{MSE}{MSR} \sim F(n-p, p-1) $$

이를 통해 F-stat 과 p-value를 구할 수 있다. 이 값은 (explained variance) / (unexplained variance) 가 F-stat 을 따른다는 사실을 이용한 검정 방법이라고 할 수 있다. 왜 이렇게 되는지에 대해서는 이 링크를 참고.

결정계수 (Coefficient of determination, Rsquare)

결정계수는 아래와 같이 정의되며, 전체 분산 중 모델에의해 설명되는 분산의 양을  의미한다. 

$$ R^2 = \frac{SSR}{SST} $$

adjusted 결정계수를 사용하기도 하는데, 이는 아래와 같이 정의된다. 

$$ \bar{R^2} = 1-(1-R^2) \frac {n-1}{n-p} $$ 

적합된 모델의 계수

summary(model)$coefficient

새로운 모델을 아래와 같이 구할 수 있다. sales = 3.5 + 0.045youtube + 0.187facebook.

새로운 모델

model  <- lm(sales ~ youtube + facebook, data = marketing)
summary(model)

Call:
lm(formula = sales ~ youtube + facebook, data = marketing)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.5572  -1.0502   0.2906   1.4049   3.3994 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.50532    0.35339   9.919   <2e-16 ***
youtube      0.04575    0.00139  32.909   <2e-16 ***
facebook     0.18799    0.00804  23.382   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.018 on 197 degrees of freedom
Multiple R-squared:  0.8972,	Adjusted R-squared:  0.8962 
F-statistic: 859.6 on 2 and 197 DF,  p-value: < 2.2e-16

신뢰구간

confint(model)

위와 같이 계수의 신뢰구간을 출력할 수도 있다. 이는 t 분포를 이용한다. 

로지스틱 회귀분석의 최대우도추정법 직관적 이해

Logistic regression의 Random component 

$$ Y \sim B(n, \pi) $$ 

Logistic regression은 반응변수 y가 위와 같은 이항 분포를 따른다고 가정한다.  

Model 

$$ \pi = \frac{1}{1+e^{-{(X\beta+\alpha)}}} $$ 

이 때, X를 통해 모수 pi 를 위와 같이 추정하는 것이다. 또는 일반적으로 아래와같이 표현되기도 한다. 

$$ logit(\pi) = X\beta + \alpha $$ 

Logistic regression에서 하고자하는건 결국, Random component에 주어진 pi를 추정하는 것인데, 이것이 0-1 사이의 값을 갖기 때문에 logit 변환을 수행하는 것이다. logit link function은 0~1사이의 값을 음의 무한대에서 양의 무한대로 바꾸어준다. 

Random component가 이항 분포이므로 likelihood는 아래와 같이 계산된다.

이곳에 Model을 집어넣음으로써, beta와 alpha 를 maximum likelihood 추정법으로 계산할 수 있게 된다. 근데 여기서 n과 y는 n번 시도했을 때 성공횟수 y를 나타낸다. 근데, 실제 각각의 데이터를 보면 시행횟수 n=1이며, y=0 or 1이다. 따라서 위 likelihood 는 베르누이 likelihood 와 같게된다. 

결국, 위의 likelihood를 최대화하면 계수 추정값을 얻을 수 있는데, log likelihood를 최대화하는 방식으로 쉽게 계수를 추정할 수 있다. 로지스틱 회귀분석에서는 보통 iteratively reweighted least squares (IRLS) 라는 방법을 통해 계수를 추정한다. 

$$ \sum_{i=1}^{N} y_i log(\pi_i) + (1-y_i) log (1-\pi_i) $$

Logistic regression의 log likelihood이다. 이는 negative binary cross entropy 이다. 결국, log likelihood를 최대화 하는 것은 negative binary cross entropy를 최대화 하는 것과 같고 이는 binary cross entropy를 최소화화는 것과 같다. 

고급통계 - Net reclassification improvements

Net reclassification index란 예측 모형에서 새로운 모델이 과거 모델보다 얼마나 얼마나 더 나아졌는지를 측정하는 지표이다. AUC 값을 기준으로 이를 판단하는게 과거에는 일반적이었으나, 그 한계점이 지적되외서 NRI가 등장하게 되었다. Framingham Heart Study에서 HDL 관련 연구를 수행할 때, AUC의 한계를 볼 수 있었는데, HDL을 예측 모형에 추가했을 때, AUC값이 크게 증가하지 않는 것이었다. 하지만 AUC 측면에서가 아니라, 모형의 계수를 통해 HDL의 effect를 파악할 때, HDL은 심장질환에 중요한 예측 인자였다. 

AUC의 한계는 다음과 같이 제시되었다. 1) 임상적인 의미와 연결시키기 힘들다. 2) 작은 변화는 잡아내기가 힘들다.

하지만 NRI는 이전 모델과 비교해 새로운 모델에서 얼마나 많은 사람들이 제대로 재분류되었는가를 파악함으로써 모델의 개선을 측정한다. 

출처 - https://en.wikipedia.org/wiki/Net_reclassification_improvement

만약 환자 30명, 비환자 70명에 대해서 환자-비환자 여부를 판단하는 모델 Test1, Test2를 만들었다고 하자. Test 2가 Test1에 비해 얼마나 우월한지를 판단하고자 한다. 위 테이블에서 여기서 검은색 글씨는 Test 1, Test 2 둘 다 맞춘거고, 하얀색 글씨는 Test1, Test 2 둘 다 틀린 경우, 초록색 글씨는 Test 2에서만 맞추고, Test 1에서는 틀린 것, 빨간색 글씨는 Test 1에서 맞추고, Test 2 에서 틑린 것이다. 

이 때, NRI의 계산은 매우 간단한데, Test 2를 통해서 얻은 이득이 얼마나 더 많은지를 계산해주면 된다. 이것은 Test 2를 도입했을 때, 추가적으로 맞춘것-오히려 틀린것 을 해줌으로써 계산되며, Event, Non-event 두 그룹의 대상자에 대해서 계산한 다음에 더함으로써 최종 NRI를 계산하게 된다. 

즉, Event가 발생한 사람들에 대한 NRI는 아래와 같이 계산된다.

$$ NRI_e = (4-2)/30 = 0.067 $$ 

다음으로 Non-event에 대한 NRI는 아래와 같이 계산된다.

$$ NRI_{ne} = (8-4)/70 = 0.057 $$ 

따라서 NRI는 아래와 같이 계산된다.

$$ NRI = NRI_e + NRI_{ne}= 0.124 $$ 

Gradient Boosting Algorithm의 직관적인 이해

실패를 통해 성공을 발전시켜라. 낙담과 실패는 성공으로 가는 가장 확실한 두 개의 디딤돌이다. -데일 카네기

Gradient Boosting Algorithm (GBM)은 회귀분석 또는 분류 분석을 수행할 수 있는 예측모형이며 예측모형의 앙상블 방법론 중 부스팅 계열에 속하는 알고리즘입니다. Gradient Boosting Algorithm은 Tabular format 데이터 (엑셀형태와 같이 X-Y Grid로 되어있는 데이터)에 대한 예측에서 엄청난 성능을 보여주고, 머신러닝 알고리즘 중에서도 가장 예측 성능이 높다고 알려진 알고리즘입니다. 그렇기 때문에 Gradient Boosting Algorithm을 구현한 패키지들이 많습니다. LightGBM, CatBoost, XGBoost 같은 파이썬 패키지들이 모두 Gradient Boosting Algorithm을 구현한 패키지들입니다. GBM은 계산량이 상당히 많이 필요한 알고리즘이기 때문에, 이를 하드웨어 효율적으로 구현하는 것이 필요한데, 위 패키지들은 모두 GBM을 효율적으로 구현하려고한 패키지들이라고 볼 수 있습니다.

Boosting 이란?

주요 앙상블 알고리즘은 bagging과 boosting으로 나눌 수 있고, Gradient boosting은 boosting 계열의 앙상블 알고리즘입니다. Gradient Boosting 알고리즘은 Gradient 를 이용하여 Boosting하는 알고리즘입니다. Gradient boosting 외의 boosting 계열의 알고리즘의 예로는 Adaptive boosting을 들 수 있습니다. 

먼저 Boosting의 개념부터 살펴봅시다. Boosting이란 약한 분류기를 결합하여 강한 분류기를 만드는 과정입니다. 분류기 A, B, C 가 있고, 각각의 0.3 정도의 accuracy를 보여준다고 합시다. A, B, C를 결합하여 더 높은 정확도, 예를 들어 0.7 정도의 accuracy를 얻는 게 앙상블 알고리즘의 기본 원리입니다. Boosting은 이 과정을 순차적으로 실행합니다. A 분류기를 만든 후, 그 정보를 바탕으로 B 분류기를 만들고, 다시 그 정보를 바탕으로 C 분류기를 만듭니다. 그리고 최종적으로 만들어진 분류기들을 모두 결합하여 최종 모델을 만드는 것이 Boosting의 원리입니다. 

GBM의 직관적인 이해

GBM을 이해하는 가장 쉬운 방법은 Residual fitting으로 이해하는 것입니다. 아주 간단한 모델 A를 통해 y를 예측하고 남은 잔차 (residual)을 다시 B라는 모델을 통해 예측하고 A+B 모델을 통해 y를 예측한다면 A보다 나은 B 모델을 만들 수 있게 되죠. 이러한 방법을 계속하면 잔차는 계속해서 줄어들게되고, training set을 잘 설명하는 예측 모형을 만들 수 있게 됩니다. 하지만 이러한 방식은 bias는 상당히 줄일 수 있어도, 과적합이 일어날 수도 있다는 단점이 있습니다. 따라서 실제로 GBM을 사용할 때는 sampling, penalizing 등의 regularization 테크닉을 이용하여 더 advanced 된 모델을 이용하는 것이 보편적입니다. 하지만 본 포스팅의 목적은 GBM 에 대하여 직관적인 이해를 해보는 것이기 때문에 이 부분은 다음에 기회가 있을 때 정리해보도록 하겠습니다. 

위 그림을 보시면 tree 1을 통해 예측하고 남은 잔차를 tree2를 통해 예측하고, 이를 반복함으로서 점점 잔차를 줄여나가는 것을 볼 수 있습니다. 이 때, 각각의 모델 tree1,2,3 을약한 분류기 (weak learner), 이를 결합한 분류기를 강한 분류기 (strong learner)라고도 합니다. 보통 약한 분류기로는 간단한 의사결정나무 (decision tree)를 많이 사용합니다. 이를 Gradient boosting tree라고도 하는데, 구현한 대표적인 라이브러리로 XGboost를 들 수 있습니다. XGBoost는 python, R로 이용할 수 있습니다. 

Gradient란 무엇인가?

그렇다면 residual fitting과 gradient가 무슨 관계인지 궁금하실텐데요. residual은 loss function을 squared error로 설정하였을 때, negative gradient 입니다. 따라서 residual에 fitting해서 다음 모델을 순차적으로 만들어 나가는 것은 negative gradient를 이용해 다음 모델을 순차적으로 만들어 나가는 것으로 볼 수 있습니다. 우선 아래 식을 통해 이를 확인해봅시다. 

loss function을 아래와 같이 정의하면,

$$ j(y_i, f(x_i)) = \frac{1}{2}(y_i - f(x_i))^2 $$ 

이 때의 negative gradient는 residual이 됩니다.

residual이 $$ y_i - f(x_i) $$ 이기 때문에 negative gradient가 residual이 된다는 것을 알 수 있습니다. 

따라서 gradient boosting은 다음 모델을 만들 때, negative gradient를 이용해 만들기 때문에 gradient boosting 이라고 할 수 있습니다. 또한 residual fitting model은 gradient boosting 모델의 한 가지 종류입니다. 왜냐하면 다른 loss function을 이용해서 gradient boosting을 할 수 있기 때문입니다. 예를 들어 회귀 문제가 아닌 분류 문제라면 loss function으로 squared error 를 이용할 수 없기 때문에 새로운 loss function을 만들어야하는데, 이 경우에도 negative gradient를 이용해서 새로운 model 을 fitting하고 이를 합산해 나가시는 식으로 최종 모델을 만들게 됩니다. 

직관적으로 이것은 무엇을 뜻할까요? 왜 negative gradient를 이용해서 새로운 모델을 만드는 걸까요? negative gradient는 pseudo-residual이라고도 불리며, 이것은 어떤 데이터 포인트에서 loss function이 줄어들기 위해 f(x)가 가려고하는 방향입니다. 이 방향에 새로운 모델을 fitting해서 이것을 이전 모델과 결합하면, f(x) 는 loss function이 줄어드는 방향으로 업데이트가 되겠죠. 이것이 gradient boosting의 아이디어입니다. Gradient boosting을 gradient descent + boosting 이라고 하기도 합니다. 어쨌든, loss function을 줄이는 방향의 negative gradient를 얻고, 이를 활용해 boosting을 하는 것이기 때문에 gradient descent와 boosting이 결합된 방법이다. 라고 이해하셔도 괜찮습니다. 

위 그림은 위키피디아에서 gradient boosting에 대하여 설명한 알고리즘도입니다. 이제 위 내용을 이해할 수 있습니다. gradient boosting에서는 learning rate를 통하여 pseudo-residual에 fitting된 모델을 어느정도로 update할지에 대해서 정하게 되는데, 이 값은 위 알고리즘도의 3번처럼, loss를 최소화하는 값을 optimization 과정을 통해 얻을 수도 있겠지만, 일반적으로 0.001~0.01 정도로 매우 작은 값으로 설정하는 것이 보편적입니다. 왜냐하면 작은 값으로 설정하여야 굉장히 세밀한 분류기를 얻을 수 있습니다. learning rate가 높을 수록 빠르게 모델의 bias를 줄여나가지만, learning rate가 적으면 디테일한 부분을 놓칠 수 있다고 이해하시면 됩니다. 이와 관련해서는 gradient boosting 과정을 simulation 해볼 수 있는 사이트를 소개합니다. 이 곳에서 직접 gradient boosting 방법이 어떻게 함수를 fitting 해나가는지를 시각적으로 확인할 수 있습니다. 

참고

https://machinelearningmastery.com/gentle-introduction-gradient-boosting-algorithm-machine-learning/

Tensorflow GPU 버전을 설치

설치버전

Cuda 9.0v

Cudnn 7.0.5v

tensorflow-gpu 1.9.0

GPU 

1080TI

운영체제

Ubuntu 16.04

설치방법

아래포스트를 기본으로 오류를 수정해나가면서 설치함

https://medium.com/@zhanwenchen/install-cuda-and-cudnn-for-tensorflow-gpu-on-ubuntu-79306e4ac04e

설치 팁

- Cuda, Cudnn을 먼저 설치하고, (이 때, Cuda Cudnn의 버전이 tensorflow 에서 지원하는지를 확인해야 함. 최신버전의 경우 tensorflow가 지원을 안하는 경우가 있을 수도 있음.) 이에 맞는 tensorflow-gpu를 설치하는 것이 좋다.

- Cuda, Cudnn, tensorflow-gpu는 서로 의존적이므로, 버전에 특히 신경 써야하고, 함부로 업그레이드 하지 않아야하고 현재 버전을 기록해 두어야한다. 

- 굳이 하나하나 Tool 별로 Document를 찾아가면서 설치하는 것보다, 먼저 설치한 사람의 설치 방법을 보고 설치하는 것이 훨씬 빠르고 안전하다. 

- 위 포스트에서 cuda와 cudnn이 설치 되었다는 것을 확인하는 방법으로 sample 코드를 실행하는 방법을 사용하였는데, 확실한 설치 여부를 알 수 있으므로 좋다. 

설치할 때 발생한 오류

1. libstdc++6 의 문제 

해결 : https://askubuntu.com/questions/671791/lib32stdc6-package-depends-on-gcc-base-but-my-installed-version-is-newer

2. tensorflow-gpu 버전 문제

pip install tensorflow-gpu를 통해 설치하니 tensorflow-gpu 1.12.0 버전이 설치됨. 근데, 이 버전에서는 구버전의 cudnn을 쓰지 않는 문제가 발생. 따라서 tensorflow-gpu를 1.9.0 버전으로 다운그레이드 함. 

설치 완료 후, Multi GPU 환경에서 default로 0번 gpu만 사용하는 문제

tensorflow를 import하기 전에 아래 문구를 삽입함으로써, 특정 gpu를 사용하도록 강제할 수 있다.

import os

os.environ["CUDA_VISIBLE_DEVICES"]= "2"   

import os

os.environ["CUDA_VISIBLE_DEVICES"]= ""  

tensorflow의 이전 방법에서는 아래 문구를 통해 특정 gpu를 사용할 수 있었는데 1.9.0 버전에서는 아래 문구가 안 먹힌다. 

R igraph 설치 오류 해결

문제 : libgfortran.so.4 오류로 인해 igraph 패키지가 설치되지 않음

해결 : 

sudo apt-get install gcc-4.8-base libgfortran4

library(devtools)

install_github("igraph/rigraph")

R Default library path 바꾸기 (공용 라이브러리 경로 설정)

.libPaths() 를 하면 현재 사용하고 있는 라이브러리 경로들이 나온다.

[1] "C:/Users/사용자이름/Documents/R/win-library/3.3"

[2] "C:/Program Files/R/R-3.3.2/library"

R은 여기서 첫 번째 element를 Default로 쓰는데, 이것은 보통 개인 폴더로 잡혀 있다. 

하지만 다중 사용자 환경에서는 내가 설치한 라이브러리를 남이 사용하도록 해야할 경우도 있는데 이때는 아래와 같이 하면 된다. 

1. R_LIBS 환경변수를 공용 라이브러리 경로로 설정함 

2. R에서 Sys.getenv("R_LIBS")를 통해 환경변수가 잘 잡혔는지를 볼 수 있다.

3. .libPaths()를 다시 입력하면 공용 라이브러리 경로가 1번으로 잡혀있는 것을 확인할 수 있다. 

[1] "C:/Program Files/R/R-3.3.2/library" 

[2] "C:/Users/사용자이름/Documents/R/win-library/3.3"

이후에는 패키지를 설치하면 모든 사용자가 이용할 수 있게 된다. 

Python PDPbox 패키지 설치시 문제 해결

문제

 PDPbox 패키지를 설치하던 도중 아래 에러가 발생하면서 설치가 완료되지 않음

/usr/bin/ld: cannot find -lpython3.5m

 collect2: error: ld returned 1 exit status

 error: command 'gcc' failed with exit status 1

해결

sudo apt install python3.5-dev

구글링해도 안나와서 올립니다.

도움을 얻은 페이지

https://github.com/giampaolo/psutil/issues/1143

생존분석 - 모수적 모형

Cox의 비례위험 모형

일반적으로 생존분석에 사용되는 모델인 Cox의 비례위험모형의 경우 준모수적 (Semi-parametric) 방법입니다. 이 방법은 생존함수 (Survival function)을 추정하기 위해 Hazard function에 대한 모델링을 합니다. 시간에 따라 변하는 Baseline hazard가 존재하고, Covariate 들의 effect에 따라 hazard가 달라진다고 모델링합니다. 그리고 hazard의 effect인 Hazard ratio가 시간에 따라 동일하다는 비례위험을 가정하여 Hazard function과 Survival function 을 추정하게 됩니다. 이 방법의 장점은 Survival function에 분포에 대해서 가정을하지 않기 때문에 강건하며, Covariate를 보정할 수 있다는 것입니다. 또한 Baseline hazard function이 어떤 함수인지를 알 필요가 없이 계수 추정을 통해 Hazard ratio만 구할 수 있다면 Covariate의 effect를 파악할 수 있습니다. 

모수적 모형

반면 모수적 모형은 Hazard function과 Survival function의 모양에 대해 Cox 보다 더 강한 '가정'을 합니다. 이 가정은 Hazard function과 Survival function이 어떤 분포를 따를 것이라는 믿음에서 출발합니다. 가장 많이 사용되는 모수적 모델은 Weibull 모형입니다. 

Weibull 모형의 Hazard function과 Survival function 은 다음과 같이 정의됩니다.

$$ h(t) = \lambda p t^{p-1} $$

$$ S(t) = exp(-\lambda t^p) $$ 

h(t)와 S(t)에 정의에 따라 h(t)가 저렇게 정의되면 S(t)가 저렇게 계산되는 것을 확인할 수 있습니다.

$$ \lambda = exp(\beta_0 + \beta_1 TRT) $$ 

만약 Treatment 변수에 대해 Weibull model을 만들 경우, 이 때 위와 같이 식을 놓고, beta0, beta1, p 3개의 계수를 추정하게 됩니다. 

또는 Hazard를 기준으로 계수를 추정하는 것이 아니라 Survival time을 기준으로 추정할 수 있는데, 이 경우가 모수적 방법에서 많이 사용되는 방법입니다. 이를 Acceleration Failure Time Model이라고 합니다. 

$$ S(t) = exp(-\lambda t^p) $$

$$ t = [-lnS(t)]^{1/p} \times \frac{1}{\lambda^{\frac{1}{p}}} $$

$$ \frac{1}{\lambda^{\frac{1}{p}}} = exp(\alpha_0 + \alpha_1 TRT) $$

AFT 모델에서는 위와 같이 놓고 alpha0, alpha1, p를 추정하게 됩니다. 

이렇게 두 가지 방법으로 구한 계수 beta와 alpha의 의미는 다른데, 만약 p가 고정되있다고 한다면, beta는 exponential을 취하면 Hazard ratio를 뜻하며, alpha는 exponential을 취하면 Accerleration factor가 됩니다. Acceleration factor에 대해서는 본 포스팅에서는 다루지 않을 예정이므로 이 자료를 참고 바랍니다. 

Strata 별로 p가 동일하다고 가정한 경우, Weibull model은 Proportional Hazard 가정과 Acceleration Failure Time 가정을 만족합니다. (https://en.wikipedia.org/wiki/Accelerated_failure_time_model)

PBC 데이터를 통한 Weibull model의 구축

pbc.dat

library(survival)

data <- read.table("pbc.dat")

# Column rename
colnames(data) <- c('id', 'futime', 'status', 'drug', 'age', 'sex', 'ascites', 'hepato', 'spiders', 'edema', 'bili', 'chol', 'albumin', 'copper', 'alk_phos', 'sgot', 'trig', 
                   'platelet', 'protime', 'stage')

head(data)

data <- data[data$drug != '.', ]
data$drug <- as.character(data$drug)
data$drug[(data$drug == 2)] = 0
data$drug <- factor(data$drug)

Drug = 0 플라시보, Drug = 1 페니실린에 대해 Death에 대한 Time to event 분석을 수행합니다.

Weibull Model Fitting

library(survival)
library(eha)

# AFT  
swei = survreg(s~drug, dist = "weibull",data = data) 

# PH  
sweiph = phreg(s~drug, dist = "weibull",data = data)

summary(swei) 
summary(sweiph)

Call:
survreg(formula = s ~ drug, data = data, dist = "weibull")
              Value Std. Error     z      p
(Intercept)  8.3167     0.1094 75.99 <2e-16
drug1       -0.0438     0.1431 -0.31  0.759
Log(scale)  -0.1532     0.0729 -2.10  0.036

Scale= 0.858 

Weibull distribution
Loglik(model)= -1348.2   Loglik(intercept only)= -1348.3
	Chisq= 0.09 on 1 degrees of freedom, p= 0.76 
Number of Newton-Raphson Iterations: 5 
n= 312 

Call:
phreg(formula = s ~ drug, data = data, dist = "weibull")

Covariate          W.mean      Coef Exp(Coef)  se(Coef)    Wald p
drug 
               0    0.491     0         1           (reference)
               1    0.509     0.051     1.052     0.167     0.759 

log(scale)                    8.317               0.109     0.000 
log(shape)                    0.153               0.073     0.036 

Events                    144 
Total time at risk        625985 
Max. log. likelihood      -1348.2 
LR test statistic         0.09 
Degrees of freedom        1 
Overall p-value           0.759387

Median Survival Time Comparison

predict(swei,newdata=data.frame(drug='0'),type=c("quantile"),se.fit=FALSE,p=0.5)
predict(swei,newdata=data.frame(drug='1'),type=c("quantile"),se.fit=FALSE,p=0.5)

Weibull Model을 통해 구한 Survival Function과 Kaplen-Meier estimation의 비교

data$event_all=ifelse(data$status==0,0,1)
data$futime = ifelse(data$futime==0,0.1,data$futime)
s = Surv(data$futime,data$event_all)

sKM <- survfit(s ~ drug, data=data)
plot(sKM)
lines(predict(swei,newdata=list(drug='0'),type="quantile", p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col='red')
lines(predict(swei,newdata=list(drug='1'),type="quantile", p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col='red')

위 그림에서 검은색 직선이 Kaplen-Meier Survival function estimation이며, 붉은색과 파란색이 Weibull 모델로 추정한 Survival function 입니다. Weibull 모델의 S(t)의 모양에 대한 가정이 맞다면 두 직선은 가까울 것입니다. 그리고 위 직선들은 p가 같다는 가정하에 구해진 직선들입니다. 즉, PH, AFT 가정을 한 후에 그려진 직선입니다. 그렇기 때문에 Weibull 모형의 Survival function에 대한 가정, 그리고, PH, AFT 가정을 해도 될지에 대해서 체크를 해보아야합니다.

Weibull 가정의 체크 

(https://stat.ethz.ch/education/semesters/ss2011/seminar/contents/handout_9.pdf)

Weibull 모형에서 가정을 체크하는 방법은 X 축을 log(t), Y 축을 log(-log(S(t)) 로 그래프를 그려보는 것입니다. 이 때, S(t)는 Kaplen-Meier 방법으로 구한 Survival function 입니다. 이 때, 그룹간의 그래프의 모양을 보고 Weibull 생존모형 가정이 맞는지, 또한 두 그룹에서 PH, AFT 가정이 맞는지를 체크할 수 있습니다. 

먼저, 두 그룹의 직선이 평행하게 나오면 Weibull 가정을 만족합니다. 그렇지 않으면 Weibull 가정이 만족되지 못하며 AFT 가정도 만족되지 못합니다. 또한 두 그룹의 기울기의 비가 일정하면 PH 가정을 만족합니다. 

sKM <- survfit(s~ drug, data=data)
kmfit.TRT1 = summary(sKM)
kmfit.TRT2 = data.frame(kmfit.TRT1$strata, kmfit.TRT1$time, kmfit.TRT1$surv)
names(kmfit.TRT2) = c("TRT","time","survival")
TRT0=kmfit.TRT2[kmfit.TRT2$TRT=="drug=0", ]
TRT1=kmfit.TRT2[kmfit.TRT2$TRT=="drug=1", ]
plot(log(TRT0$time), log(-log(TRT0$survival)), pch=16, xlab="time in days using log-arithmic scale", ylab="log-log survival",type="p", xlim=c(3,10), ylim=c(-6,1), font=2, lwd=2, font.lab=2, col="Black",main="log-log curves by treatment")
par(new=T)
plot(log(TRT1$time), log(-log(TRT1$survival)), pch=3, col="grey50", ylab="",xlab="",  xlim=c(3,10), ylim=c(-6,1), font=2, lwd=2, font.lab=2)

이 그래프의 경우 우선 그룹에서 거의 직선이라는 것을 확인할 수 있지만, 두 그래프가 평행하지는 않습니다. (교차가 일어납니다.) 그렇기 때문에 Weibull 가정을 만족하지만 PH 가정과 AFT 가정을 만족하지 않으며, 두 그룹에서 p가 다르다는 알 수 있습니다.