Survival Data Analysis
Seo Young Park, PhD
\(~\)
Survival Data Analysis란 time to event 타입의 변수를 결과변수로 하는 분석을 말합니다. time to event 타입의 변수란, (1) 이벤트 여부와 (2) 이벤트까지 걸린 시간 (또는 이벤트가 없었다는 것을 마지막으로 확인한 시간) 두가지로 구성된 변수를 말합니다. 임상연구에서는 보통 사망이나 재발과 같은 이벤트의 위험에 대해 분석할 때 쓰게 되죠. 그러면 R로 Survival Data Analysis를 하는 방법을 살펴보겠습니다.
\(~\)
Kaplan-Meier Curve 그리기
일단 R에서 생존자료를 다루려면 survival 이라는 패키지를 이용해야합니다. 이 패키지의 Surv 함수를 이용하여 time to event 타입의 변수를 인식시킬 수 있습니다. 이 강의에서 사용하고 있는 예제 데이터를 불러와서, time to event 타입의 변수를 만들어봅시다. Surv 함수 안에 time과 event 두가지 값을 넣어줘야 하는데, time에는 이벤트까지 걸린 시간 (또는 이벤트가 없었다는 것을 마지막으로 확인한 시간), event에는 이벤트 여부를 넣어줍니다.
load("clean_data.RData")
library(survival)
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 함수 안에 넣어주면 그래프가 그려집니다.
library(survminer)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) ## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2.
\(~\)
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로지스틱 회귀분석을 할 때 OR는 자동으로 계산되지 않았던 것처럼, Cox regression을 할 때도 HR (Hazard Ratio)가 자동으로 계산되지 않습니다. 로지스틱 회귀분석의 결과는 lmtest 패키지를 이용해서 결과를 정리하는 법을 보여드렸죠. Cox regression의 경우 broom 이라는 패키지를 쓰면 결과를 쉽게 뽑아내고 정리할 수 있습니다. (broom 패키지는 로지스틱 회귀분석 결과에도 적용 가능하긴 하지만, confidence interval을 계산하는 방식이 정확도가 떨어져서 OR의 CI가 1 값을 포함하는지 여부가 p-value와 모순 될 때가 가끔 있어서 추천하지 않습니다. ) broom 패키지를 로드하고 tidy() 함수 안에 Cox regression 결과를 넣어주면 됩니다.
library(broom)
cox.obj<-coxph(Surv.obj~age, data=dat3)
tidy(cox.obj)## # A tibble: 1 x 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 옵션을 써주는 것이 좋습니다.
cox.obj<-coxph(Surv.obj~age+CCI.b+log.CA19.9+adc_decrease, 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.0100567 0.9707602 1.050944 0.6211417649
## 2 CCI.b2 1.8591860 0.6408277 5.393919 0.2538182232
## 3 CCI.b3 1.0821170 0.2872995 4.075807 0.9071472163
## 4 CCI.b>=4 0.9771283 0.1133027 8.426802 0.9832076611
## 5 log.CA19.9 1.6386821 1.2638068 2.124754 0.0001941537
## 6 adc_decrease 1.0275045 0.9976835 1.058217 0.0709761761카테고리가 4개인 범주형 변수 CCI.b의 overall p-value를 알아내는 법도 로지스틱 회귀분석과 똑같습니다. update() 함수가 Cox 모델에도 똑같이 적용됩니다.
cox.obj0<-update(cox.obj, .~.-CCI.b)
anova(cox.obj, cox.obj0)## Analysis of Deviance Table
## Cox model: response is Surv.obj
## Model 1: ~ age + CCI.b + log.CA19.9 + adc_decrease
## Model 2: ~ age + log.CA19.9 + adc_decrease
## loglik Chisq Df P(>|Chi|)
## 1 -59.719
## 2 -60.474 1.5097 3 0.68이 모델의 Harrell’s c-index를 알아내려면 Cox model의 결과물을 summary 함수에 넣어도 되지만, 딱 c-index 값만 뽑아내려면 concordance() 함수를 사용합니다.
summary(cox.obj)## Call:
## coxph(formula = Surv.obj ~ age + CCI.b + log.CA19.9 + adc_decrease,
## data = dat3, na.action = na.exclude)
##
## n= 48, number of events= 21
## (2 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.01001 1.01006 0.02025 0.494 0.621142
## CCI.b2 0.62014 1.85919 0.54345 1.141 0.253818
## CCI.b3 0.07892 1.08212 0.67662 0.117 0.907147
## CCI.b>=4 -0.02314 0.97713 1.09928 -0.021 0.983208
## log.CA19.9 0.49389 1.63868 0.13254 3.727 0.000194 ***
## adc_decrease 0.02713 1.02750 0.01503 1.806 0.070976 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.0101 0.9900 0.9708 1.051
## CCI.b2 1.8592 0.5379 0.6408 5.394
## CCI.b3 1.0821 0.9241 0.2873 4.076
## CCI.b>=4 0.9771 1.0234 0.1133 8.427
## log.CA19.9 1.6387 0.6102 1.2638 2.125
## adc_decrease 1.0275 0.9732 0.9977 1.058
##
## Concordance= 0.799 (se = 0.056 )
## Likelihood ratio test= 20.72 on 6 df, p=0.002
## Wald test = 20.49 on 6 df, p=0.002
## Score (logrank) test = 25.78 on 6 df, p=2e-04concordance(cox.obj)$concordance## [1] 0.7990431종종 아래와 같이 warning이 뜨고, 결과물에 무한대를 뜻하는 Inf가 뜰 때가 있습니다. 왜 그럴까요?
cox.obj<-coxph(Surv.obj~age+TNM.b+log.CA19.9+adc_decrease, data=dat3, na.action = na.exclude)## Warning in fitter(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.048997e+00 1.008673 1.090934 0.0167687524
## 2 TNM.b3 2.960113e+08 0.000000 Inf 0.9975277502
## 3 TNM.b4 1.074234e+09 0.000000 Inf 0.9973643826
## 4 TNM.b5 1.092070e+09 0.000000 Inf 0.9973622955
## 5 log.CA19.9 1.835940e+00 1.373204 2.454607 0.0000412487
## 6 adc_decrease 1.058245e+00 1.019714 1.098231 0.0027750530이런 경우 대부분, 너무 작은 카테고리가 있어서 그 카테고리에서 이벤트가 없어서 그렇습니다. 너무 작은 카테고리를 사용하면 안되는 이유죠.
table(dat3$TNM.b, dat3$Recur)##
## 0 1
## 1or2 7 0
## 3 7 5
## 4 12 13
## 5 2 4Cox model에서 backward elimination을 시행하는 방법도 똑같습니다. 그런데 아래와 같은 에러가 뜨는 경우는, 결측이 있을 경우입니다. 결측이 있더라도 항상 이런 문제가 생기는 것은 아니고, 그 결측이 있었던 변수가 중간에 탈락하게 되면 갑자기 분석에 들어가는 환자 명수가 달라지면서 이렇게 오류가 납니다. (로지스틱 회귀나 선형 회귀분석에서도 똑같습니다.)
cox.obj<-coxph(Surv.obj~age+CCI.b+log.CA19.9+adc_decrease+adc_1_mean, data=dat3, na.action = na.exclude)
step(cox.obj, direction="backward")## Start: AIC=129.78
## Surv.obj ~ age + CCI.b + log.CA19.9 + adc_decrease + adc_1_mean## Error in drop1.default(fit, scope$drop, scale = scale, trace = trace, : number of rows in use has changed: remove missing values?이런 경우, 결측이 몇명인가, 왜 결측이 있었나를 따져보고 결측인 환자를 분석에서 제외해도 문제가 없다는 판단이 내려질 경우, 결측을 제외하고 backward elimination을 실행하는 코드입니다. backward elimination을 통해 최종 선택된 모델을 cox.obj.b에 저장하고 그 결과물을 뽑아냈습니다.
cox.obj<-coxph(Surv.obj~age+CCI.b+log.CA19.9+adc_decrease+adc_1_mean, data=dat3,
subset=!is.na(adc_decrease) & !is.na(adc_1_mean))
cox.obj.b<-step(cox.obj, direction="backward", trace=0)
summary(cox.obj.b)## Call:
## coxph(formula = Surv.obj ~ log.CA19.9 + adc_decrease, data = dat3,
## subset = !is.na(adc_decrease) & !is.na(adc_1_mean))
##
## n= 46, number of events= 21
##
## coef exp(coef) se(coef) z Pr(>|z|)
## log.CA19.9 0.49146 1.63470 0.12128 4.052 5.07e-05 ***
## adc_decrease 0.02273 1.02299 0.01362 1.669 0.0951 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## log.CA19.9 1.635 0.6117 1.289 2.073
## adc_decrease 1.023 0.9775 0.996 1.051
##
## Concordance= 0.764 (se = 0.06 )
## Likelihood ratio test= 17.13 on 2 df, p=2e-04
## Wald test = 16.98 on 2 df, p=2e-04
## Score (logrank) test = 18.77 on 2 df, p=8e-05data.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 log.CA19.9 1.634703 1.2888690 2.073333 5.068886e-05
## 2 adc_decrease 1.022990 0.9960479 1.050661 9.508519e-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
위의 페이퍼에 나온 방법 중 CD1 방법과 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+CCI.b+log.CA19.9+adc_decrease, data=dat3, na.action = na.exclude)
cox.lp<-predict(cox.obj, type="lp")\(~\)
먼저, 위의 링크에 있는 리뷰페이퍼에 나온 방법 중 CD1 방법은 survivalROC 패키지의 survivalROC() 함수를 이용합니다. Stime은 time to event이고, status는 event 여부입니다. marker는 생존예측모델의 output을 말하는데, cox.lp를 넣어주면 됩니다. predict.time은 어떤 시점에 대해서 ROC 분석을 할지 시점을 지정해주는 옵션입니다. method="KM"은 censor된 환자를 계산에 반영할 때 Kaplan-Meier 방법으로 추정한 생존함수를 이용하라는 뜻인데, 이렇게 해야 CD1 방법이 됩니다. 이렇게 만든 ROC 커브의 area under the curve를 구하는 방법입니다.
library(survivalROC)
t.ROC<-survivalROC(Stime=dat3$RFS, status=dat3$Recur, marker=cox.lp, predict.time=365.25, method="KM")##
## 2 records with missing values dropped.t.ROC$AUC## [1] 0.8533689위의 결과를 이용해서 아래와 같이 ROC 커브를 그릴 수도 있습니다. 아주아주 자세히 보시면 이것이 monotone함수가 아니라는 것을 알 수 있습니다. 원래 ROC 커브는 이론적으로는 monotone이어야 하는데 말이죠. 이것이 CD1 방법의 큰 단점입니다.
plot(t.ROC$FP, t.ROC$TP, type="l", xlab="1 - Specificity", ylab="Sensitivity")
\(~\)
다음은 위 논문에 나온 방법들 중 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을 계산할 수 없게 됩니다.
library(timeROC)
t.ROC2<-timeROC(T=dat3$RFS, delta=as.double(as.character(dat3$Recur)), marker=cox.lp, times=365.25, cause=1, iid=TRUE)
t.ROC2$AUC## t=0 t=365.25
## NA 0.8662456confint(t.ROC2)## $CI_AUC
## 2.5% 97.5%
## t=365.25 71.43 101.82
##
## $CB_AUC
## 2.5% 97.5%
## t=365.25 71.94 101.31
##
## $C.alpha
## 95%
## 1.894311위와 같이 confint() 함수를 이용해서 AUC의 신뢰구간 계산이 가능합니다. (CI_AUC 아래 나온 값입니다.) 그리고 plot 함수에 넣으면 ROC 커브를 쉽게 그릴 수 있는데, 단, 꼭 시점을 지정해주셔야 합니다.
plot(t.ROC2, time=365.25)
이 CD5 방법은 항상 monotone인 ROC 커브를 만들어냅니다.
\(~\)
Survival IDI, NRI
Cox 모델을 이용한 예측모델에서도 새롭게 더해진 설명변수의 added value를 NRI, IDI를 이용해 수치화할 수 있습니다. 이것은 survIDINRI 라는 패키지를 통해서 계산할 수 있는데요, 단, 모든 변수에 결측이 없는 경우만 가능합니다. 따라서 여기에서는 결측이 없는 age, CCI.b 두개만 설명변수로 가지고 있는 모델에,log.CA19.9가 설명변수로 추가될 경우의 added value를 계산해보겠습니다.
일단 이 패키지의 IDI.INF() 함수를 이용할 건데요, 재미있는 것은 이 함수에 두개의 Cox model의 결과물을 넣는 것이 아니라, 두개의 Cox model에 들어가는 설명변수들의 데이터를 행렬 형태로 넣어야 한다는 점입니다. (함수 내부에서 Cox model을 자동으로 fitting합니다.) data frame 형태는 받아들이지 않고 모든 값이 numeric type인 행렬만 받아들이기 때문에, model.matrix() 라는 함수를 써서 행렬을 만드는 것이 가장 간단합니다. 단, 이 model.matrix() 함수는 원래 회귀분석의 design matrix를 만드는 목적으로 쓰이는 함수라서, 제 1열의 모든 값이 1인, 우리에게 필요없는 부분을 가지고 있습니다. 그래서 제 1열을 제거해준 값을 사용해야합니다. (왜 그런지 설명하자니 너무 복잡한 디테일인 것 같네요… 강의에서 간단히 말씀드리겠습니다.)
아래에서 yy는 이벤트까지 걸린 시간, 이벤트 여부의 행렬입니다. xx0는 age, CCI.b 두개의 설명변수의 행렬, xx1는 age, CCI.b,log.CA19.9 세 개의 행렬입니다. 이것들을 각각 indata, covs0, covs1에 넣고, t0로 관심있는 특정 시점을 지정해줍니다. seed1은 신뢰구간을 계산할 때 쓰는 resampling의 random seed입니다. 이 값은 아무 값으로 해도 상관없는데, 값을 꼭 지정해주어야만 다음번에 같은 결과를 reproduce할 수 있습니다. 결과물을 저장한 객체를 IDI.INF.OUT() 함수 안에 넣어야 결과를 볼 수 있습니다.
library(survIDINRI)## Loading required package: survC1yy<-cbind(dat3$RFS, as.double(as.character(dat3$Recur)))
xx0<-model.matrix(~age+CCI.b, data=dat3)[, -1]
xx1<-model.matrix(~age+CCI.b+log.CA19.9, data=dat3)[, -1]
obj<-IDI.INF(indata=yy, covs0=xx0, covs1=xx1, t0=365.25, seed1=2021)
IDI.INF.OUT(obj)## Est. Lower Upper p-value
## M1 0.153 -0.025 0.387 0.113
## M2 0.257 -0.253 0.620 0.306
## M3 0.077 -0.090 0.508 0.452M1이 IDI, M2가 NRI입니다.
\(~\)
Optimisim-corrected c-index
저번 로지스틱 회귀모형을 돌린 후 optimism-corrected AUC를 구하는 방법을 제가 만든 긴 코드로 보여드렸습니다. Cox 모델을 돌린 후 optimisim-corrected c-index (Harrell’s c-index)를 구하는 것은 다행히 R 패키지로 구현되어있어서 간단하게 구할 수 있습니다. c-index 자체를 계산해주는 것은 아니라 약간의 후처리가 필요하지만 직접 프로그램하는 것보다는 훨씬 간단합니다.
일단 rms라는 패키지를 설치하고 로드한 후, 관심있는 Cox regression 모델을 coxph() 함수가 아닌, rms 패키지의 함수인 cph()를 이용해서 돌려주어야 합니다. cph()는 coxph() 와 사용법이 거의 같습니다. 단, optimism-corrected c-index 계산을 위해서는 x=T, y=T 라는 구문을 꼭 넣어주어야 합니다.
cox.obj<-coxph(Surv.obj~age+log.CA19.9+adc_decrease, data=dat3, na.action = na.exclude)
summary(cox.obj)## Call:
## coxph(formula = Surv.obj ~ age + log.CA19.9 + adc_decrease, data = dat3,
## na.action = na.exclude)
##
## n= 48, number of events= 21
## (2 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.01469 1.01480 0.01977 0.743 0.4573
## log.CA19.9 0.52030 1.68253 0.12418 4.190 2.79e-05 ***
## adc_decrease 0.02500 1.02532 0.01453 1.721 0.0852 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.015 0.9854 0.9762 1.055
## log.CA19.9 1.683 0.5943 1.3190 2.146
## adc_decrease 1.025 0.9753 0.9965 1.055
##
## Concordance= 0.791 (se = 0.054 )
## Likelihood ratio test= 19.21 on 3 df, p=2e-04
## Wald test = 19.35 on 3 df, p=2e-04
## Score (logrank) test = 21.78 on 3 df, p=7e-05concordance(cox.obj)$concordance## [1] 0.7910686library(rms)
cph.obj<-cph(Surv.obj~age+log.CA19.9+adc_decrease, data=dat3, x=T, y=T)
cph.obj## Frequencies of Missing Values Due to Each Variable
## Surv.obj age log.CA19.9 adc_decrease
## 0 0 0 2
##
## Cox Proportional Hazards Model
##
## cph(formula = Surv.obj ~ age + log.CA19.9 + adc_decrease, data = dat3,
## x = T, y = T)
##
##
## Model Tests Discrimination
## Indexes
## Obs 48 LR chi2 19.21 R2 0.349
## Events 21 d.f. 3 Dxy 0.582
## Center 3.154 Pr(> chi2) 0.0002 g 1.180
## Score chi2 21.78 gr 3.254
## Pr(> chi2) 0.0001
##
## Coef S.E. Wald Z Pr(>|Z|)
## age 0.0147 0.0198 0.74 0.4573
## log.CA19.9 0.5203 0.1242 4.19 <0.0001
## adc_decrease 0.0250 0.0145 1.72 0.0852
## 자 이제 결과물을 validate() 함수에 넣어주면 되는데요, 주의할 점은 validate() 함수를 돌리기 전에 꼭 random seed 번호를 지정해주어야 한다는 것입니다. 왜냐하면 optimism-corrected c-index를 계산하는 과정에 bootstrap resampling은 random하게 데이터를 resampling하는 것이기 때문에, random seed를 지정하지 않으면 이 계산을 실행할 때마다 다른 결과를 얻게되기 때문입니다. reproducible research를 위해서 항상 random seed를 지정하고 저장하는 것이 좋습니다.
set.seed(2021)
vobj<-validate(cph.obj, B=500)
vobj## index.orig training test optimism index.corrected n
## Dxy 0.5821 0.6083 0.5562 0.0521 0.5300 500
## R2 0.3487 0.3970 0.3151 0.0819 0.2668 500
## Slope 1.0000 1.0000 0.8135 0.1865 0.8135 500
## D 0.1299 0.1612 0.1144 0.0468 0.0832 500
## U -0.0143 -0.0146 0.0606 -0.0752 0.0609 500
## Q 0.1442 0.1758 0.0538 0.1220 0.0222 500
## g 1.1798 1.4428 1.0889 0.3539 0.8259 500여기서 index.orig은 apparent 버전, 즉 원 데이터에 그대로 적용했을 때의 값이고, index.corrected가 바로 optimism-corrected 버전입니다.
자, 이 output 중 Dxy 라는 값을 이용해서 c-index를 계산할 수 있습니다. Dxy 값에 대해서 도움말을 보시면 이렇게 나와있습니다
“The values corresponding to the row Dxy are equal to 2 * (C - 0.5) where C is the C-index or concordance probability”
즉 Dxy = 2*(C - 0.5)라는 거죠. 이것을 C에 대해서 풀면, C=Dxy/2 + 0.5 입니다. 따라서 위의 결과에 나온 Dxy 값을 2로 나누고 0.5를 더하면 그게 Harrell’s c-index라는 뜻입니다. 일단, Dxy의 index.orig 값을 가지고, 원 데이터에서 계산한 c-index와 일치하는지 볼까요?
concordance(cox.obj)$concordance## [1] 0.7910686vobj["Dxy", "index.orig"]/2+0.5## [1] 0.7910686너무나 아름답게 완벽히 일치하죠!
그러면, optimism-corrected c-index는 아래와 같이 구할 수 있습니다.
vobj["Dxy", "index.corrected"]/2+0.5## [1] 0.7650143이렇게 rms 패키지의 validate() 함수를 이용하면 직접 bootstrap을 구현하는 프로그램을 짜지 않고도 쉽게 optimism-corrected c-index를 구할 수 있습니다.