Evaluating performance for the Cox model
Load the required libraries:
# load packages
library(knitr)
library(tidyverse)
library(lubridate)
library(survival)
library(survminer)
library(gtsummary)
library(ISwR)
library(KMsurv)
library(timeROC)
library(rms)Load the dataset melanom:
data("melanom", package="ISwR")
str(melanom)## 'data.frame': 205 obs. of 6 variables:
## $ no : int 789 13 97 16 21 469 685 7 932 944 ...
## $ status: int 3 3 2 3 1 1 1 1 3 1 ...
## $ days : int 10 30 35 99 185 204 210 232 232 279 ...
## $ ulc : int 1 2 2 2 1 1 1 1 1 1 ...
## $ thick : int 676 65 134 290 1208 484 516 1288 322 741 ...
## $ sex : int 2 2 2 1 2 2 2 2 1 1 ...head(melanom)## no status days ulc thick sex
## 1 789 3 10 1 676 2
## 2 13 3 30 2 65 2
## 3 97 2 35 2 134 2
## 4 16 3 99 2 290 1
## 5 21 1 185 1 1208 2
## 6 469 1 204 1 484 2Discrimination for the survival data : time-dependent ROC curve
First of all we evaluate the discrimination accuracy of the tumor thickness as a prognostic biomarker for death from malignant melanoma:
ROC.thick<-timeROC(T=melanom$days,delta=melanom$status,marker=melanom$thick,cause=1, times=c(1800,2000,2200),iid=TRUE)
ROC.thick## Time-dependent-Roc curve estimated using IPCW (n=205, with competing risks).
## Cases Survivors Other events Censored AUC_1 (%) se_1 AUC_2 (%) se_2
## t=1800 45 124 36 0 75.39 4.10 76.51 3.85
## t=2000 46 103 56 0 75.60 4.14 76.11 3.82
## t=2200 50 83 72 0 72.55 4.40 74.33 3.81
##
## Method used for estimating IPCW:marginal
##
## Total computation time : 0.12 secs.Note that since we are here in a competing risk setting (two causes of death, but we are interested in cause=1), we obtain two AUCs estimated each one for a risk. We are now interested in evaluating AUC_1. We can also evaluate the confidence intervals for the estimated AUCs:
confint(ROC.thick)## $CI_AUC_1
## 2.5% 97.5%
## t=1800 67.36 83.41
## t=2000 67.49 83.72
## t=2200 63.92 81.19
##
## $CB_AUC_1
## 2.5% 97.5%
## t=1800 66.47 84.30
## t=2000 66.59 84.61
## t=2200 62.97 82.14
##
## $C.alpha.1
## 95%
## 2.175971
##
## $CI_AUC_2
## 2.5% 97.5%
## t=1800 68.96 84.05
## t=2000 68.63 83.60
## t=2200 66.86 81.79
##
## $CB_AUC_2
## 2.5% 97.5%
## t=1800 68.47 84.55
## t=2000 68.14 84.08
## t=2200 66.38 82.28
##
## $C.alpha.2
## 95%
## 2.087521Then, we estimate a Cox model by adding sex to the tumor thickness:
mod.cov <- coxph(Surv(days,status==1)~sex+log(thick), data=melanom)
summary(mod.cov)## Call:
## coxph(formula = Surv(days, status == 1) ~ sex + log(thick), data = melanom)
##
## n= 205, number of events= 57
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sex 0.4580 1.5809 0.2687 1.705 0.0883 .
## log(thick) 0.7809 2.1834 0.1573 4.963 6.94e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex 1.581 0.6326 0.9337 2.677
## log(thick) 2.183 0.4580 1.6040 2.972
##
## Concordance= 0.749 (se = 0.033 )
## Likelihood ratio test= 33.45 on 2 df, p=5e-08
## Wald test = 31 on 2 df, p=2e-07
## Score (logrank) test = 32.52 on 2 df, p=9e-08We now extract the linear predictor from the Cox model (X*beta): note that the linear predictor of this cause-specific Cox model is evaluated at baseline and it is independent of time.
lin.pred <- predict(mod.cov)And we now compute the time-dependent discrimination from this model:
ROC.cox <- timeROC(T=melanom$days,delta=melanom$status,marker=lin.pred,cause=1, times=c(1800,2000,2200),iid=TRUE)
ROC.cox## Time-dependent-Roc curve estimated using IPCW (n=205, with competing risks).
## Cases Survivors Other events Censored AUC_1 (%) se_1 AUC_2 (%) se_2
## t=1800 45 124 36 0 75.76 4.09 76.94 3.85
## t=2000 46 103 56 0 76.17 4.15 76.39 3.83
## t=2200 50 83 72 0 74.33 4.34 75.13 3.75
##
## Method used for estimating IPCW:marginal
##
## Total computation time : 0.1 secs.along with the estimated confidence intervals:
confint(ROC.cox)## $CI_AUC_1
## 2.5% 97.5%
## t=1800 67.75 83.77
## t=2000 68.04 84.30
## t=2200 65.82 82.83
##
## $CB_AUC_1
## 2.5% 97.5%
## t=1800 66.92 84.60
## t=2000 67.19 85.15
## t=2200 64.94 83.71
##
## $C.alpha.1
## 95%
## 2.16351
##
## $CI_AUC_2
## 2.5% 97.5%
## t=1800 69.40 84.49
## t=2000 68.87 83.90
## t=2200 67.77 82.48
##
## $CB_AUC_2
## 2.5% 97.5%
## t=1800 68.84 85.05
## t=2000 68.32 84.45
## t=2200 67.23 83.02
##
## $C.alpha.2
## 95%
## 2.10383Now we can compare the AUCs: is there a significant improvement in discrimination by adding sex?
pp1 <- compare(ROC.thick,ROC.cox)
round(pp1$p_values_AUC_1,digits=4)## t=1800 t=2000 t=2200
## 0.7713 0.6594 0.1936No significant increment in discrimination is observed from this test.
Now we add also ulcer in the Cox model but this time as a covariate (not as a strata variable):
mod.cov2 <- coxph(Surv(days,status==1)~sex+log(thick)+ulc, data=melanom)
summary(mod.cov2)## Call:
## coxph(formula = Surv(days, status == 1) ~ sex + log(thick) +
## ulc, data = melanom)
##
## n= 205, number of events= 57
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sex 0.3813 1.4641 0.2706 1.409 0.15879
## log(thick) 0.5756 1.7782 0.1794 3.209 0.00133 **
## ulc -0.9389 0.3911 0.3243 -2.895 0.00379 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex 1.4641 0.6830 0.8615 2.4883
## log(thick) 1.7782 0.5624 1.2511 2.5273
## ulc 0.3911 2.5571 0.2071 0.7384
##
## Concordance= 0.768 (se = 0.032 )
## Likelihood ratio test= 42.66 on 3 df, p=3e-09
## Wald test = 36.99 on 3 df, p=5e-08
## Score (logrank) test = 42.81 on 3 df, p=3e-09We now extract also from this model the linear predictor (X*beta):
lin.pred2 <- predict(mod.cov2)And now we compare the performance of this model with the tumor thickness:
ROC.cox2 <- timeROC(T=melanom$days,delta=melanom$status,marker=lin.pred2,cause=1, times=c(1800,2000,2200),iid=TRUE)
ROC.cox2## Time-dependent-Roc curve estimated using IPCW (n=205, with competing risks).
## Cases Survivors Other events Censored AUC_1 (%) se_1 AUC_2 (%) se_2
## t=1800 45 124 36 0 78.74 3.95 79.36 3.71
## t=2000 46 103 56 0 78.59 4.06 78.58 3.73
## t=2200 50 83 72 0 75.88 4.27 76.66 3.69
##
## Method used for estimating IPCW:marginal
##
## Total computation time : 0.12 secs.confint(ROC.cox2)## $CI_AUC_1
## 2.5% 97.5%
## t=1800 71.00 86.47
## t=2000 70.64 86.54
## t=2200 67.51 84.25
##
## $CB_AUC_1
## 2.5% 97.5%
## t=1800 70.18 87.29
## t=2000 69.80 87.38
## t=2200 66.62 85.13
##
## $C.alpha.1
## 95%
## 2.167015
##
## $CI_AUC_2
## 2.5% 97.5%
## t=1800 72.10 86.62
## t=2000 71.27 85.88
## t=2200 69.43 83.89
##
## $CB_AUC_2
## 2.5% 97.5%
## t=1800 71.63 87.09
## t=2000 70.80 86.35
## t=2200 68.97 84.35
##
## $C.alpha.2
## 95%
## 2.08551pp1 <- compare(ROC.thick,ROC.cox2)
round(pp1$p_values_AUC_1,digits=4)## t=1800 t=2000 t=2200
## 0.1473 0.2036 0.1781Again, no significant difference in discrimination is observed.
Finally, we can plot the estimated time-dependent AUCs:
plotAUCcurve(ROC.thick, FP = 2, add = FALSE,conf.int = FALSE,conf.band = FALSE, col = "black")
plotAUCcurve(ROC.cox, FP = 2, add =TRUE, conf.int = FALSE,conf.band = FALSE, col = "blue")
plotAUCcurve(ROC.cox2, FP = 2, add =TRUE, conf.int = FALSE,conf.band = FALSE, col = "red")
title(main="Time-dependent AUCs")
legend("topleft", legend=c("Tumor thickness","Thick + sex", "Thick + sex+ulc"), lty=c(1,1,1),lwd=c(2,2,2),
col=c("black","blue", "red"), bty="n")
We can also fix a specific time-horizon :
plot(ROC.thick,time=2200,add = FALSE, col = "black",title=FALSE, lwd=2)
plot(ROC.cox,time=2200,add = TRUE, col = "blue",title=FALSE, lwd=2)
plot(ROC.cox2,time=2200,add = TRUE, col = "red",title=FALSE, lwd=2)
legend("bottomright", legend=c("Tumor thickness","Thick + sex", "Thick + sex + ulc"), lty=c(1,1,1),lwd=c(2,2,2),
col=c("black","blue", "red"), bty="n")
title(main="Time-dependent ROC curves at 2200 days")
Note that discrimination in this context is focused exclusively on the performance of the linear predictor extracted from the Cox model. We are not considering here the survival estimates.
Calibration for the Cox model
Just for didactic reasons, we are in this example IGNORING the competing risk present in these data, and act as we have an unique event. Therefore the method is illustrated is valid only in the standard survival setting, not in the competing risk setting (more advanced approaches are needed in that case, as for example estimating cumulative incidence curves using the multi-state approach). We now evaluate calibration for the Cox model: first of all, we determine summaries of variables for effect and plotting ranges, values to adjust to, and overall ranges for the dataset Melanom, since this is required in order to use the R functions developed in the rms library.
options(datadist='dd')
dd <-datadist(melanom)
dd## no status days ulc thick sex
## Low:effect 222.0000 1 1525.0000 1 97.0000 1
## Adjust to 469.0000 2 2005.0000 1 194.0000 1
## High:effect 731.0000 3 3042.0000 2 356.0000 2
## Low:prediction 11.0000 1 294.2195 1 32.0000 1
## High:prediction 943.0488 3 4119.2439 2 838.7805 2
## Low 2.0000 1 10.0000 1 10.0000 1
## High 992.0000 3 5565.0000 2 1742.0000 2
##
## Values:
##
## status : 1 2 3
## ulc : 1 2
## sex : 1 2Then, we re-estimate the Cox model by adding some options useful to fix the time horizon to evaluate calibration, in this case fixed at 2000 days:
mod.cov2 <- cph(Surv(days,status==1)~sex+log(thick)+ulc, data=melanom,x=TRUE , y=TRUE,surv=TRUE,time.inc=2000)Finally, we use the function calibrate by comparing predicted probabilities from the Cox model with the observed ones estimated by the Kaplan-Meier:
calkm <- calibrate(mod.cov2, u=2000, m=40, cmethod= 'KM', B=200)## Using Cox survival estimates at 2000 Dayscalkm## calibrate.cph(fit = mod.cov2, cmethod = "KM", u = 2000, m = 40,
## B = 200)
##
##
## n=205 B=200 u=2000 Day
##
## index.orig training test mean.optimism mean.corrected n
## [1,] 0.01505666 0.005391453 0.026202605 -0.0208111521 0.035867808 194
## [2,] 0.01466150 0.013204855 0.003653485 0.0095513694 0.005110134 200
## [3,] 0.01195721 0.004408990 0.001133622 0.0032753681 0.008681838 200
## [4,] -0.01622951 0.009643115 -0.009066033 0.0187091476 -0.034938658 200
## [5,] 0.01160306 0.005764159 0.005642487 0.0001216718 0.011481386 199
## mean.predicted KM KM.corrected std.err
## [1,] 0.4394888 0.4545455 0.4753566 0.17797217
## [2,] 0.6603385 0.6750000 0.6654486 0.10971343
## [3,] 0.8233362 0.8352934 0.8320181 0.07384917
## [4,] 0.8978469 0.8816174 0.8629083 0.06477229
## [5,] 0.9396165 0.9512195 0.9510978 0.03536639plot(calkm, cex.lab=0.7)
Another more general approach to estimating calibration curves for
survival models is to use the flexible adaptive hazard regression
approach (not covered here). For survival models, “predicted” means
predicted survival probability at a single time point, and “observed”
refers to the corresponding Kaplan-Meier survival estimate. The number
of intervals of predicted survival probabilities in which subjects are
ranked depend on the parameter m that correspond to the average
number of subjects in each group evaluated at u-time.
As an internal validation tool, the bootstrap is used to de-bias the
estimates to correct for overfitting, allowing estimation of the
likely future calibration performance of the fitted model.