Time-dependent variable의 정의
Definition : 생존 분석에 있어서 시간에 따라 변수의 값이 바뀌는 변수를 time-dependent variable이라고한다.
예 : 흡연 여부, 직장에서의 지위 등
Extended cox model : Time-dependent variable을 생존분석에서 사용하는 방법
개념 : 만약 시간에 따라 변수의 값이 변화하는 변수가 있다면 extended cox model을 수행하여야한다.
- 기존 데이터의 형태
ID | clinic | event | time | prison | dose |
---|---|---|---|---|---|
1 | 1 | 1 | 428 | 0 | 50 |
2 | 1 | 1 | 275 | 1 | 55 |
3 | 1 | 1 | 262 | 0 | 55 |
4 | 1 | 1 | 183 | 0 | 30 |
5 | 1 | 1 | 259 | 1 | 65 |
6 | 1 | 1 | 714 | 0 | 55 |
- Extended cox model을 위한 데이터의 형태
ID | clinic | prison | dose | ID2 | start | time | event |
---|---|---|---|---|---|---|---|
1 | 1 | 0 | 50 | 1 | 0 | 7 | 0 |
1 | 1 | 0 | 50 | 1 | 7 | 13 | 0 |
1 | 1 | 0 | 50 | 1 | 13 | 17 | 0 |
1 | 1 | 0 | 50 | 1 | 17 | 19 | 0 |
1 | 1 | 0 | 50 | 1 | 19 | 26 | 0 |
1 | 1 | 0 | 50 | 1 | 26 | 29 | 0 |
모든 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 가정을 체크할 수 있다.
ID | clinic | event | time | prison | dose |
---|---|---|---|---|---|
1 | 1 | 1 | 428 | 0 | 50 |
2 | 1 | 1 | 275 | 1 | 55 |
3 | 1 | 1 | 262 | 0 | 55 |
4 | 1 | 1 | 183 | 0 | 30 |
5 | 1 | 1 | 259 | 1 | 65 |
6 | 1 | 1 | 714 | 0 | 55 |
만약에 위와 같은 데이터에 대해서 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")
ID | clinic | prison | dose | ID2 | start | time | event |
---|---|---|---|---|---|---|---|
1 | 1 | 0 | 50 | 1 | 0 | 365 | 0 |
1 | 1 | 0 | 50 | 1 | 365 | 428 | 1 |
2 | 1 | 1 | 55 | 2 | 0 | 275 | 1 |
3 | 1 | 0 | 55 | 3 | 0 | 262 | 1 |
4 | 1 | 0 | 30 | 4 | 0 | 183 | 1 |
5 | 1 | 1 | 65 | 5 | 0 | 259 | 1 |
data.cp365$hv1 <- data.cp365$clinic * (data.cp365$start < 365) data.cp365$hv2 <- data.cp365$clinic * (data.cp365$start >= 365)
위와 같이 두 개의 heavy side function을 원래 clinic의 값에 곱해준 것을 두 개의 변수로 정의한다.
ID | clinic | prison | dose | ID2 | start | time | event | hv1 | hv2 |
---|---|---|---|---|---|---|---|---|---|
1 | 1 | 0 | 50 | 1 | 0 | 365 | 0 | 1 | 0 |
1 | 1 | 0 | 50 | 1 | 365 | 428 | 1 | 0 | 1 |
2 | 1 | 1 | 55 | 2 | 0 | 275 | 1 | 1 | 0 |
3 | 1 | 0 | 55 | 3 | 0 | 262 | 1 | 1 | 0 |
4 | 1 | 0 | 30 | 4 | 0 | 183 | 1 | 1 | 0 |
5 | 1 | 1 | 65 | 5 | 0 | 259 | 1 | 1 | 0 |
6 | 1 | 0 | 55 | 6 | 0 | 365 | 0 | 1 | 0 |
6 | 1 | 0 | 55 | 6 | 365 | 714 | 1 | 0 | 1 |
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에서 보이는 생존률에 대한 경향성이 시간이 지날 수록 강해지고 있는 것이다. 이는 위의 카플란 마이서 생존곡선에서도 확인할 수 있다.
'Data science > Statistics' 카테고리의 다른 글
고급통계 - Net reclassification improvements (0) | 2019.01.23 |
---|---|
생존분석 - 모수적 모형 (Weibull model) (0) | 2018.12.05 |
생존분석 - Cox 비례위험 모형의 가정 체크 (1) | 2018.11.29 |
생존분석 - Stratified Cox 비례위험모형 (0) | 2018.11.28 |
생존분석 - Competing Risk를 고려한 생존 분석 (3) | 2018.11.25 |