생존분석 실습

Kaplan-meier curve, 비례위험가정 확인, Time-dependent analysis 그리고 모수적 생존분석을 중심으로 R 코드를 정리했습니다. 본 내용은 성균관의대 사회의학교실 특강에서 실습할 예정입니다.
presentation
R
shiny
docker
Author
Affiliation

Zarathu Co.,Ltd

Published

October 31, 2020

김진섭 대표는 성균관의대 사회의학교실 김종헌 교수님 수업에 참가, Kaplan-meier curve, 비례위험가정 및 모형적합도, Time-dependent covariate 그리고 모수적 생존분석을 중심으로 R 코드를 실습할 예정입니다.

요약

  • 자체 개발한 jskm 패키지로 kaplan-meier 그림을 그린다.

  • 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 으로 구간별 생존율을 구하자.

library(DT);library(survival);library(jskm)
datatable(veteran, rownames = F, caption = "Example data", options = list(scrollX = T))
sfit <- survfit(Surv(time, status) ~ trt, data = veteran)
summary(sfit, times = c(100, 200, 300, 365), extend = T)
Call: survfit(formula = Surv(time, status) ~ trt, data = veteran)

                trt=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
  100     34      34   0.5020  0.0606       0.3962        0.636
  200     12      19   0.1947  0.0501       0.1176        0.322
  300      5       6   0.0885  0.0371       0.0390        0.201
  365      4       1   0.0708  0.0336       0.0279        0.180

                trt=2 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
  100     21      45    0.333  0.0578       0.2367        0.467
  200     13       7    0.216  0.0517       0.1354        0.345
  300      8       4    0.146  0.0454       0.0797        0.269
  365      6       2    0.110  0.0407       0.0530        0.227

trt 1 은 “Standard”, 2 는 “Test” 이며 jskm 을 적용하면 아래와 같다.

jskm(sfit)

라벨을 수정하고, risk table 과 log-rank p-value 를 추가하자.

jskm(sfit, ystrataname = "Treat", ystratalabs = c("Standard", "Test"), table = T, pval = T)

십자가 무늬는 실제 censoring 이 발생한 부분이며 mark = F 로 숨길 수 있다. 생존율이 아닌 누적발생률을 %로 보는 코드는 아래와 같다.

jskm(sfit, ystrataname = "Treat", ystratalabs = c("Standard", "Test"), table = T, pval = T, 
     marks = F, cumhaz = T, surv.scale = "percent" )

p-value 위치는 pval.coord legend 위치는 legendposition 옵션을 이용한다. 선을 흑백으로 바꾸려면 linecols = "black" 을 추가한다. legendposition 은 x,y 값 모두 0~1 scale 임을 주의하자.

jskm(sfit, ystrataname = "Treat", ystratalabs = c("Standard", "Test"), table = T, pval = T, 
     marks = F, pval.coord = c(100, 0.1), legendposition = c(0.85, 0.6), linecols = "black")

마지막으로 특정 시간을 기준으로 나누어보는 landmark analysis 옵션을 소개한다.

jskm(sfit, ystrataname = "Treat", ystratalabs = c("Standard", "Test"), table = T, pval = T, 
     marks = F, cut.landmark = 365)
Warning: Removed 2 row(s) containing missing values (geom_path).

연속변수의 최적 cut-off 구하기

maxstat 패키지를 이용한다.

library(maxstat)
mtest <- maxstat.test(Surv(time, status) ~ karno, data = veteran, smethod = "LogRank")
mtest

Maximally selected LogRank statistics using none

data:  Surv(time, status) by karno
M = 4.6181, p-value = NA
sample estimates:
estimated cutpoint 
                40 
cut <- mtest$estimate
veteran$karno_cat <- factor(as.integer(veteran$karno >= cut))

sfit2 <- survfit(Surv(time, status) ~ karno_cat, data = veteran)
jskm(sfit2, ystrataname = "Karno", ystratalabs = paste(c("<", "≥"), cut), table = T, pval = T)

비례위험가정 확인

Logrank test, Cox model 로 추정할 때 비례위험을 가정하므로 이것이 깨지면 큰일이다. 본 글에서는 비례위험가정을 확인하는 그림 2개와 테스트를 소개한다. 자세한 내용은 https://3months.tistory.com/357?category=743476 를 참고하기 바란다.

Log-log plot

\(\log(t)\)\(\log(-\log(S(t)))\) 관계를 그림으로 보는 방법이다. 왜 로그를 이용하는지는 모수적 생존분석에서 이야기하겠다.

plot(sfit, fun="cloglog", lty=1:2, col=c("Black", "Grey50"), lwd=2, font.lab=2, main="Log-log KM curves by Treat", 
     ylab="log-log survival", xlab="Time (log scale)")
