Log-log plot, Observed-expected plot 으로 비례위험가정을 확인 후, cox.zph 함수로 p-value 를 구한다.
anova 로 여러 모형의 log-likelohood 를 비교하고, step 으로 AIC 기반 최적모형을 고를 수 있다.
Time-dependent analysis 는 (1) 비례위험가정이 깨졌을 때, (2) 반복측정 공변량이 있을 때 수행한다.
모수적 생존분석은 생존함수 \(S(t)\) 를 구할 수 있어 예측모형을 만들 수 있다.
Kaplan-meier plot
Kaplan-meier plot 은 R 기본 plot에서도 제공하지만, survminer 패키지의 ggsurvplot 함수에서 다양한 옵션을 제공한다. 본 실습에서는 본사가 개발한 jskm 패키지의 jskm 함수를 survival 패키지 내장 데이터 veteran 에 적용하겠다. 우선 패키지를 불러온 후 survfit 으로 구간별 생존율을 구하자.
nft<-sample(1:10,N,replace=T)#number of follow up time pointscrp<-round(abs(rnorm(sum(nft)+N, mean=100,sd=40)),1)time<-NAid<-NAi=0for(ninnft){i=i+1time.n<-sample(1:500,n)time.n<-c(0,sort(time.n))time<-c(time,time.n)id.n<-rep(i,n+1)id<-c(id,id.n)}df.td<-cbind(data.frame(id,time)[-1,],crp)datatable(df.td, rownames =F, caption ="df.td: Time dependent CRP", options =list(scrollX =T))
df.tf 는 기본정보가 담긴 데이터, df.td 는 time-dependent covariate 가 담긴 데이터이다. tmerge 함수를 2번 실행하면 두 정보를 합칠 수 있다. 먼저 df.tf 만 이용해서 tstart, tstop 변수를 만들자.
tmerge 함수의 첫번째는 baseline data, 둘째는 time-dependent covariate 가 담긴 데이터가 들어가지만, tstart, tstop 를 만들기 위해 모두 df.tf 를 넣었다. status1 이라는 변수를 event(time, status) 로 지정함으로서 tstart, tstop 을 인식할 수 있다. status1 변수 자체는 status 와 동일하다. 이렇게 만든 df 에 time-dependent 정보가 담긴 df.td 를 결합하면 원하는 데이터를 얻을 수 있다. tmerge 의 자세한 내용은 https://ww2.amstat.org/meetings/sdss/2018/onlineprogram/ViewPresentation.cfm?file=304494.pdf 를 참고하기 바란다.
crp 변수를 tdc(time, crp) 로 만들었다. 이제 cox model 을 실행할 수 있는데, 반복측정정보를 cluster 옵션에 넣는 것을 잊지 말자.
model.td<-coxph(Surv(tstart, tstop, status1)~grp+age+crp, data =df2, cluster =id)summary(model.td)
Call:
coxph(formula = Surv(tstart, tstop, status1) ~ grp + age + crp,
data = df2, cluster = id)
n= 376, number of events= 67
coef exp(coef) se(coef) robust se z Pr(>|z|)
grp 0.5022750 1.6524764 0.2525914 0.2555150 1.966 0.0493 *
age 0.0005535 1.0005536 0.0081077 0.0072342 0.077 0.9390
crp 0.0007922 1.0007925 0.0027391 0.0023373 0.339 0.7347
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
grp 1.652 0.6052 1.0015 2.727
age 1.001 0.9994 0.9865 1.015
crp 1.001 0.9992 0.9962 1.005
Concordance= 0.554 (se = 0.04 )
Likelihood ratio test= 4.21 on 3 df, p=0.2
Wald test = 4.34 on 3 df, p=0.2
Score (logrank) test = 4.18 on 3 df, p=0.2, Robust = 4.55 p=0.2
(Note: the likelihood ratio and score tests assume independence of
observations within a cluster, the Wald and robust score tests do not).
모수적(parametric) 생존분석
Cox model 은 baseline hazard 없이도 HR 을 구할 수 있는 장점이 있다. 아래 식
\[h(t) = h_0(t) \cdot \exp(\sum \beta_i x_i)\]
에서 \(h_0(t)\) 를 몰라도 \(\beta\) 들을 구할 수 있다는 뜻이고, cox model 이 준모수적(semi-parametric) 모형으로 불리는 이유이기도 하다. 그러나 Cox model 로 예측모형을 만들 때 이것은 단점이 된다. \(t\) 년 생존율을 구할 수 없기 때문이다. 생존함수 \(S(t)\) 는 아래처럼 계산하는데
\[S(t) = \int_{0}^{t} h(u) \,du\] baseline hazard 를 모르므로 \(h(t)\) 도 알 수 없고 따라서 \(S(t)\) 도 수식으로 표현할 수 없다. Cox model 로 예측모형을 만든 연구는 (1) 데이터에서 시간 \(t\) 마다 \(S(t)\) 의 값을 직접 구해 이용하거나, (2) 인구집단통계에서 \(S(t)\) 를 얻어온다.
그러면 baseline hazard 가 어떤 형태라고 가정하면 어떨까? 이것이 모수적 생존분석이며 cox model 과 장단점을 비교하면 아래와 같다.
Cox model
– distribution of survival time unkonwn
– Less consistent with theoretical \(S(t)\) (typically step function)
+ Does not rely on distributional assumptions
+ Baseline hazard not necessary for estimation of hazard ratio
Parametric Survival Model
+ Completely specified \(h(t)\) and \(S(t)\)
+ More consistent with theoretical \(S(t)\)
+ time-quantile prediction possible
– Assumption on underlying distribution
아래는 대표적인 분포들이며 본 글에서는 흔히 쓰는 weibull 을 다루려 한다.
아까 비례위험가정 얘기할 때 weibull 모형은 log-log 그래프가 직선인지도 확인해야 한다고 했는데, 그 이유는 아래 식에 나와있듯이 \(\log(-\log(S(t)))\) 와 \(\log(t)\) 가 정비례관계이기 때문이다.
\(p\) 를 scale parameter 라 하며 \(p = 1\) 이면 baseline hazard 가 시간에 따라 일정함을 의미하며, 자세한 내용은 https://stat.ethz.ch/education/semesters/ss2011/seminar/contents/handout_9.pdf 를 참고하자. R의 survreg 함수를 이용하며, 결과해석은 cox model 과 동일한데 scale parameter 값이 추가로 나온다(scale parameter를 미리 정할 수도 있다).
model.weibull<-survreg(Surv(time, status)~trt, data =veteran)summary(model.weibull)
Call:
survreg(formula = Surv(time, status) ~ trt, data = veteran)
Value Std. Error z p
(Intercept) 4.7218 0.3275 14.42 <2e-16
trt 0.0478 0.2079 0.23 0.818
Log(scale) 0.1585 0.0673 2.35 0.019
Scale= 1.17
Weibull distribution
Loglik(model)= -748.1 Loglik(intercept only)= -748.1
Chisq= 0.05 on 1 degrees of freedom, p= 0.82
Number of Newton-Raphson Iterations: 5
n= 137
Scale = 1.17 임을 확인할 수 있고, trt 그룹별 \(S(t)\) 를 그려보면 아래와 같다.
pcut<-seq(0.01, 1, by =0.01)## 1%-99%ptime<-predict(model.weibull, newdata =data.frame(trt =1), type ="quantile", p =pcut, se =T)matplot(cbind(ptime$fit, ptime$fit+1.96*ptime$se.fit, ptime$fit-1.96*ptime$se.fit), 1-pcut, xlab ="Days", ylab ="Survival", type ='l', lty =c(1, 2, 2), col=1)
