1. Hands-on excerise building a predictive model for time-to-event data using COX PH model
1.1 Descriptive analysis
For this tutorial, we will be using the HBV data. This data set contains information on the follow-up of untreated hepatitis B patients (natural disease course) with advanced fibrosis/cirrhosis.
The time-to-event outcome is defined using HCC and HCC_T,
predictors <-colnames(hbv2)[c(2:11, 14:15)]#making formulasuniv_formulas <-sapply(predictors,function(x)as.formula(paste('Surv(HCC_t,HCC)~',x)))#making a list of models;univ_models <-lapply(univ_formulas,function(x){coxph(x,data=hbv2)})#extract summary statistics;univ_results <-lapply(univ_models,function(x){return(summary(x)$coef)})univ_results
based on the univariate analysis the following predictors are statistically significant: age (and agec), spleno, bilirub, plat (and platc), albumin, gamma.
perform variable selection using LASSO
library(glmnet)y <-Surv(hbv2$HCC_t,hbv2$HCC)x <-model.matrix( ~ AGE + SEX + HBEAG + SPLENO + BILIRUB + AST + ALT + PLAT + ALBUMIN + GAMMA, hbv2)set.seed(123)cv.fit <-cv.glmnet(x, y, family="cox", type.measure ="C")plot(cv.fit)
coef(cv.fit, s ="lambda.min")
11 x 1 sparse Matrix of class "dgCMatrix"
1
(Intercept) .
AGE 0.027952059
SEX1 -0.006688095
HBEAG1 0.475707548
SPLENO1 0.294741080
BILIRUB 0.011933437
AST .
ALT -0.028187790
PLAT -0.007172029
ALBUMIN -0.081434515
GAMMA 0.020888938
predictor selection kept the following variables: age, sex, HBEAG, SPLENO, BILIRUB, ALT, PLAT, ALBUMIN, GAMMA
We do have better C-statistics under model 2. The Toronto HCC risk index is developed using categorized lab values.
We are also interested in building a risk score in this example. Given that logging the lab values only improved the model very slightly, we will stay with model m1 for interpretability and simplicity.
1.4 Model predictive performance
1.4.1 Regression on PI in the validation data
The prognostic index (PI): in cox model is the effect of the covariates and \(log(h(t)) = log(h_0(t)) = exp(PI)\).
#1. obtain model estimated log hazard function;logh <-model.matrix(m1) %*%coef(m1)#2. centre the PI by using the mean;m1$means #returns mean value of all covaraites;
AGE factor(SEX)1 factor(HBEAG)1 factor(SPLENO)1 BILIRUB
46.577520 0.000000 0.000000 0.000000 16.227115
ALT PLAT ALBUMIN GAMMA
3.827273 143.653846 39.487179 18.837179
we can use the distribution of PI to define risk groups, e.g., by quartiles.
when we have a validation dataset, we can fit COX PH model on PI in the validated data set to assess model discrimination and model fit.
since we don’t have a validation dataset, to demonstrate PI, I will use a mock validation dataset generated via simple random sample of the original data.
# generating mock validation data;set.seed(123)sample_id <-sample(x =1:length(hbv2$SubjectID), size=100, replace = T)hbv2_val <- hbv2[sample_id,]# fit model in validation data;m1_val <-coxph(Surv(HCC_ty, HCC) ~ AGE+factor(SEX)+factor(HBEAG)+factor(SPLENO)+BILIRUB+ALT+PLAT+ALBUMIN+GAMMA, data=hbv2_val)#1. obtain model estimated log hazard function(updated!);#using coef estimated from m1 using the original data;logh_val <-model.matrix(m1_val) %*%coef(m1)logh_mean_val <-sum(coef(m1)*m1_val$means)PI_val <- logh_val - logh_mean_val # centred PI;hbv2_val$PI_val <- PI_valfit_val <-coxph(Surv(HCC_ty, HCC) ~ PI_val +offset(PI_val), data=hbv2_val)summary(fit_val)
Call:
coxph(formula = Surv(HCC_ty, HCC) ~ PI_val + offset(PI_val),
data = hbv2_val)
n= 100, number of events= 34
coef exp(coef) se(coef) z Pr(>|z|)
PI_val -0.02858 0.97183 0.17885 -0.16 0.873
exp(coef) exp(-coef) lower .95 upper .95
PI_val 0.9718 1.029 0.6845 1.38
Concordance= 0.753 (se = 0.041 )
Likelihood ratio test= 0.03 on 1 df, p=0.9
Wald test = 0.03 on 1 df, p=0.9
Score (logrank) test = 0.03 on 1 df, p=0.9
perfect fit in the validation data (this is expected since the validation data is drawn from the original data)! the estimated slope is 0.9718 (se 0.18), p-value = 0.873
1.4.2 Measures of discrimination using C-stat
Harrell’s C quantifies the degree of concordance as the proportion of evaluate pairs where the patient with a longer survival time has better predicted survival;
Plot discrimination measures over time (dynamic AUC)
library(survivalROC)## Define a functionfun.survivalROC <-function(lp, t) { res <-with(hbv2,survivalROC(Stime = HCC_ty,status = HCC,marker =get(lp),predict.time = t,method ="KM")) # KM method without smoothing;## Plot ROCswith(res, plot(TP ~ FP, type ="l", main =sprintf("Y = %.0f, C = %.2f", t, AUC)))abline(a =0, b =1, lty =2) res}## 2 x 5 layoutlayout(matrix(1:10, byrow = T, ncol =5))res.survivalROC.m1 <-lapply(1:10, function(t) {fun.survivalROC(lp ="PI", t)})
# Save PI;lp_dev <-predict(m1, newdata = hbv2)mtime <-seq(0.5, 15, by=0.5) #6 month to 15 years by 6 month;cstat <-matrix(0, length(mtime), 2) for (i in1:length(mtime)) {# Harrell C c <-concordance(m1, ymax = mtime[i]) # Save resultscstat[i,] <-c(coef(c), sqrt(c(vcov(c))))}colnames(cstat) <-c("Harrell_C","Harrell_C_SE")cstat <-data.frame(cstat)alpha <- .05plot(mtime, cstat[, c("Harrell_C")], type ='l', col =1, ylim =c(0, 1),xlim =c(0, 15),lty =1, lwd =3,xlab ="Years", ylab ="Concordance",bty ="n")polygon(c(mtime, rev(mtime)),c(cstat[, c("Harrell_C")] -qnorm(1- alpha /2) * cstat[, c("Harrell_C_SE")],rev(cstat[, c("Harrell_C")] +qnorm(1- alpha /2) * cstat[, c("Harrell_C_SE")])),col =rgb(160, 160, 160, maxColorValue =255, alpha =100),border =FALSE)title("Discrimination over the time \n Harrell C-statistics",cex.main =0.95)
1.4.3 Assessing Calibration
Two overall performance measures are proposed for prediction models with a survival outcome:
Brier score: it is the mean squared difference between observed event indicators and predicted risks at a fixed time point (e.g. at 5 years), lower is better;
library(riskRegression)
riskRegression version 2023.09.08
Score(list("Model 1"= m1),formula =Surv(HCC_ty, HCC) ~1, data = hbv2, conf.int =TRUE, times =4.99,cens.model ="cox",metrics ="brier")
Metric Brier:
Results by model:
model times Brier lower upper
1: Null model 4.99 18.3 14.6 22.0
2: Model 1 4.99 15.0 11.3 18.7
Results of model comparisons:
times model reference delta.Brier lower upper p
1: 4.99 Model 1 Null model -3.3 -6.3 -0.3 0.03060707
NOTE: Values are multiplied by 100 and given in %.