legend("bottomright",lty=1:2,legend=c("Standard", "Test"), bty="n", lwd=2, col=c("Black", "Grey50"))

두 선이 평행한지 확인하면 되고 직선인지 곡선인지는 상관없다. 모수적 생존분석에서 다룰 weibull 모형에서는 직선인지도 확인해야 한다.

Observed-expected plot

비례위험을 가정하는 cox model 예상과 비교하는 방법이다.

plot(sfit, lty="dashed", col=c("Black", "Grey50"), lwd=2, font=2, font.lab=2, main="Observed Versus Expected Plots by Treat", 
     ylab="Survival probability", xlab="Time")
par(new = T)

#expected
exp <- coxph(Surv(time, status) ~ trt, data = veteran)
new_df <- data.frame(trt = c(1, 2))
kmfit.exp <- survfit(exp, newdata = new_df)
plot(kmfit.exp, lty = "solid", col=c("Blue", "Red"), lwd=2, font.lab=2)

Goodness of fit

cox.zph 함수로 통계검정을 수행한다.

cox.zph(exp)
       chisq df    p
trt     3.54  1 0.06
GLOBAL  3.54  1 0.06
plot(cox.zph(exp), var = "trt")
abline(h = 0, lty = 3)

선이 시간 상관없이 일정할수록, 즉 x축과 평행할수록 비례위험가정을 만족한다고 판단한다. 위 그림은 x축과 평행은 아니지만 경향성이 있다고 볼수도 없는 애매한 느낌이며 p 는 0.06 이다.

모형 비교

Cox 모형에서 얻은 log-likelihood 값으로 여러 모형을 비교할 수 있다. 모형들은 n수가 전부 동일 해야 비교 가능하므로, 에러 나올땐 먼저 결측치를 확인하자.

exp$loglik
[1] -505.4491 -505.4442
exp2 <- coxph(Surv(time, status) ~ trt + age, data = veteran)
exp3 <- coxph(Surv(time, status) ~ trt + age + celltype, data = veteran)

anova(exp, exp2, exp3)
Analysis of Deviance Table
 Cox model: response is  Surv(time, status)
 Model 1: ~ trt
 Model 2: ~ trt + age
 Model 3: ~ trt + age + celltype
   loglik   Chisq Df P(>|Chi|)    
1 -505.44                         
2 -505.14  0.6162  1    0.4325    
3 -492.43 25.4161  3 1.264e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

step 함수를 이용, AIC 기반 최적 모형을 고를 수 있다. scope 옵션으로 빠지면 안 될 변수를 미리 정한다.

step(exp3, scope = list(lower = ~ 1))
Start:  AIC=994.86
Surv(time, status) ~ trt + age + celltype

           Df     AIC
- age       1  993.04
- trt       1  993.65
<none>         994.86
- celltype  3 1014.27

Step:  AIC=993.04
Surv(time, status) ~ trt + celltype

           Df     AIC
- trt       1  992.05
<none>         993.04
- celltype  3 1012.89

Step:  AIC=992.05
Surv(time, status) ~ celltype

           Df     AIC
<none>         992.05
- celltype  3 1010.90
Call:
coxph(formula = Surv(time, status) ~ celltype, data = veteran)

                    coef exp(coef) se(coef)     z        p
celltypesmallcell 1.0013    2.7217   0.2535 3.950 7.83e-05
celltypeadeno     1.1477    3.1510   0.2929 3.919 8.90e-05
celltypelarge     0.2301    1.2588   0.2773 0.830    0.407

Likelihood ratio test=24.85  on 3 df, p=1.661e-05
n= 137, number of events= 128 

Time-dependent analysis

자세한 내용은 https://cran.r-project.org/web/packages/survival/vignettes/timedep.pdf 를 참고하기 바란다.

비례위험가정 깨졌을 때 (time-dependent coefficients)

어떤 공변량이 비례위험가정을 만족하지 않을 경우, 먼저 survSplit 으로 time 을 쪼개 몇 개의 그룹으로 나눈다.

vet2 <- survSplit(Surv(time, status) ~ ., data = veteran, cut=c(90, 180), episode = "tgroup", id = "id")
datatable(vet2, rownames = F, caption = "Time split data", options = list(scrollX = T))

이제 공변량의 계수를 시간그룹 별로 따로 구한다.

vfit2 <- coxph(Surv(tstart, time, status) ~ trt + prior + karno:strata(tgroup), data=vet2)
summary(vfit2)
Call:
coxph(formula = Surv(tstart, time, status) ~ trt + prior + karno:strata(tgroup), 
    data = vet2)

  n= 225, number of events= 128 

                                  coef exp(coef)  se(coef)      z Pr(>|z|)    
