고급통계 - 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가 다르다는 알 수 있습니다.

Time-dependent variable의 정의 

Definition : 생존 분석에 있어서 시간에 따라 변수의 값이 바뀌는 변수를 time-dependent variable이라고한다. 

예 : 흡연 여부, 직장에서의 지위 등

Extended cox model : Time-dependent variable을 생존분석에서 사용하는 방법 

개념 : 만약 시간에 따라 변수의 값이 변화하는 변수가 있다면 extended cox model을 수행하여야한다. 

- 기존 데이터의 형태

- Extended cox model을 위한 데이터의 형태

모든 time에 대해 데이터를 쪼개서 start와 time으로 정의한다. 하나의 row를 시간에 따라 쪼개므로 한 사람당 여러개의 row가 생기게 된다. 만약 해당 시간 내에서 event가 발생한 경우, event는 1이된다.

library(data.table)
library(survival)
data <- read.csv("addict.csv")

- Extended cox model을 위한 데이터의 형태 만드는 법

data.cp <- survSplit(data, cut=data$time[data$event==1], end="time", event="event", start="start", id="ID2")

만약 아래와 같이 logdose라는 time-dependent variable이 있다고 하면 

data.cp$logdose <- data.cp$dose * log(data.cp$time)

Extended cox model 은 아래와 같이 수행할 수 있다.

Y = Surv(data.cp$start, data.cp$time, data.cp$event)
coxph(Y ~ clinic + logdose +prison+ cluster(ID2), data=data.cp)

Time-dependent variable의 활용 1 : PH (proportional hazard) Assumption check

개념 : PH Assumption에 위반이 되는 것이 의심되는 변수가 있을 때, 이 변수에 대한 time-dependent variable을 정의함으로서 PH 가정을 체크할 수 있다. 

만약에 위와 같은 데이터에 대해서 dose에 대해 PH Assumption 가정이 위반된다고 해보자. 즉, 이말은 시간에 따라 dose의 hazard ratio가 달라진다는 것이다. 이를 테스트하기위해 dose*g(t)라는 새로운 변수를 정의하고 이 변수의 계수가 유의한지 확인하는 방법을 통해 해당 변수에 대한 PH Assumption을 체크할 수가 있다. 

Y = Surv(data.cp$start, data.cp$time, data.cp$event)
coxph(Y ~ clinic + dose + logdose +prison+ cluster(ID2), data=data.cp)

Call:
coxph(formula = Y ~ clinic + dose + logdose + prison + cluster(ID2), 
    data = data.cp)

            coef exp(coef) se(coef) robust se     z       p
clinic  -1.01988   0.36064  0.21542   0.23637 -4.31 1.6e-05
dose    -0.08262   0.92070  0.03598   0.02960 -2.79  0.0053
logdose  0.00862   1.00865  0.00645   0.00525  1.64  0.1007
prison   0.34063   1.40584  0.16747   0.15972  2.13  0.0329

Likelihood ratio test=66.35  on 4 df, p=1e-13
n= 18708, number of events= 150 
   (1379 observations deleted due to missingness)

위 결과에서 logdose의 p-value > 0.05이므로, dose가 PH Assumption을 만족한다고 결론내릴 수 있다. 위 방법은 Wald test를 한 것이며, 다른 방법으로는 Likelihood ratio test를 통해 p-value를 구할 수도 있다. 

Likelihood ratio test

Y = Surv(data.cp$start, data.cp$time, data.cp$event)
mod1 <- coxph(Y ~ clinic + dose + logdose +prison+ cluster(ID2), data=data.cp)
mod2 <- coxph(Y ~ clinic + dose +prison+ cluster(ID2), data=data.cp)

print(mod1$loglik[2])
print(mod2$loglik[2])

LRT <- (-2) * (mod2$loglik[2]-mod1$loglik[2])
pvalue <- 1 - pchisq(LRT, 2) 

이는 logdose를 포함한 full model과 logdose를 포함하지 않은 restricted model의 likelihood ratio를 구하고 이것이 카이제곱 분포를 따른다는 것을 이용한 검정이다. p-value를 계산한 결과 0.4 정도로 위의 wald test와 마찬가지로 logdose는 의 계수는 유의하지 않았다. 

따라서 dose에 대해 PH Assumption이 어긋나지 않는다고 결론 내릴 수 있다.

Time-dependent variable의 활용 2 : PH Assumption 만족 안하는 경우, 시간에 따라 달라지는 Hazard ratio 모델링

개념 : 어떤 변수에 대하여 heavy side function을 이용한 time-dependent variable을 정의함으로서 시간에 따라 달라지는 hazard ratio를 모델링할 수 있다. 

heavy side function은 어떤 값을 기준으로 그 이상일 때는 1, 미만일 때는 0을 돌려주는 함수를 말한다. 

예제 : 만약 clinic에 대해 시간에 따라 달라지는 HR을 모델링 하고 싶다고 해보자. 

fit <- survfit(Surv(data$time, data$event)~data$clinic)
plot(fit)

이것이 clinic=1과 clinic=2일 때 각각 그린 카플란 마이어 생존함수 추정값이다. 위 그림에서 hazard ratio가 시간에 따라 바뀐다는 것이 의심스러울 수 있다. 이 때, 어떤 시점을 기준으로 heavy side function을 이용한 time-dependent variable을 정의함으로써 이를 모델링할 수 있다. 

365를 기준으로 데이터를 쪼갠다. 

data.cp365 <- survSplit(data, cut=365, end="time", event="event", start="start", id="ID2")

data.cp365$hv1 <- data.cp365$clinic * (data.cp365$start < 365)
data.cp365$hv2 <- data.cp365$clinic * (data.cp365$start >= 365)

