Survival Data Analysis
Seo Young Park
2024-02-05
Survival Data Analysis란 time to event 타입의 변수를 결과변수로 하는 분석을 말합니다. time to event 타입의 변수란, (1) 이벤트 여부와 (2) 이벤트까지 걸린 시간 (또는 이벤트가 없었다는 것을 마지막으로 확인한 시간) 두가지로 구성된 변수를 말합니다. 임상연구에서는 보통 사망이나 재발과 같은 이벤트의 위험에 대해 분석할 때 쓰게 되죠. 그러면 R로 Survival Data Analysis를 하는 방법을 살펴보겠습니다.
\(~\)
Kaplan-Meier Curve 그리기
먼저, working directory를 설정하고 load() 함수를
이용해서 Session 1, 2에서 정리해둔 데이터를 불러옵니다.
setwd("C:\\Users\\biostat81\\Dropbox\\Teaching\\AMC Radiology R 2024")
load("clean_data.RData")오늘 데이터 분석에 필요한 다섯가지 패키지를 불러옵니다. 이 중 아직
설치하지 않은 패키지가 있다면 install.packages() 함수를
이용하여 먼저 패키지들을 설치한 후에 아래 코드를 실행하세요.
library(survival)
library(survminer)
library(broom)
library(lmtest)
library(timeROC)R에서 생존자료를 다루려면 survival 이라는 패키지를
이용해야합니다. 이 패키지의 Surv 함수를 이용하여 time to
event 타입의 변수를 인식시킬 수 있습니다. 이 강의에서 사용하고 있는 예제
데이터를 불러와서, time to event 타입의 변수를 만들어봅시다.
Surv 함수 안에 time과 event
두가지 값을 넣어줘야 하는데, time에는 이벤트까지 걸린 시간
(또는 이벤트가 없었다는 것을 마지막으로 확인한 시간),
event에는 이벤트 여부를 넣어줍니다.
Surv.obj<-Surv(time=dat3$RFS, event=dat3$Recur==1)이렇게 해서 Surv.obj라는 time to event 타입의 변수가
만들어졌습니다. 이것을 이용해서 먼저 Kaplan-Meier curve를 그려봅시다.
아래의 코드는 전체 환자의 생존을 하나의 곡선으로 표현하는 Kaplan-Meier
curve를 그리는 가장 간단한 방식입니다. survfit 함수에 이미
만든 time to event 타입 변수를 써주고 ~1을 써줍니다. 그
결과물을 ggsurvplot 함수 안에 넣어주면 그래프가
그려집니다.
fit<-survfit(Surv.obj~1, data=dat3)
ggsurvplot(fit)
그룹별로 Kaplan-Meier curve를 그리고 싶으면 survfit
함수안에 물결 표시 뒤에 1 대신, 그 그룹을 나타내는 factor type 변수를
써주면 됩니다.
fit<-survfit(Surv.obj~CA19.9.group, data=dat3)
ggsurvplot(fit)
이 그래프의 시간 변수 단위를 보기좋게 연도로 바꾸고, risk table을 넣고, 제목을 넣는 등의 변화도 줄 수 있습니다.
ggsurvplot(fit, title="Kaplan-Meier Curves", risk.table=T, xscale=365.25, break.x.by=365.25, xlab="Year",
legend=c(0.7, 0.5), legend.labs=c("Normal CA19-9", "Abnormal CA19-9"), legend.title="",
conf.int=T, pval=T, tables.y.text=F) 
\(~\)
Log-rank Test
위의 그래프를 만들 때 pval=T 옵션을 이용해서 log-rank
test 결과를 그래프에 나타나게 했습니다. 그래프를 그리지 않고 log-rank
test를 하려면 survdiff 함수 안에 ‘비교하려는 time to event
타입 변수 ~ 비교하려는 그룹’ 형식의 formula와 data 이름을 넣어주면
됩니다.
survdiff(Surv.obj~CA19.9.group, data=dat3)## Call:
## survdiff(formula = Surv.obj ~ CA19.9.group, data = dat3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## CA19.9.group=0 30 4 13.83 6.99 19
## CA19.9.group=1 20 18 8.17 11.83 19
##
## Chisq= 19 on 1 degrees of freedom, p= 1e-05survdiff(Surv.obj~sex, data=dat3)## Call:
## survdiff(formula = Surv.obj ~ sex, data = dat3)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=0 17 4 9.53 3.20 5.96
## sex=1 33 18 12.47 2.45 5.96
##
## Chisq= 6 on 1 degrees of freedom, p= 0.01\(~\)
Cox regression
Cox proportional hazards model을 fitting하는 방법을 살펴보겠습니다.
coxph() 함수안에 ‘time to event 타입 결과변수 ~ 설명변수’
형식의 formula를 넣어주고, data 이름을 넣어주면 됩니다. 간단하게 결과를
보려면 summary() 함수 안에 결과물을 넣으면 됩니다.
cox.obj<-coxph(Surv.obj~age, data=dat3)
summary(cox.obj)## Call:
## coxph(formula = Surv.obj ~ age, data = dat3)
##
## n= 50, number of events= 22
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.01944 1.01963 0.01981 0.982 0.326
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.02 0.9807 0.9808 1.06
##
## Concordance= 0.598 (se = 0.067 )
## Likelihood ratio test= 1 on 1 df, p=0.3
## Wald test = 0.96 on 1 df, p=0.3
## Score (logrank) test = 0.97 on 1 df, p=0.3로지스틱 회귀분석의 결과는 lmtest 패키지를 이용해서
결과를 정리하는 법을 보여드렸죠. Cox regression의 경우
broom 이라는 패키지를 쓰면 결과를 쉽게 뽑아내고 정리할 수
있습니다. (broom 패키지는 로지스틱 회귀분석 결과에도 적용 가능하긴
하지만, confidence interval을 계산하는 방식이 정확도가 떨어져서 OR의
CI가 1 값을 포함하는지 여부가 p-value와 모순 될 때가 가끔 있어서
추천하지 않습니다. ) broom 패키지를 로드하고
tidy() 함수 안에 Cox regression 결과를 넣어주면 됩니다.
tidy(cox.obj)## # A tibble: 1 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 age 0.0194 0.0198 0.982 0.326지수함수를 취해서 beta coefficient가 아닌 HR을 뽑아내고, confidence interval을 구하고, 결과를 보기좋게 정리하는 코드입니다.
data.frame(tidy(cox.obj, exponentiate = T, conf.int = T)[, c("term", "estimate", "conf.low", "conf.high", "p.value")])## term estimate conf.low conf.high p.value
## 1 age 1.019634 0.9808067 1.059998 0.3262994Multivariable Cox model도 같은 방식으로 돌리고 결과를 정리하면
됩니다. 로지스틱 회귀분석 때와 마찬가지로,
na.action = na.exclude 옵션을 써주는 것이 좋습니다. (사실,
아래 코드로 돌리는 모델은 이 데이터셋의 사이즈에 비해 너무 복잡한
모델이라서 overfitting될 우려가 있는 모델입니다. 바람직한 모델이
아니지만, 코드 시연을 위해서 돌려보는 모델이라는 것을 참고하기
바랍니다.)
cox.obj<-coxph(Surv.obj~age+TNM.b+log.CA19.9, data=dat3, na.action = na.exclude)## Warning in coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
## Loglik converged before variable 2,3,4 ; coefficient may be infinite.data.frame(tidy(cox.obj, exponentiate = T, conf.int = T)[,
c("term", "estimate", "conf.low", "conf.high", "p.value")])## term estimate conf.low conf.high p.value
## 1 age 1.039392e+00 0.9975642 1.082973 0.065242299
## 2 TNM.b3 8.729025e+07 0.0000000 Inf 0.997930801
## 3 TNM.b4 3.280817e+08 0.0000000 Inf 0.997780967
## 4 TNM.b5 4.789020e+08 0.0000000 Inf 0.997738165
## 5 log.CA19.9 1.473379e+00 1.1758672 1.846166 0.000757892그런데 warning이 뜨고, 결과물에 무한대를 뜻하는 Inf가
나왔습니다. 왜 그럴까요? 이런 경우 대부분, 너무 작은 카테고리가 있어서
그 카테고리에서 이벤트가 없어서 그렇습니다. 너무 작은 카테고리를
사용하면 안되는 이유죠.
xtabs(~TNM.b + Recur, data=dat3)## Recur
## TNM.b 0 1
## 1or2 7 0
## 3 7 5
## 4 12 13
## 5 2 4TNM.b 변수의 1or2 카테고리에, 이벤트를 겪은 환자가 전혀 없었음을 알
수 있습니다. 이 문제를 해결하기 위해서 1or2 카테고리와 3 카테고리를
합친, TNM.c라는 변수를 만들어봅시다. 첫번째 세션에서 배웠던
mapvalues() 함수를 쓰겠습니다. 이 함수는 plyr
패키지에 들어있는 함수이므로 패키지를 먼저 로드해준 후 사용합니다.
library(plyr)##
## Attaching package: 'plyr'## The following object is masked from 'package:ggpubr':
##
## mutatedat3$TNM.c<-mapvalues(dat3$TNM.b, from=c("1or2", "3", "4", "5"),
to=c("3.or.lower", "3.or.lower", "4", "5"))
xtabs(~TNM.c + Recur, data=dat3)## Recur
## TNM.c 0 1
## 3.or.lower 14 5
## 4 12 13
## 5 2 4이제 모든 카테고리에 이벤트가 적어도 한개 이상 있습니다. 아까 돌렸던 모델에서 TNM.b 대신에 TNM.c 변수를 이용하여 다시 모델을 fitting해 봅시다.
cox.obj<-coxph(Surv.obj~age+TNM.c+log.CA19.9, data=dat3, na.action = na.exclude)
data.frame(tidy(cox.obj, exponentiate = T, conf.int = T)[,
c("term", "estimate", "conf.low", "conf.high", "p.value")])## term estimate conf.low conf.high p.value
## 1 age 1.039306 0.9974887 1.082877 6.577421e-02
## 2 TNM.c4 5.340499 1.6084383 17.732064 6.215741e-03
## 3 TNM.c5 7.991658 1.7895391 35.688853 6.485198e-03
## 4 log.CA19.9 1.547024 1.2444648 1.923143 8.507445e-05로지스틱 회귀분석에서, 카테고리가 3개 이상인 변수에 대해 overall p-value를 구했던 것 기억나시나요? 같은 방법으로 likelihood ratio test를 이용해서 TNM.c의 overall p-value를 구해봅시다.
cox.obj0<-update(cox.obj, .~.-TNM.c)
lrtest(cox.obj, cox.obj0)## Likelihood ratio test
##
## Model 1: Surv.obj ~ age + TNM.c + log.CA19.9
## Model 2: Surv.obj ~ age + log.CA19.9
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 4 -61.511
## 2 2 -66.840 -2 10.659 0.004847 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Pr(>Chisq) 밑에 있는 값 0.004847이 TNM.c의 overall p-value 입니다.
로지스틱 회귀분석에서 c-index를 이용하여 로지스틱 회귀모형의
분류성능을 요약했던 것처럼, survival model에 대해서도 모델이 얼마나
생존을 잘 예측하는지 평가하는 지표가 여러가지 있습니다. 그 중에서도 가장
많이 쓰이는 것이 Harrell’s c-index 입니다. 이것을 얻어내려면 Cox model의
결과물을 summary() 함수에 넣어도 되지만, 딱 c-index 값만
뽑아내려면 concordance() 함수를 사용합니다.
summary(cox.obj)## Call:
## coxph(formula = Surv.obj ~ age + TNM.c + log.CA19.9, data = dat3,
## na.action = na.exclude)
##
## n= 50, number of events= 22
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.03855 1.03931 0.02095 1.840 0.06577 .
## TNM.c4 1.67532 5.34050 0.61228 2.736 0.00622 **
## TNM.c5 2.07840 7.99166 0.76350 2.722 0.00649 **
## log.CA19.9 0.43633 1.54702 0.11104 3.930 8.51e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.039 0.9622 0.9975 1.083
## TNM.c4 5.340 0.1872 1.6084 17.732
## TNM.c5 7.992 0.1251 1.7895 35.689
## log.CA19.9 1.547 0.6464 1.2445 1.923
##
## Concordance= 0.765 (se = 0.055 )
## Likelihood ratio test= 25.96 on 4 df, p=3e-05
## Wald test = 22.37 on 4 df, p=2e-04
## Score (logrank) test = 25.58 on 4 df, p=4e-05concordance(cox.obj)$concordance## [1] 0.7652174Cox model에서 backward elimination을 시행하는 방법도 똑같습니다. 그런데 아래와 같은 에러가 뜨는 경우는, 결측이 있을 경우입니다. 결측이 있더라도 항상 이런 문제가 생기는 것은 아니고, 그 결측이 있었던 변수가 중간에 탈락하게 되면 갑자기 분석에 들어가는 환자 명수가 달라지면서 이렇게 오류가 납니다. (로지스틱 회귀나 선형 회귀분석에서도 똑같습니다.)
cox.obj.full<-coxph(Surv.obj~age+TNM.c+log.CA19.9+sex+adc_decrease,
data=dat3, na.action = na.exclude)
step(cox.obj.full, direction="backward")## Start: AIC=120.71
## Surv.obj ~ age + TNM.c + log.CA19.9 + sex + adc_decrease## Error in drop1.default(fit, scope$drop, scale = scale, trace = trace, : number of rows in use has changed: remove missing values?이런 경우, 결측이 몇명인가, 왜 결측이 있었나를 따져보고 결측인 환자를 분석에서 제외해도 문제가 없다는 판단이 내려질 경우, 결측을 제외하고 backward elimination을 실행하는 코드입니다.
cox.obj.full<-coxph(Surv.obj~age+TNM.c+log.CA19.9+sex+adc_decrease,
subset=!is.na(adc_decrease),
data=dat3, na.action = na.exclude)
step(cox.obj.full, direction="backward")## Start: AIC=120.71
## Surv.obj ~ age + TNM.c + log.CA19.9 + sex + adc_decrease
##
## Df AIC
## - sex 1 119.59
## <none> 120.71
## - adc_decrease 1 122.75
## - age 1 123.19
## - TNM.c 2 123.97
## - log.CA19.9 1 136.46
##
## Step: AIC=119.59
## Surv.obj ~ age + TNM.c + log.CA19.9 + adc_decrease
##
## Df AIC
## <none> 119.59
## - age 1 122.44
## - adc_decrease 1 123.42
## - TNM.c 2 126.95
## - log.CA19.9 1 137.46## Call:
## coxph(formula = Surv.obj ~ age + TNM.c + log.CA19.9 + adc_decrease,
## data = dat3, subset = !is.na(adc_decrease), na.action = na.exclude)
##
## coef exp(coef) se(coef) z p
## age 0.04505 1.04608 0.02072 2.174 0.02967
## TNM.c4 2.03454 7.64873 0.69276 2.937 0.00332
## TNM.c5 2.08350 8.03251 0.86496 2.409 0.01601
## log.CA19.9 0.62327 1.86502 0.14429 4.319 1.56e-05
## adc_decrease 0.03904 1.03981 0.01634 2.389 0.01689
##
## Likelihood ratio test=30.57 on 5 df, p=1.139e-05
## n= 48, number of events= 21backward elimination을 통해 최종 선택된 모델을
cox.obj.b에 저장하고 그 결과물을 뽑아냅시다.
cox.obj.b<-step(cox.obj.full, direction="backward")## Start: AIC=120.71
## Surv.obj ~ age + TNM.c + log.CA19.9 + sex + adc_decrease
##
## Df AIC
## - sex 1 119.59
## <none> 120.71
## - adc_decrease 1 122.75
## - age 1 123.19
## - TNM.c 2 123.97
## - log.CA19.9 1 136.46
##
## Step: AIC=119.59
## Surv.obj ~ age + TNM.c + log.CA19.9 + adc_decrease
##
## Df AIC
## <none> 119.59
## - age 1 122.44
## - adc_decrease 1 123.42
## - TNM.c 2 126.95
## - log.CA19.9 1 137.46data.frame(tidy(cox.obj.b, exponentiate = T, conf.int = T)[,
c("term", "estimate", "conf.low", "conf.high", "p.value")])## term estimate conf.low conf.high p.value
## 1 age 1.046084 1.004453 1.089439 2.967373e-02
## 2 TNM.c4 7.648725 1.967486 29.734895 3.315486e-03
## 3 TNM.c5 8.032505 1.474332 43.762953 1.600608e-02
## 4 log.CA19.9 1.865020 1.405598 2.474604 1.563906e-05
## 5 adc_decrease 1.039810 1.007037 1.073651 1.689015e-02\(~\)
time-dependent ROC analysis
생존을 예측하는 모델을 만들고 나서 그 예측모델의 성능을 평가하는 방법 중 대표적인 것이 위에서 다룬 Harrell’s c-index가 있고, 또 time-dependent ROC analysis라는 것도 많이 쓰입니다. 이분형 결과변수일 때와 마찬가지로 예측모델에 내놓는 위험도의 가능한 모든 cutoff에 대해서 sensitivity와 specificity를 계산해서 요약한다는 점에서는 일반적인 ROC analysis와 같습니다. 그런데, 일단 가장 큰 다른 점은, 생존분석에서는 시간이 지남에 따라 환자별 이벤트여부가 달라진다는 점에 있습니다. 그래서 일단 특정 시점을 지정해줘야만, 그 시점에서의 생존/사망 여부에 대한 ROC analysis를 할 수가 있습니다. 두번째로 다른 점은, censoring 때문에 sensitivity, specificity 계산이 어렵다는 점입니다. 특정 시점에서 어떤 환자들은 생존했고 어떤 환자들은 사망했고 어떤 환자들은 censor되어서 생존/사망 여부를 모른다는 점입니다. 그래서 이렇게 censor된 환자 데이터를 어떻게 sensitivity, specificity 계산에 반영하느냐에 따라 여러가지 time-dependent ROC analysis 방법이 개발되어있습니다. 자세한 것은 이 리뷰페이퍼에 잘 정리되어있습니다. https://bmcmedresmethodol.biomedcentral.com/articles/10.1186/s12874-017-0332-6
위의 페이퍼에 나온 방법 중 CD5 방법으로 time-dependent ROC analysis를
하는 코드를 소개드리려고 합니다. 먼저 로지스틱 회귀모델을 ROC 분석할
때처럼, Cox model이 예측한 각 환자별 위험도를 뽑아내야 합니다. 로지스틱
회귀분석 때와 같이 predict() 함수를 사용하면 되는데,
로지스틱 회귀분석과는 달리 이벤트 확률(=사망확률)을 계산해주지는
않습니다. 이것은 Cox proportional hazard model이 baseline hazard 또는
baseline survival probability를 추정하지 않은채로 “각 환자의 생존확률이
[baseline 생존확률] (=설명변수 값이 모두 0인 환자의 생존확률)의
몇제곱인가”만 추정한다는 특성 때문입니다. 자세한 이야기는 이 강의의
범위를 벗어나니 생존자료분석 강의를 들으시기를 추천합니다. 아무튼 그런
이유로, linear predictor (beta coefficients와 설명변수들의 선형결합) 를
뽑아내서 time-dependent ROC analysis를 합니다. 아래의 코드는 모델이
추정한 각 환자의 linear predictor를 cox.lp에 저장합니다.
cox.obj<-coxph(Surv.obj~age+TNM.c+log.CA19.9, data=dat3, na.action = na.exclude)
cox.lp<-predict(cox.obj, type="lp")위 논문에 나온 방법들 중 CD5 방법으로 time-dependent ROC analysis를
하는 방법을 알아보겠습니다. 이 방법은 Inverse Probability of Censoring
방법으로 데이터에 웨이트를 주어서 sensitivity와 specificity를 계산하는
방법입니다. timeROC 패키지의 timeROC() 함수를
쓰면 되는데, T는 time to event, delta는 이벤트
여부입니다. 단, delta에는 factor type을 넣으면 오작동
하네요. 꼭 numeric type으로 바꿔서 넣어주셔야겠습니다.
marker는 역시 Cox 모델에서 나온 linear predictor를 넣으시면
되고, times는 ROC 분석을 하는 관심 시점을 넣어줍니다.
cause는 delta에 넣은 값 중 어떤 것이 이벤트를
의미하는지 넣어줍니다. iid=TRUE는 AUC를 계산할 때
iid-representation을 계산한다는 뜻인데, 이 옵션을 쓰지 않으면 디폴트로
iid=FALSE가 되고 그러면 confidence interval을 계산할 수
없게 됩니다.
t.ROC<-timeROC(T=dat3$RFS, delta=as.double(as.character(dat3$Recur)),
marker=cox.lp, times=365.25, cause=1, iid=TRUE)
t.ROC$AUC## t=0 t=365.25
## NA 0.7657352confint(t.ROC)## $CI_AUC
## 2.5% 97.5%
## t=365.25 60.45 92.69
##
## $CB_AUC
## 2.5% 97.5%
## t=365.25 60.42 92.73
##
## $C.alpha
## 95%
## 1.964532위와 같이 confint() 함수를 이용해서 AUC의 신뢰구간
계산이 가능합니다. (CI_AUC 아래 나온 값입니다.) 그리고 plot
함수에 넣으면 ROC 커브를 쉽게 그릴 수 있는데, 단, 꼭 시점을 지정해주셔야
합니다.
plot(t.ROC, time=365.25)
이렇게 생존자료를 카플란 마이어 커브로 요약하고, Cox regression를 수행하고 성능을 평가하는 방법을 알아보았습니다. Cox regression의 중요한 가정이 있는데 바로 proportional hazards assumption 입니다. 이 가정이 무슨 뜻인지, 어떻게 체크하는지, 이 가정이 어긋날 경우 어떻게 해야하는지는 이 강의에서 다루지 않았지만 꼭 공부해보시기 바랍니다.