trt                          -0.011025  0.989035  0.189062 -0.058    0.953    
prior                        -0.006107  0.993912  0.020355 -0.300    0.764    
karno:strata(tgroup)tgroup=1 -0.048755  0.952414  0.006222 -7.836 4.64e-15 ***
karno:strata(tgroup)tgroup=2  0.008050  1.008083  0.012823  0.628    0.530    
karno:strata(tgroup)tgroup=3 -0.008349  0.991686  0.014620 -0.571    0.568    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                             exp(coef) exp(-coef) lower .95 upper .95
trt                             0.9890      1.011    0.6828    1.4327
prior                           0.9939      1.006    0.9550    1.0344
karno:strata(tgroup)tgroup=1    0.9524      1.050    0.9409    0.9641
karno:strata(tgroup)tgroup=2    1.0081      0.992    0.9831    1.0337
karno:strata(tgroup)tgroup=3    0.9917      1.008    0.9637    1.0205

Concordance= 0.725  (se = 0.024 )
Likelihood ratio test= 63.04  on 5 df,   p=3e-12
Wald test            = 63.7  on 5 df,   p=2e-12
Score (logrank) test = 71.33  on 5 df,   p=5e-14

반복측정 공변량이 있을 때

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6015946/pdf/atm-06-07-121.pdf 예제를 이용하였다.

library(survsim)
Loading required package: eha
Loading required package: statmod
N=100 #number of patients
set.seed(123)
df.tf<-simple.surv.sim(#baseline time fixed
 n=N, foltime=500,
 dist.ev=c('llogistic'),
 anc.ev=c(0.68), beta0.ev=c(5.8),
 anc.cens=1.2,
 beta0.cens=7.4,
 z=list(c("unif", 0.8, 1.2)),
 beta=list(c(-0.4),c(0)),
 x=list(c("bern", 0.5),
 c("normal", 70, 13)))

for (v in 4:7){
  df.tf[[v]] <- round(df.tf[[v]])
}

names(df.tf)[c(1,4,6,7)]<-c("id", "time", "grp","age")
df.tf <- df.tf[, -3]

datatable(df.tf, rownames = F, caption = "df.tf: Original data", options = list(scrollX = T))
 nft<-sample(1:10,
 N,replace=T)#number of follow up time points
crp<-round(abs(rnorm(sum(nft)+N,
 mean=100,sd=40)),1)
time<-NA
id<-NA
i=0
for(n in nft){
i=i+1
time.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 변수를 만들자.

df <- tmerge(df.tf, df.tf, id = id, status1 = event(time, status))

datatable(df, rownames = F, caption = "df: add tstart/tstop", options = list(scrollX = T))

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 를 참고하기 바란다.

df2 <- tmerge(df, df.td, id = id, crp = tdc(time, crp))

datatable(df2, rownames = F, caption = "df2: final", options = list(scrollX = T))

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)\) 가 정비례관계이기 때문이다.

\[ \begin{align} S(t) &= \exp(-\lambda t^p) \\ -\log(S(t)) &= \lambda t^p \\ \log(-\log(S(t))) &= \log(\lambda) + p\log(t) \\ \log(-\log(S(t))) &\propto \log(t) \end{align} \]

\(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)

\(S(t)\) 를 구할 수 없는 cox model 의 그림과 비교해보자.

model.cox <- exp
kmfit.exp <- survfit(exp, newdata = data.frame(trt = 1))
plot(kmfit.exp, lty = c(1, 2, 2), col=1, lwd=2, xlab = "Days", ylab = "Survival")

지금까지 생존분석 때 고려할 내용을 다루었으며 처음의 요약을 반복하면 아래와 같다.

  • 자체 개발한 jskm 패키지로 kaplan-meier 그림을 그린다.

  • Log-log plot, Observed-expected plot 으로 비례위험가정을 확인 후, cox.zph 함수로 p-value 를 구한다.

  • anova 로 여러 모형의 log-likelohood 를 비교하고, step 으로 AIC 기반 최적모형을 고를 수 있다.

  • Time-dependent analysis 는 (1) 비례위험가정이 깨졌을 때, (2) 반복측정 공변량이 있을 때 수행한다.

  • 모수적 생존분석은 생존함수 \(S(t)\) 를 구할 수 있어 예측모형을 만들 수 있다.

자세한 내용은 중간중간 링크한 자료들을 참고하기 바란다.

Citation

BibTeX citation:
@online{kim2020,
  author = {Jinseob Kim},
  title = {생존분석 {실습}},
  date = {10/31/2020},
  url = {https://blog.zarathu.com/posts/2020-10-29-survivalpractice},
  langid = {en}
}
For attribution, please cite this work as:
Jinseob Kim. 10AD–31AD. “생존분석 실습.” 10AD–31AD. https://blog.zarathu.com/posts/2020-10-29-survivalpractice.