위와 같이 두 개의 heavy side function을 원래 clinic의 값에 곱해준 것을 두 개의 변수로 정의한다.  

Extended cox model

Y365 <- Surv(data.cp365$start, data.cp365$time, data.cp365$event)
coxph(Y365 ~ dose + hv2 + hv1 + prison + cluster(ID2), data=data.cp365)

Call:
coxph(formula = Y365 ~ dose + hv2 + hv1 + prison + cluster(ID2), 
    data = data.cp365)

           coef exp(coef) se(coef) robust se     z       p
dose   -0.03548   0.96514  0.00643   0.00652 -5.44 5.3e-08
hv2    -1.83052   0.16033  0.38595   0.39838 -4.59 4.3e-06
hv1    -0.45937   0.63168  0.25529   0.25998 -1.77   0.077
prison  0.37795   1.45929  0.16841   0.16765  2.25   0.024

Likelihood ratio test=74.25  on 4 df, p=3e-15
n= 360, number of events= 150 

이 때 결과를 보면 hv1의 HR은 0.63, hv2의 HR은 0.16이다. 이 말은 clinic 1과 clinic2의 hazard ratio가 365미만에서는 0.63이었는데, 365 이상에서는 0.16이라고 해석할 수 있다. 즉, 두 clinic에서 보이는 생존률에 대한 경향성이 시간이 지날 수록 강해지고 있는 것이다. 이는 위의 카플란 마이서 생존곡선에서도 확인할 수 있다. 

Cox 비례위험 모형의 가정 체크

개념 : Cox 비례위험 모형은 "비례위험 가정" 을한다. 

비례위험가정은 시간에 따라 어떤 변수에서의 Hazard ratio가 일정해야한다는 것이다.

비례위험가정이 어긋나면, Stratified cox regression 또는 time-dependant variable 정의를 통한 extended cox regression을 수행할 수 있다. 

비례위험 가정의 체크는 아래방식으로 할 수 있다. 

1. Log-log plot

2. Observed-expected plot

3. Goodness of fit test

본 포스팅에서는 R의 survival, survminer 패키지를 통한 실습을 통해 cox 비례 위험 모형의 가정을 체크하는 방법을 알아본다. 

addict.csv

library(data.table)
library(survival)
data <- read.csv("addict.csv")

head(data)

먼저 데이터의 형태는 위와 같이 생겼다. 

관심변수 : clinic 에 대하여 Kaplan-Meier 플롯 그리기

# Treatment 별로 그룹지어서 KM Estimation
fit = survfit(y~clinic, data=data)
par(mai=c(1,1,1,1))
plot(fit,main="KM Curve", xlab="Time(Week)", ylab=expression(paste(hat(S),"(t)")), col=c("blue", "red"))

우리는 이것이 비례위험 가정을 만족하는지를 확인하고 싶다. 비례위험 가정을 만족하려면 모든 timepoints에서 Hazard ratio가 같아야한다. hazard란 f(t)/S(t) (참고 : https://3months.tistory.com/349?category=743476) 로 위 그림만 봐서는 가정을 만족하는지 한 눈에 알 수 없다.  

로그로그 플롯

개념 : 로그로그 플롯을 그리고 그룹간에 평행하게 나오면 비례위험가정을 만족한다. 

plot(fit, fun="cloglog", lty=1:2, col=c("Black", "Grey50"), lwd=2, font=2, ylim=c(-3,2), font.lab=2, main="log-log KM curves by Rx", ylab="log-log survival", xlab="time in weeks using logarithmic scale")
legend("bottomright",lty=1:2,legend=c("Group1","Group2"),bty="n", text.font=2, lwd=2, col=c("Black", "Grey50"))

관찰값-기댓갓 플롯

개념 : Kaplen-Meier 생존 함수와 Cox regression을 통해 구한 생존함수를 동시에 그려본다.

이 때, Expected plot은 Cox 비례위험 모형에서 나오는 생존함수로 관찰한 생존함수와 모양의 차이가 많이나면 비례위험 가정을 만족하지 않는 것이다. 

kmfit <- survfit(y~clinic, data=data)
#observed
plot(kmfit, lty="dashed", col=c("Black", "Grey50"), lwd=2, font=2, font.lab=2, main="Observed Versus Expected Plots by Clinic", ylab="Survival probability", xlab="Survival time")
par(new=T)
#expected
exp <- coxph(y~clinic, data=data)
new_df <- data.frame(clinic=c(1,2))
kmfit.exp <- survfit(exp, newdata=new_df)
plot(kmfit.exp, lty="solid", col=c("Black", "Grey50"), lwd=2, font=2, font.lab=2)

Goodness of fit test

개념 : Goodness of fit test(GoF Test) 는 관찰값-기댓값의 차이를 통계적으로 검정하는 방법이다. (카이스퀘어 분포 이용)

# Gooodness of fit test
mod1 <- coxph(y~clinic, data=data)
cox.zph(mod1, transform=rank)

plot(cox.zph(mod1, transform=rank), se=F, var="clinic")
abline(h=0, lty=3)
plot(cox.zph(mod1, transform=rank), se=F, var="prison")
abline(h=0, lty=3)

        rho chisq       p
clinic -0.255  10.5 0.00118

p-value가 0.001이므로 clinic 변수는 비례위험 가정을 만족하지 않는다. 이는 그림을 통해서도 확인할 수 있다. 위 그림에서 직선이 평행하지를 확인하면되는데, 직선은 추정된 beta값을 의미하는데 직선이 평행하지 않고 어떠한 경향성이 보이면 beta 값이 시간에 따라 바뀐다는 의미로 비례위험 가정에 어긋난다. 

참고 : https://www.r-bloggers.com/cox-model-assumptions/