생존분석 - Time-dependent variable의 활용

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에서 보이는 생존률에 대한 경향성이 시간이 지날 수록 강해지고 있는 것이다. 이는 위의 카플란 마이서 생존곡선에서도 확인할 수 있다.