Tutorial 2: Validate the prediction model

Overview

In Tutorial 1 you developed a prediction model. A model always performs better on the data it was developed on than on new patients, so in this tutorial you will check how well it really works. It has three parts:

1. Uniform shrinkage and optimism-corrected C-index: use the bootstrap to calculate the uniform shrinkage factor and calculate optimism.

2. Model performance: calculate predicted risks, make calibration plots and a decision curve.

3. Internal-external cross-validation: develop the model on some centers and test it on the center that was left out.

This tutorial can be run on its own: if you have not run Tutorial 1, the objects it needs are rebuilt automatically (see Load the data and models).

Download the tutorial

Do not just read this page: run the code yourself. You only learn R by doing it. Download the R Markdown file of this tutorial (Tutorial_2_Validate_Prediction_Model.Rmd) by clicking the Code button at the top right of this page and choosing Download Rmd, and run every chunk in RStudio, one by one. Read the output, change the code, break it and fix it, and try the exercises before you open the answers. Reading the results on this page is not a substitute.

To make sure everything works:

1. Create a folder for the course, for example Predictive_models/, and put the downloaded .Rmd file in it. Tutorial 1 is not required.

2. Open the .Rmd file in RStudio. When you knit or run chunks from an Rmd file, R uses the folder that contains the file as the working directory. Check this by running:

   getwd()

   The output should be the path of the folder from step 1. If it is not, set it with setwd("path/to/your/folder") or in RStudio via Session > Set Working Directory > To Source File Location.

3. Do not use file paths in the code. All files are read from and saved to the working directory, so keep all files in this one folder.

4. Running the tutorial saves its figures (calibration plots, decision curve and the internal-external cross-validation) in this same folder.

The folder will then look like this:

Predictive_models/
├── Tutorial_1_Develop_Prediction_Model.Rmd
├── Tutorial_2_Validate_Prediction_Model.Rmd
├── model_performance.png           (created by Tutorial 2)
├── dca.png                         (created by Tutorial 2)
└── Internal_external_cross_validation.png  (created by Tutorial 2)

Setup

Clear the environment and load libraries

As in Tutorial 1, everything is saved in one folder: the folder that contains this Rmd file (the working directory), so no file paths are needed. If you also do Tutorial 1, keep both tutorials in the same folder.

A package is a collection of functions that someone wrote and shared. library(package) loads a package, so that you can use its functions directly. We only load rms, because we use some of its functions (such as cph, validate and fastbw) without the package:: prefix. All other packages are called with the :: notation (see below), so they do not have to be loaded, but they have to be installed once. We use survival, PredictionTools (calibration plots), dcurves (decision curve analysis), ggplot2 and ggpubr (making and combining figures), and dplyr, mice and Hmisc (only to rebuild the models below).

# clear the environment: remove all objects that R currently remembers
rm(list=ls())

# load the rms package: after this, all its functions can be used directly
library(rms)

# The other packages must be installed once on your computer (not run here):
# install.packages(c("survival", "dcurves", "ggplot2", "ggpubr", "dplyr", "mice", "Hmisc"))
# remotes::install_github("CHMMaas/PredictionTools")  # documentation: https://github.com/CHMMaas/PredictionTools

The :: notation

In this tutorial you will often see code like rms::validate(...). The notation package::function() means “use the function function from the package package”.

You do not have to write it this way once you have loaded a package with library(): after library(rms), validate(...) works as well. We still write it out because it is nice to see where each function comes from, and it avoids mistakes when two packages have a function with the same name (R then uses the one from the package that was loaded last).

Set the seed

We set the random seed once, here at the top. Every random step below (the bootstrap samples) then continues from this starting point, so the whole tutorial is reproducible without resetting the seed before each call. If you change the order of the code or add extra random steps, the numbers will change.

# set the seed once for the whole tutorial
set.seed(1)

Load the data and models

This tutorial needs the data, the imputed data sets and the models from Tutorial 1. The code below rebuilds them (a condensed copy of Tutorial 1, which takes a minute), so you can start with this tutorial directly.

Show the code that rebuilds the Tutorial 1 objects

# data: pbc from the survival package, with an artificial center variable
data(pbc, package="survival")
data <- pbc
data$center <- rep(1:2, length.out=nrow(data))
data <- dplyr::select(data, time, status,
               trt, age, sex, ascites, hepato, spiders, edema,
               bili, chol, albumin, copper, alk.phos, ast, trig,
               platelet, protime, stage, center)
data$trt     <- factor(data$trt, levels=c(1, 2),
                       labels=c("D-penicillamine", "Placebo"))
data$sex     <- factor(data$sex, levels=c("m", "f"),
                       labels=c("Male", "Female"))
data$ascites <- factor(data$ascites, levels=c(0, 1),
                       labels=c("No", "Yes"))
data$hepato  <- factor(data$hepato, levels=c(0, 1),
                       labels=c("No", "Yes"))
data$spiders <- factor(data$spiders, levels=c(0, 1),
                       labels=c("No", "Yes"))
data$edema   <- factor(data$edema, levels=c(0, 0.5, 1),
                       labels=c("No edema", "Untreated or successfully treated",
                                "Despite diuretics"))
data$stage   <- factor(data$stage, levels=1:4,
                       labels=c("1", "2", "3", "4"))
data$center  <- factor(data$center, levels=c(1, 2),
                       labels=c("Center 1", "Center 2"))

# outcome: death within 10 years, transplant is censored
S <- survival::Surv(time=data$time/365, event=as.numeric(data$status==2))
horizon <- 10
S_10 <- S
S_10[S_10[, 1]>horizon, 2] <- 0
S_10[S_10[, 1]>horizon, 1] <- horizon

# multiple imputation with event indicator and Nelson-Aalen estimate
data_mi <- data
data_mi$time <- NULL
data_mi$status <- NULL
data_mi$event_10 <- S_10[, "status"]
data_mi$cum_hazard <- mice::nelsonaalen(data.frame(time=S_10[, "time"],
                                                   status=S_10[, "status"]),
                                        time, status)
m <- 5
imputed_data <- mice::mice(data_mi, m=m, seed=1, print=FALSE)

# models
dd <- rms::datadist(data)
options(datadist="dd")
fit_args <- list(x=TRUE, y=TRUE, surv=TRUE)
full_non_linear_model <- Hmisc::fit.mult.impute(
  S_10 ~ rcs(age, 3) + rcs(bili, 3) + rcs(albumin, 3) + sex + ascites +
    hepato + stage, rms::cph,
  fitargs=fit_args, xtrans=imputed_data, data=data, fit.reps=TRUE)

## 
## Wald Statistic Information
## 
## Variance Inflation Factors Due to Imputation:
## 
##         age        age'        bili       bili'     albumin    albumin' 
##        1.01        1.01        1.01        1.00        1.00        1.01 
##  sex=Female ascites=Yes  hepato=Yes     stage=2     stage=3     stage=4 
##        1.00        1.15        1.09        1.09        1.08        1.09 
## 
## Fraction of Missing Information:
## 
##         age        age'        bili       bili'     albumin    albumin' 
##        0.01        0.01        0.01        0.00        0.00        0.01 
##  sex=Female ascites=Yes  hepato=Yes     stage=2     stage=3     stage=4 
##        0.00        0.13        0.08        0.08        0.07        0.08 
## 
## d.f. for t-distribution for Tests of Single Coefficients:
## 
##         age        age'        bili       bili'     albumin    albumin' 
##    70992.09    29693.89   106704.51   199948.55   205552.84    32084.48 
##  sex=Female ascites=Yes  hepato=Yes     stage=2     stage=3     stage=4 
##   358897.26      242.90      630.61      651.50      799.20      602.05 
## 
## The following fit components were averaged over the 5 model fits:
## 
##   linear.predictors means stats center

form_full <- S_10 ~ rcs(age, 3) + rcs(bili, 3) + albumin + sex +
  ascites + hepato + stage
full_model <- Hmisc::fit.mult.impute(form_full, rms::cph,
  fitargs=fit_args, xtrans=imputed_data, data=data, fit.reps=TRUE)

## 
## Wald Statistic Information
## 
## Variance Inflation Factors Due to Imputation:
## 
##         age        age'        bili       bili'     albumin  sex=Female 
##        1.01        1.01        1.01        1.00        1.00        1.00 
## ascites=Yes  hepato=Yes     stage=2     stage=3     stage=4 
##        1.14        1.10        1.08        1.08        1.09 
## 
## Fraction of Missing Information:
## 
##         age        age'        bili       bili'     albumin  sex=Female 
##        0.01        0.01        0.01        0.00        0.00        0.00 
## ascites=Yes  hepato=Yes     stage=2     stage=3     stage=4 
##        0.13        0.09        0.08        0.07        0.08 
## 
## d.f. for t-distribution for Tests of Single Coefficients:
## 
##         age        age'        bili       bili'     albumin  sex=Female 
##    66673.51    25006.17   103533.09   203798.25 31013743.69   337198.38 
## ascites=Yes  hepato=Yes     stage=2     stage=3     stage=4 
##      251.90      492.12      655.56      784.71      576.02 
## 
## The following fit components were averaged over the 5 model fits:
## 
##   linear.predictors means stats center

form_bw <- S_10 ~ rcs(age, 3) + rcs(bili, 3) + albumin + stage
bw_model <- Hmisc::fit.mult.impute(form_bw, rms::cph,
  fitargs=fit_args, xtrans=imputed_data, data=data, fit.reps=TRUE)

## 
## Wald Statistic Information
## 
## Variance Inflation Factors Due to Imputation:
## 
##     age    age'    bili   bili' albumin stage=2 stage=3 stage=4 
##    1.00    1.00    1.00    1.00    1.01    1.09    1.06    1.07 
## 
## Fraction of Missing Information:
## 
##     age    age'    bili   bili' albumin stage=2 stage=3 stage=4 
##    0.00    0.00    0.00    0.00    0.01    0.08    0.06    0.06 
## 
## d.f. for t-distribution for Tests of Single Coefficients:
## 
##          age         age'         bili        bili'      albumin      stage=2 
##  20946309.50 131064158.78    381202.54    364582.19     97205.48       592.59 
##      stage=3      stage=4 
##      1199.64      1069.91 
## 
## The following fit components were averaged over the 5 model fits:
## 
##   linear.predictors means stats center

# set data distributions
dd <- rms::datadist(data)
options(datadist="dd")

1 Uniform shrinkage and optimism-corrected C-index

A model fitted on a data set always performs a bit better on that same data set than on new patients. This is called optimism or overfitting, and it is worse when the model was built using variable selection. In this part you use bootstrapping to:

1. estimate a uniform shrinkage factor that shrinks the coefficients towards zero, which makes predictions less extreme and better calibrated on new patients, and
2. estimate the optimism-corrected C-index, a measure of how well the model separates patients who die early from patients who survive.

The bootstrap repeats the entire backward selection procedure in every bootstrap sample, so the optimism includes the effect of variable selection.

1.1 Bootstrap validation on every imputed data set

For each imputed data set we validate the full model with 200 bootstrap samples. In each bootstrap sample validate() repeats backward selection (bw=TRUE, p < 0.05), refits the model, and tests it on the original data. Three quantities are stored:

• the calibration slope in the test data (Slope, column test), which is our shrinkage factor,
• the apparent Dxy of the final model on the development data (Dxy, column index.orig), and
• the optimism in Dxy (Dxy, column optimism).

Dxy is converted to the C-index with C = Dxy/2 + 0.5, so the optimism in the C-index is the optimism in Dxy divided by 2.

# calculate uniform shrinkage factor across all imputations
# start with empty vectors, c() with nothing in it, to collect the results
slope <- c()
apparent_C <- c()
opt_C <- c()

# the for loop repeats the code between { } for every imputed data set i = 1, ..., m
for (i in 1:m){
  # validate() arguments:
  #   full_model$fits[[i]]: the model fitted on the i-th imputed data set
  #   B=200: the number of bootstrap samples
  #   bw=TRUE, rule="p", sls=0.05, type="individual": repeat the backward
  #     selection (p < 0.05) in every bootstrap sample
  #   tol: how strict R is about nearly identical predictors (here: less strict)
  val_model <- rms::validate(full_model$fits[[i]],
                        B=200,
                        bw=TRUE,
                        rule="p",
                        sls=0.05,
                        type="individual", tol=1e-13)

  # val_model["Slope", "test"] takes the cell in row "Slope" and column "test";
  # c(old, new) adds the new value to the end of the vector
  slope <- c(slope, val_model["Slope", "test"])
  # Dxy is converted to the C-index: C = Dxy/2 + 0.5
  apparent_C <- c(apparent_C, val_model["Dxy", "index.orig"]/2 + 0.5)
  opt_C <- c(opt_C, val_model["Dxy", "optimism"]/2)
}

## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  sex     0.11   1    0.7352 0.11     1    0.7352 -1.89
##  hepato  0.57   1    0.4521 0.68     2    0.7119 -3.32
##  ascites 2.16   1    0.1416 2.84     3    0.4170 -3.16
## 
## Approximate Estimates after Deleting Factors
## 
##             Coef    S.E.  Wald Z         P
## age      0.06839 0.02381  2.8723 4.074e-03
## age'    -0.03714 0.02236 -1.6610 9.672e-02
## bili     0.80138 0.10722  7.4745 7.749e-14
## bili'   -2.33364 0.37180 -6.2767 3.459e-10
## albumin -0.83717 0.19953 -4.1956 2.721e-05
## stage=2  0.37999 0.74155  0.5124 6.084e-01
## stage=3  0.79708 0.72348  1.1017 2.706e-01
## stage=4  1.25171 0.72724  1.7212 8.522e-02
## 
## Factors in Final Model
## 
## [1] age     bili    albumin stage

## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  hepato  0.02   1    0.9005 0.02     1    0.9005 -1.98
##  sex     0.19   1    0.6648 0.20     2    0.9033 -3.80
##  ascites 2.89   1    0.0891 3.09     3    0.3774 -2.91
## 
## Approximate Estimates after Deleting Factors
## 
##             Coef    S.E.  Wald Z         P
## age      0.06744 0.02387  2.8257 4.717e-03
## age'    -0.03699 0.02244 -1.6490 9.916e-02
## bili     0.79170 0.10823  7.3147 2.580e-13
## bili'   -2.29741 0.37576 -6.1140 9.717e-10
## albumin -0.82819 0.20013 -4.1383 3.499e-05
## stage=2  0.31366 0.74301  0.4222 6.729e-01
## stage=3  0.78342 0.72401  1.0821 2.792e-01
## stage=4  1.29211 0.72657  1.7784 7.534e-02
## 
## Factors in Final Model
## 
## [1] age     bili    albumin stage

## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  hepato  0.01   1    0.9156 0.01     1    0.9156 -1.99
##  sex     0.09   1    0.7694 0.10     2    0.9526 -3.90
## 
## Approximate Estimates after Deleting Factors
## 
##                 Coef    S.E.  Wald Z         P
## age          0.07121 0.02401  2.9652 3.024e-03
## age'        -0.04477 0.02273 -1.9701 4.882e-02
## bili         0.79408 0.10861  7.3113 2.647e-13
## bili'       -2.34493 0.37721 -6.2165 5.084e-10
## albumin     -0.72450 0.20710 -3.4983 4.682e-04
## ascites=Yes  0.58661 0.23893  2.4552 1.408e-02
## stage=2      0.34249 0.74112  0.4621 6.440e-01
## stage=3      0.83950 0.72398  1.1596 2.462e-01
## stage=4      1.21614 0.72765  1.6713 9.466e-02
## 
## Factors in Final Model
## 
## [1] age     bili    albumin ascites stage

## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  sex     0.10   1    0.7493 0.10     1    0.7493 -1.90
##  hepato  0.46   1    0.4971 0.56     2    0.7546 -3.44
##  ascites 3.23   1    0.0721 3.80     3    0.2842 -2.20
## 
## Approximate Estimates after Deleting Factors
## 
##             Coef    S.E.  Wald Z         P
## age      0.06761 0.02398  2.8190 4.817e-03
## age'    -0.03617 0.02249 -1.6079 1.079e-01
## bili     0.78418 0.10885  7.2039 5.850e-13
## bili'   -2.27420 0.37683 -6.0350 1.589e-09
## albumin -0.87070 0.20159 -4.3192 1.566e-05
## stage=2 -0.07336 0.62166 -0.1180 9.061e-01
## stage=3  0.43475 0.59812  0.7269 4.673e-01
## stage=4  0.91167 0.60111  1.5166 1.294e-01
## 
## Factors in Final Model
## 
## [1] age     bili    albumin stage

## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  sex     0.09   1    0.7581 0.09     1    0.7581 -1.91
##  hepato  0.11   1    0.7349 0.21     2    0.9005 -3.79
##  ascites 3.17   1    0.0751 3.38     3    0.3371 -2.62
## 
## Approximate Estimates after Deleting Factors
## 
##             Coef    S.E.  Wald Z         P
## age      0.06739 0.02384  2.8261 4.712e-03
## age'    -0.03690 0.02235 -1.6510 9.873e-02
## bili     0.78875 0.10923  7.2208 5.167e-13
## bili'   -2.28764 0.37968 -6.0252 1.689e-09
## albumin -0.84314 0.20032 -4.2090 2.565e-05
## stage=2  0.33196 0.74363  0.4464 6.553e-01
## stage=3  0.77505 0.72399  1.0705 2.844e-01
## stage=4  1.29600 0.72649  1.7839 7.444e-02
## 
## Factors in Final Model
## 
## [1] age     bili    albumin stage

1.2 Shrinkage factor

The shrinkage factor is the mean slope across imputations. A value below 1 indicates overfitting: the model’s predictions are too extreme, and multiplying the coefficients by the factor corrects this.

First look at the slope of every imputed data set. If they differ a lot, a single average would hide that:

# slope of each imputed data set
round(slope, 3)

## [1] 0.905 0.916 0.912 0.913 0.896

# how much do they vary around the mean?
round(c(mean=mean(slope), sd=sd(slope), min=min(slope), max=max(slope)), 3)

##  mean    sd   min   max 
## 0.908 0.008 0.896 0.916

The variation between imputations comes from the different imputed values and from the random bootstrap samples. Here the slopes are close to their mean, so averaging over the imputations gives one sensible shrinkage factor:

# shrinkage factor: the mean of the slopes over the imputed data sets
shrinkage_factor <- mean(slope)
# cat() prints text and values; "\n" starts a new line
cat("Shrinkage factor:", shrinkage_factor, "\n")

## Shrinkage factor: 0.9083665

# shrunk coefficients: multiply every coefficient of the final model by the factor
shrinkage_factor*bw_model$coefficients

##         age        age'        bili       bili'     albumin     stage=2 
##  0.06139414 -0.03347546  0.71964386 -2.09131974 -0.76150109  0.24212356 
##     stage=3     stage=4 
##  0.65159136  1.08803671

1.3 Optimism-corrected C-index

The C-index (also called the c-statistic) measures discrimination: how well the model distinguishes between low- and high-risk patients.

• 0.5 means that the model performs no better than random guessing in distinguishing between low- and high-risk patients.
• 1 means that the model perfectly distinguishes between low- and high-risk patients.

The validate() object already contains the C-index of the final model, so we do not have to calculate it ourselves. Here is the output for the last imputed data set (val_model from the loop above):

# output of validate() for the last imputed data set
val_model

##       index.orig training   test optimism index.corrected   n
## Dxy       0.6711   0.6868 0.6665   0.0203          0.6507 200
## R2        0.4309   0.4508 0.4133   0.0375          0.3934 200
## Slope     1.0000   1.0000 0.8962   0.1038          0.8962 200
## D         0.1338   0.1430 0.1266   0.0164          0.1174 200
## U        -0.0012  -0.0012 0.0039  -0.0051          0.0039 200
## Q         0.1349   0.1442 0.1227   0.0215          0.1134 200
## g         1.6010   1.7179 1.5335   0.1844          1.4166 200
## 
## Factors Retained in Backwards Elimination
## 
##  age bili albumin sex ascites hepato stage
##  *   *    *                          *    
##  *   *    *       *                  *    
##  *   *    *       *   *              *    
##  *   *    *                          *    
##  *   *    *           *              *    
##  *   *                *              *    
##  *   *    *           *              *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *           *       *      *    
##  *   *                *              *    
##  *   *    *           *       *      *    
##  *   *                *              *    
##  *   *    *       *                  *    
##  *   *    *                          *    
##  *   *                *              *    
##  *   *    *           *              *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                   *           
##  *   *    *           *              *    
##  *   *    *                          *    
##  *   *                *              *    
##  *   *    *       *   *              *    
##  *   *    *           *                   
##  *   *    *                          *    
##  *   *    *       *                  *    
##  *   *    *           *              *    
##  *   *    *                          *    
##  *   *    *           *       *      *    
##  *   *    *           *              *    
##  *   *    *                          *    
##  *   *                               *    
##  *   *    *           *              *    
##  *   *    *                          *    
##  *   *    *           *              *    
##  *   *    *                          *    
##  *   *    *           *              *    
##  *   *    *                          *    
##  *   *    *           *                   
##  *   *    *       *   *                   
##  *   *    *                          *    
##  *   *    *           *              *    
##  *   *    *       *   *              *    
##  *   *    *           *       *           
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *           *                   
##  *   *    *           *              *    
##  *   *    *           *              *    
##  *   *    *       *   *              *    
##  *   *                *              *    
##  *   *    *           *       *      *    
##  *   *    *           *              *    
##  *   *    *           *       *           
##  *   *    *           *              *    
##  *   *    *           *       *           
##  *   *    *                          *    
##  *   *    *           *              *    
##  *   *    *           *              *    
##  *   *    *           *              *    
##  *   *    *                          *    
##      *    *           *                   
##  *   *    *           *              *    
##  *   *    *           *              *    
##  *   *    *           *              *    
##  *   *    *                          *    
##  *   *    *           *              *    
##  *   *    *           *              *    
##  *   *    *           *              *    
##  *   *    *           *              *    
##  *   *    *       *                  *    
##  *   *    *                          *    
##  *   *    *                   *      *    
##  *   *    *           *              *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *                *       *      *    
##  *   *    *           *                   
##  *   *    *           *              *    
##  *   *    *                          *    
##  *   *                *              *    
##  *   *    *       *                  *    
##  *   *                               *    
##  *   *    *       *   *              *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *       *   *              *    
##  *   *    *                               
##  *   *    *                          *    
##  *   *    *       *   *                   
##  *   *            *   *              *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *           *              *    
##  *   *    *           *                   
##  *   *    *           *              *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *           *                   
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *                *              *    
##  *   *    *           *              *    
##      *    *       *   *              *    
##  *   *    *           *              *    
##  *   *    *           *              *    
##      *    *           *              *    
##  *   *    *                          *    
##  *   *    *                          *    
##      *    *           *              *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *       *   *              *    
##  *   *    *                   *      *    
##  *   *                *              *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *           *              *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *       *   *              *    
##  *   *    *       *   *                   
##  *   *    *       *                  *    
##  *   *    *                          *    
##  *   *    *           *              *    
##  *   *    *           *              *    
##  *   *    *       *   *              *    
##  *   *    *           *              *    
##  *   *    *       *   *              *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                   *      *    
##  *   *    *           *              *    
##  *   *    *           *       *           
##  *   *    *       *                  *    
##  *   *    *           *              *    
##  *   *    *                               
##  *   *    *                          *    
##  *   *    *           *       *           
##  *   *    *                          *    
##  *   *                               *    
##  *   *    *                          *    
##  *   *    *                          *    
##      *    *                          *    
##  *   *    *           *                   
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *           *              *    
##  *   *    *                          *    
##  *   *    *           *              *    
##  *   *    *       *                  *    
##  *   *    *           *              *    
##  *   *    *       *   *              *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *           *              *    
##  *   *    *           *              *    
##  *   *                               *    
##  *   *    *           *              *    
##      *    *           *              *    
##  *   *    *                               
##  *   *    *                          *    
##  *   *    *           *              *    
##  *   *    *           *              *    
##  *   *    *           *              *    
##  *   *    *           *              *    
##  *   *    *                          *    
##  *   *    *       *   *              *    
##  *   *    *           *              *    
##      *    *       *   *              *    
##  *   *    *       *                  *    
##  *   *                        *      *    
##  *   *    *           *       *      *    
##  *   *    *       *                  *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *                          *    
##  *   *    *           *              *    
##  *   *    *                          *    
##  *   *    *           *              *    
##  *   *    *           *              *    
##  *   *                *              *    
##  *   *    *                               
## 
## Frequencies of Numbers of Factors Retained
## 
##  3  4  5  6 
## 10 97 76 17

Each row is a performance measure (Dxy is the C-index on another scale: C = Dxy/2 + 0.5, Slope is the calibration slope). The columns are:

• index.orig: the apparent performance of the final model on the development data.
• training: the average performance of the model developed in each bootstrap sample, evaluated on that same bootstrap sample.
• test: the average performance of the model developed in each bootstrap sample, evaluated on the original data.
• optimism: training minus test, how much the apparent performance is too optimistic.
• index.corrected: index.orig minus optimism, the optimism-corrected performance.
• n: the number of bootstrap samples that were used.

For the C-index we take index.orig and the optimism of Dxy from every imputed data set (stored in apparent_C and opt_C above), average them over the imputations and subtract the mean optimism:

# optimism-corrected C-index: mean(x) averages a vector over the imputed data sets
optimism_C <- mean(opt_C)
original_C <- mean(apparent_C)
cat("Original C-index:", round(original_C, 3), "\n",
    "Optimism:", round(optimism_C, 3), "\n",
    "Optimism-corrected C-index:", round(original_C - optimism_C, 3), "\n")

## Original C-index: 0.836 
##  Optimism: 0.009 
##  Optimism-corrected C-index: 0.827

Note on the confidence interval. validate() only returns the averages over the bootstrap samples, not a variance or the C-index of every bootstrap sample, so it does not give a confidence interval. If you want one, you have to code the bootstrap yourself: repeat the whole procedure (imputation, variable selection, fitting, validation) in bootstrap samples and take the percentiles of the optimism-corrected C-index.

1.4 Exercises: Part 1

Try each exercise yourself first. Then open the field below it to see the answer, the code and its results. The answers use new object names. Where they use random numbers, they restore the state of the random number generator afterwards, so the rest of the tutorial gives the same numbers with or without them.

Exercise 1. Compare the original C-index (original_C) with the optimism-corrected one. How large is the optimism, and is this a lot for this model?

Show answer, code and results

# original C-index, optimism and corrected C-index
round(c(original=original_C,
        optimism=optimism_C,
        corrected=original_C - optimism_C), 3)

##  original  optimism corrected 
##     0.836     0.009     0.827

# how much information do we have compared with the number of parameters?
n_events <- sum(S_10[, "status"])
n_param <- length(full_model$coefficients)
c(events=n_events, parameters=n_param, events_per_parameter=round(n_events/n_param, 1))

##               events           parameters events_per_parameter 
##                156.0                 11.0                 14.2

The optimism is the difference between the original and the corrected C-index. Optimism is larger when there are few events compared with the number of candidate parameters, so the last line gives a feeling for how much overfitting to expect. Whether the optimism is “a lot” depends on how it compares with the C-index you need for your purpose.

Exercise 2. Change B=200 to B=50 and back to B=500. How stable is the shrinkage factor?

Show answer, code and results

To keep this quick we use only the first imputed data set.

# save the state of the random number generator, so that this exercise does not
# change the random numbers used in the rest of the tutorial
old_seed <- .Random.seed

B_values <- c(50, 200, 500)
slope_B <- c()
for (B in B_values){
  val <- rms::validate(full_model$fits[[1]], B=B, bw=TRUE,
                       rule="p", sls=0.05, type="individual", tol=1e-13)
  slope_B <- c(slope_B, val["Slope", "test"])
}

## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  sex     0.11   1    0.7352 0.11     1    0.7352 -1.89
##  hepato  0.57   1    0.4521 0.68     2    0.7119 -3.32
##  ascites 2.16   1    0.1416 2.84     3    0.4170 -3.16
## 
## Approximate Estimates after Deleting Factors
## 
##             Coef    S.E.  Wald Z         P
## age      0.06839 0.02381  2.8723 4.074e-03
## age'    -0.03714 0.02236 -1.6610 9.672e-02
## bili     0.80138 0.10722  7.4745 7.749e-14
## bili'   -2.33364 0.37180 -6.2767 3.459e-10
## albumin -0.83717 0.19953 -4.1956 2.721e-05
## stage=2  0.37999 0.74155  0.5124 6.084e-01
## stage=3  0.79708 0.72348  1.1017 2.706e-01
## stage=4  1.25171 0.72724  1.7212 8.522e-02
## 
## Factors in Final Model
## 
## [1] age     bili    albumin stage  
## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  sex     0.11   1    0.7352 0.11     1    0.7352 -1.89
##  hepato  0.57   1    0.4521 0.68     2    0.7119 -3.32
##  ascites 2.16   1    0.1416 2.84     3    0.4170 -3.16
## 
## Approximate Estimates after Deleting Factors
## 
##             Coef    S.E.  Wald Z         P
## age      0.06839 0.02381  2.8723 4.074e-03
## age'    -0.03714 0.02236 -1.6610 9.672e-02
## bili     0.80138 0.10722  7.4745 7.749e-14
## bili'   -2.33364 0.37180 -6.2767 3.459e-10
## albumin -0.83717 0.19953 -4.1956 2.721e-05
## stage=2  0.37999 0.74155  0.5124 6.084e-01
## stage=3  0.79708 0.72348  1.1017 2.706e-01
## stage=4  1.25171 0.72724  1.7212 8.522e-02
## 
## Factors in Final Model
## 
## [1] age     bili    albumin stage  
## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  sex     0.11   1    0.7352 0.11     1    0.7352 -1.89
##  hepato  0.57   1    0.4521 0.68     2    0.7119 -3.32
##  ascites 2.16   1    0.1416 2.84     3    0.4170 -3.16
## 
## Approximate Estimates after Deleting Factors
## 
##             Coef    S.E.  Wald Z         P
## age      0.06839 0.02381  2.8723 4.074e-03
## age'    -0.03714 0.02236 -1.6610 9.672e-02
## bili     0.80138 0.10722  7.4745 7.749e-14
## bili'   -2.33364 0.37180 -6.2767 3.459e-10
## albumin -0.83717 0.19953 -4.1956 2.721e-05
## stage=2  0.37999 0.74155  0.5124 6.084e-01
## stage=3  0.79708 0.72348  1.1017 2.706e-01
## stage=4  1.25171 0.72724  1.7212 8.522e-02
## 
## Factors in Final Model
## 
## [1] age     bili    albumin stage

names(slope_B) <- paste0("B=", B_values)
round(slope_B, 3)

##  B=50 B=200 B=500 
## 0.898 0.893 0.908

# restore the random number generator
assign(".Random.seed", old_seed, envir=globalenv())

With few bootstrap samples the shrinkage factor is noisier; with more samples the estimates become more stable and are closer to each other. Run it again with another seed to see how much the estimate changes by chance alone.

Exercise 3. Validate the backward selected model bw_model on its own (bw=FALSE, so without repeating the variable selection), and compare it with the validation of the full model with variable selection in every bootstrap sample (bw=TRUE, as above). Use the first imputed data set. Compare the optimism in Dxy and the slope. Which one gives the larger optimism, and why?

Show answer, code and results

# save the state of the random number generator
old_seed <- .Random.seed

# final model only: the variable selection is not repeated
val_final <- rms::validate(bw_model$fits[[1]], B=200, bw=FALSE)

# full model with backward selection repeated in every bootstrap sample
val_strategy <- rms::validate(full_model$fits[[1]], B=200, bw=TRUE,
                              rule="p", sls=0.05, type="individual", tol=1e-13)

## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  sex     0.11   1    0.7352 0.11     1    0.7352 -1.89
##  hepato  0.57   1    0.4521 0.68     2    0.7119 -3.32
##  ascites 2.16   1    0.1416 2.84     3    0.4170 -3.16
## 
## Approximate Estimates after Deleting Factors
## 
##             Coef    S.E.  Wald Z         P
## age      0.06839 0.02381  2.8723 4.074e-03
## age'    -0.03714 0.02236 -1.6610 9.672e-02
## bili     0.80138 0.10722  7.4745 7.749e-14
## bili'   -2.33364 0.37180 -6.2767 3.459e-10
## albumin -0.83717 0.19953 -4.1956 2.721e-05
## stage=2  0.37999 0.74155  0.5124 6.084e-01
## stage=3  0.79708 0.72348  1.1017 2.706e-01
## stage=4  1.25171 0.72724  1.7212 8.522e-02
## 
## Factors in Final Model
## 
## [1] age     bili    albumin stage

# optimism in Dxy and the (test) slope
round(rbind(final_model_only=c(optimism_Dxy=val_final["Dxy", "optimism"],
                               slope=val_final["Slope", "test"]),
            whole_strategy=c(optimism_Dxy=val_strategy["Dxy", "optimism"],
                             slope=val_strategy["Slope", "test"])), 3)

##                  optimism_Dxy slope
## final_model_only        0.017 0.927
## whole_strategy          0.023 0.896

# restore the random number generator
assign(".Random.seed", old_seed, envir=globalenv())

Validating only the final model gives a smaller optimism and a slope closer to 1. This is because variable selection uses the outcome: the variables that remain are partly the ones that look good in this particular data set, also by chance. If you only validate the final, fixed model, you ignore this extra source of overfitting and underestimate the optimism. By repeating the selection in every bootstrap sample, the validation checks the whole modelling strategy, not only the final model.

2 Model performance

A prediction model is more than a C-index. In this part you calculate predicted 10-year risks for three models (the full model, the backward selected model, and the backward selected model with shrinkage), check how well the predicted risks agree with the observed risks (calibration), and assess the clinical usefulness of the model with a decision curve analysis.

2.1 Calculate predictions and make calibration plots

For each of the three models and for each imputed data set we do the following:

1. Calculate the linear predictor (lp) for each patient. For the shrunk model this is the linear predictor of the backward selected model multiplied by the shrinkage factor.
2. Estimate the baseline hazard at the 10-year horizon. For the shrunk model the baseline hazard has to be re-estimated, because the linear predictor changed. This is done with a Cox model that only has the linear predictor as an offset.
3. Turn these into predicted risks at 10 years with fun.event().
4. Make a calibration plot with val.surv.mi(), which pools the performance over the imputed data sets.

Important: this is apparent performance. All performance measures in this part are calculated on the same data that were used to develop the model, so they ignore optimism (overfitting) and are too optimistic. The optimism-corrected C-index in Part 1 corrects for this for the C-index only. For calibration and decision curves you would need to repeat the evaluation with a bootstrap or in new patients (see the internal-external cross-validation in Part 3).

# The three models we compare, and the title of each calibration plot.
# The shrunk model is the backward selected model with shrunk coefficients.
model_names   <- c("full", "bw", "shrunk")
model_titles  <- c(full   = "Full model",
                   bw     = "Final model without shrinkage",
                   shrunk = "Shrunk model")
fitted_models <- list(full=full_model, bw=bw_model, shrunk=bw_model)
model_formulas <- list(full=form_full, bw=form_bw)

# Open one png file with three panels next to each other: one calibration plot
# per model
png(filename="model_performance.png", width=3600, height=1400, res=300)
layout(matrix(1:3, ncol=3, byrow=TRUE))

for (model_name in model_names){
  cat("Now making plot for", model_name, "model \n")

  # predicted 10-year risks of all patients: one column per imputed data set
  predicted_risks <- c()

  for (i in 1:m){
    # the i-th imputed data set
    imputed_i <- mice::complete(imputed_data, i)

    # Step 1: linear predictor of every patient -------------------------------
    lp_i <- predict(fitted_models[[model_name]]$fits[[i]],
                    newdata=imputed_i, type="lp")
    if (model_name == "shrunk"){
      lp_bw_i <- lp_i                     # linear predictor without shrinkage
      lp_i <- shrinkage_factor*lp_bw_i    # shrink it towards zero
    }

    # Step 2: baseline hazard at the horizon ----------------------------------
    if (model_name == "shrunk"){
      # the linear predictor changed, so re-estimate the baseline hazard with a
      # Cox model that only has the linear predictor as an offset
      baseline_fit <- survival::coxph(S_10 ~ offset(lp_i), data=imputed_i)
    } else {
      # refit the Cox model on this imputed data set
      baseline_fit <- rms::cph(model_formulas[[model_name]], data=imputed_i,
                               x=TRUE, y=TRUE, se.fit=TRUE)
    }
    baseline_hazard <- survival::basehaz(baseline_fit)
    # cumulative baseline hazard at the last event time before the horizon
    last_time <- max(baseline_hazard$time[baseline_hazard$time <= horizon])
    h0_i <- baseline_hazard$hazard[baseline_hazard$time == last_time]

    # Step 3: predicted risk at the horizon -----------------------------------
    predicted_risks <- cbind(predicted_risks,
                             PredictionTools::fun.event(lp=lp_i, h0=h0_i))
  }

  # A spline is linear in the tails, so a patient with a very extreme value (for
  # example a very high bilirubin) can get a predicted risk of exactly 1, which
  # the calibration plot cannot handle: keep predictions strictly between 0 and 1
  predicted_risks <- as.matrix(predicted_risks)
  predicted_risks <- pmin(pmax(predicted_risks, 0.0001), 0.9999)

  # save under a name per model: p_full, p_bw and p_shrunk
  assign(paste0("p_", model_name), predicted_risks)

  # Step 4: calibration plot, in the next panel of the png --------------------
  PredictionTools::val.surv.mi(p=predicted_risks, y=S_10, time=horizon, g=3,
                               main=model_titles[[model_name]],
                               CI.metrics=TRUE,
                               show.metrics=c(rep(TRUE, 6), FALSE))
}

## Now making plot for full model

## Now making plot for bw model

## Now making plot for shrunk model

dev.off()

## png 
##   2

2.1.1 Look at the calibration plots

The plots show predicted risk (x-axis) against observed risk (y-axis), with patients grouped by predicted risk. Points on the diagonal mean perfect calibration. The metrics in the plot summarize this (for example calibration intercept and slope). How to read them:

• The calibration intercept < 0 indicates overestimation of the risk, and an intercept > 0 indicates underestimation.
• The calibration slope < 1 indicates overfitting, i.e., a weaker than expected strength of the association between the predictors and the outcome. A slope > 1 indicates underfitting, i.e., a stronger than expected strength of the association between the predictors and the outcome.

2.2 Do the models predict well on average?

Calibration-in-the-large: the average predicted risk should be close to the observed 10-year risk from the Kaplan-Meier curve (first number).

# check if models predict well on average
# observed 10-year risk: 1 minus the Kaplan-Meier survival at the horizon
1-summary(survival::survfit(S_10~1, data=data), time=horizon)$surv

## [1] 0.5578323

# average predicted risk: rowMeans() averages the predictions of each patient
# over the imputed data sets, mean() averages those over all patients
mean(rowMeans(p_full))

## [1] 0.5710267

mean(rowMeans(p_bw))

## [1] 0.5717457

mean(rowMeans(p_shrunk))

## [1] 0.5829371

The unshrunk models fit the development data best. The shrunk model is the furthest from the observed risk: because we shrank the coefficients, its predictions are less extreme and it underfits the development data a little. That is intended: shrinkage gives up some fit on these patients in exchange for better predictions in new patients, where the unshrunk model would be too extreme because it is overfitted. So do not expect the shrunk model to look best on the development data.

2.3 Decision curve analysis (DCA)

A decision curve shows the net benefit of using the model to make decisions (for example, treat patients with a predicted risk above a threshold), compared with treating everyone or no one. A model is clinically useful at the thresholds where its curve lies above both of these reference strategies.

# decision curve analysis (DCA)
# data.frame() puts the outcome (y) and the predicted risk of each patient (Model)
# in a table; dca() calculates the net benefit at the horizon, and the pipe |>
# passes the result on to plot(), which draws the curves (smooth=TRUE smooths them)
dca_plot <- dcurves::dca(y ~ Model,
                         time=horizon,
                         data=data.frame(y=S_10,
                                         Model=rowMeans(p_bw))) |>
  plot(smooth = TRUE)
dca_plot

The histogram shows how the predicted risks are distributed. Decision curves are most reliable in the range of thresholds where there are enough patients.

# histogram of the predicted risks: ggplot figures are built in layers that are
# added with +
hist_pred <- ggplot2::ggplot(data=data.frame(Predictions=rowMeans(p_bw)),
                             ggplot2::aes(x=Predictions)) +
  ggplot2::geom_histogram(fill="#619bff", color="grey", alpha=0.6, binwidth=0.01) +
  ggplot2::ylab("Count") +
  ggplot2::theme_minimal() +
  ggplot2::theme(panel.background = ggplot2::element_rect(fill="white", color=NA),
        plot.background = ggplot2::element_rect(fill="white", color=NA))
hist_pred

# save: ggarrange() puts the two figures below each other, ggsave() writes
# the result to a png file in the working directory
ggplot2::ggsave(file="dca.png",
       plot=ggpubr::ggarrange(dca_plot, hist_pred,
                      nrow=2, ncol=1, heights=c(2, 1), align="v"),
       width=16, height=16, units="cm", dpi=300)

2.4 Exercises: Part 2

Try each exercise yourself first. Then open the field below it to see the answer, the code and its results.

Exercise 1. Compare the three calibration plots. Which model has the best calibration slope? Does this match what you expect from shrinkage?

Show answer, code and results

We calculate the calibration slope ourselves for the last imputed data set (the objects lp_bw_i and lp_i are left over from the last round of the loop above): the slope of a Cox model with the linear predictor as the only predictor.

cal_slope_bw <- coef(survival::coxph(S_10 ~ lp_bw_i))
cal_slope_shrunk <- coef(survival::coxph(S_10 ~ lp_i))
round(c(without_shrinkage=cal_slope_bw,
        shrunk=cal_slope_shrunk,
        one_over_shrinkage_factor=1/shrinkage_factor), 3)

## without_shrinkage.lp_bw_i               shrunk.lp_i one_over_shrinkage_factor 
##                     1.000                     1.101                     1.101

On the development data, a fitted Cox model has a calibration slope of 1 by construction, so the model without shrinkage looks perfect, while the shrunk model has a slope of 1 divided by the shrinkage factor, which is above 1. You will most likely see the same pattern in the plots (the plotted slope may be calculated slightly differently, but the idea is the same). So the answer does not match what you might expect from shrinkage, and that is the point: shrinkage makes the model look worse in the data it was developed on, but it improves calibration in new patients, because the shrinkage factor was estimated with the bootstrap to correct for exactly this optimism.

Exercise 2. Is the average predicted risk close to the observed 10-year risk for all three models?

Show answer, code and results

observed_risk <- 1-summary(survival::survfit(S_10~1, data=data), time=horizon)$surv
round(c(observed=observed_risk,
        full=mean(rowMeans(p_full)),
        bw=mean(rowMeans(p_bw)),
        shrunk=mean(rowMeans(p_shrunk))), 3)

## observed     full       bw   shrunk 
##    0.558    0.571    0.572    0.583

If the average predicted risks are close to the observed risk, the models are well calibrated “in the large”. Note that for the shrunk model we re-estimated the baseline hazard after shrinking the linear predictor, but it still fits the development data a bit less well than the unshrunk models: due to the shrinkage it underfits these patients, in order to predict new patients better.

Exercise 3. At which threshold probabilities does the model have a higher net benefit than treating everyone or no one?

Show answer, code and results

Read this from the decision curve plot: the model is useful at the thresholds where its line is above both the “treat all” and the “treat none” lines. Only look at thresholds where there are patients with such predicted risks:

range(rowMeans(p_bw))

## [1] 0.05037855 0.99990000

Outside this range of predicted risks the curve is not informative, because (almost) no patients have such a predicted risk. Also note that the decision curve is based on the development data, so it is optimistic.

Exercise 4. The DCA above uses p_bw. Replace it with p_shrunk (the final model with shrinkage). Does the curve change?

Show answer, code and results

dca_shrunk <- dcurves::dca(y ~ Model,
                           time=horizon,
                           data=data.frame(y=S_10,
                                           Model=rowMeans(p_shrunk))) |>
  plot(smooth = TRUE)
dca_shrunk

Shrinkage does not change the ranking of the patients (so the C-index is the same), but it makes the predicted risks less extreme. A decision curve depends on the predicted risks, because the threshold is applied to them, so the curve can still change somewhat.

Exercise 5. Why are these performance measures optimistic, and how would you get an honest estimate?

Show answer

There is no code for this one. The models were developed and evaluated on the same patients, so the performance measures partly reflect noise in this data set that the model has fitted. For an honest estimate, evaluate the model on data that were not used to develop it: new patients (external validation), or repeat the whole modelling procedure in bootstrap samples and correct for the optimism (as we did for the C-index in Part 1), or use internal-external cross-validation (Part 3).

3 Internal-external cross-validation

Optimism-corrected performance (Part 1) uses random bootstrap samples. If your data come from several centers or regions, you can do something stronger: internal-external cross-validation. You leave out one center, develop the whole model (including variable selection and shrinkage) on the other centers, and test it on the center that was left out. You repeat this for every center. This shows how well the modelling strategy generalizes to new settings.

Note: the center variable in this tutorial is artificial (it alternates between 1 and 2), so there is no real heterogeneity between centers. This is fine for a demonstration of the code, but the results will not tell you anything about transportability.

3.1 Prepare the term names of the full model

Backward selection with fastbw() returns the names of the variables that were kept, for example age. But in the full model age and bilirubin are entered as rcs(age, 3) and rcs(bili, 3) (splines, so non-linear effects). If we would build the new formula from the variable names only, they would become linear. We therefore look up the full term specification for every variable in form_full, so the splines are kept.

# term labels of the full model, e.g. "rcs(age, 3)", "sex", ...
# (terms() splits the formula in its terms; the labels are stored as an attribute)
term_labels <- attr(terms(form_full), "term.labels")

# variable name that belongs to each term, e.g. "age", "sex", ...
# sapply() applies a function to every term label; all.vars() finds the variable
# names in a piece of R code and [1] takes the first one
term_vars <- sapply(term_labels, function(x) all.vars(parse(text=x))[1])
data.frame(term_labels, term_vars, row.names=NULL)

3.2 Run the internal-external cross-validation

For each center in turn (the left-out center) we do the following:

1. Fit the full model on the data excluding the left-out center, and perform backward selection with fastbw(). We use the same settings as in the main analysis (type="individual", p < 0.05).
2. Refit the selected model on the development centers.
3. Calculate the shrinkage factor on the development centers with the bootstrap.
4. Calculate the baseline hazard of the shrunk model on the development centers.
5. Predict the risk for the patients in the left-out center and make a calibration plot.

Note: the data were imputed once in Tutorial 1, using all centers. Strictly, you should also impute within the development regions only, so that no information from the left-out region leaks into the model. For simplicity we do not do this here.

# Open a png file with one panel per center (one calibration plot per center)
png(file="Internal_external_cross_validation.png",
    width=4000, height=2000, res=300)
layout(matrix(1:2, ncol=2, byrow=TRUE))

# Containers to keep track of what happens in every round of the loop:
# the terms of the final model and its shrinkage factor, starting with the
# model that was developed on all centers (Part 1)
selected_terms_per_model <- paste(names(bw_model$coefficients), collapse=", ")
shrinkage_per_model <- c(shrinkage_factor)

# Every center is left out once ("left_out_center" is the validation center)
for (left_out_center in unique(data$center)){

  # ---- Which patients are used for development and which for validation? ----
  is_development <- data$center != left_out_center   # TRUE for the other centers
  is_validation  <- data$center == left_out_center   # TRUE for the left-out center
  S_10_development <- S_10[is_development, ]         # outcome of development patients
  S_10_validation  <- S_10[is_validation, ]          # outcome of validation patients

  # ---- Step 1: full model on the development centers, then backward selection ----
  # the outcome in the formula must only contain the development patients, so we
  # replace S_10 by the outcome of the development patients
  full_formula_dev <- update(form_full, S_10_development ~ .)
  full_model_dev <- fit.mult.impute(full_formula_dev, cph,
                                    n.impute=m,
                                    fitargs=list(x=TRUE, y=TRUE, surv=TRUE),
                                    xtrans=imputed_data,
                                    data=data,
                                    sub=is_development,   # only the development patients
                                    fit.reps=TRUE)

  # same settings as in the main analysis: type="individual", p < 0.05
  selected_variables <- fastbw(full_model_dev,
                               type="individual",
                               rule="p", sls=0.05)$names.kept

  # ---- Step 2: refit the selected model on the development centers ----
  # use the full term specification so that rcs(age, 3) stays a spline
  selected_terms <- term_labels[term_vars %in% selected_variables]
  selected_terms_per_model <- c(selected_terms_per_model,
                                paste(selected_terms, collapse=", "))
  final_formula_dev <- as.formula(paste("S_10_development ~",
                                        paste(selected_terms, collapse=" + ")))
  final_model_dev <- fit.mult.impute(final_formula_dev, cph,
                                     n.impute=m,
                                     fitargs=list(x=TRUE, y=TRUE, surv=TRUE),
                                     xtrans=imputed_data,
                                     data=data,
                                     sub=is_development,
                                     fit.reps=TRUE)

  # ---- Steps 3-5: for every imputed data set, shrinkage factor and predictions ----
  shrinkage_per_imputation <- c()         # shrinkage factor of each imputed data set
  predicted_risks_validation <- c()       # predicted risks in the left-out center
  for (i in 1:m){
    # the i-th imputed data set, split in development and validation patients
    imputed_i <- mice::complete(imputed_data, i)
    development_i <- imputed_i[is_development, ]
    validation_i  <- imputed_i[is_validation, ]

    # Step 3: shrinkage factor from the bootstrap on the development centers
    # (backward selection is repeated in every bootstrap sample)
    validation_results <- validate(full_model_dev$fits[[i]],
                                   pr=FALSE,
                                   B=200, bw=TRUE, rule="p", sls=0.05,
                                   type="individual", tol=1e-13)
    shrinkage_i <- validation_results["Slope", "test"]
    shrinkage_per_imputation <- c(shrinkage_per_imputation, shrinkage_i)

    # Step 4: baseline hazard of the shrunk model on the development centers
    # linear predictor of the development patients, shrunk towards zero
    lp_development <- predict(final_model_dev$fits[[i]],
                              newdata=development_i, type="lp")
    lp_development_shrunk <- shrinkage_i*lp_development
    # Cox model with only the shrunk linear predictor as offset gives the baseline hazard
    baseline_fit <- survival::coxph(S_10_development ~ offset(lp_development_shrunk))
    baseline_hazard <- survival::basehaz(baseline_fit)
    # cumulative baseline hazard at the last event time before the horizon
    last_time <- max(baseline_hazard$time[baseline_hazard$time <= horizon])
    h0_i <- baseline_hazard$hazard[baseline_hazard$time == last_time]

    # Step 5: predicted risks for the patients in the left-out center
    lp_validation_shrunk <- shrinkage_i*predict(final_model_dev$fits[[i]],
                                                newdata=validation_i, type="lp")
    predicted_risks_validation <- cbind(predicted_risks_validation,
                                        PredictionTools::fun.event(lp=lp_validation_shrunk,
                                                                   h0=h0_i))
  }

  # keep the mean shrinkage factor over the imputed data sets for this round
  shrinkage_per_model <- c(shrinkage_per_model, mean(shrinkage_per_imputation))

  # keep predictions strictly between 0 and 1 (see the note above)
  predicted_risks_validation <- pmin(pmax(predicted_risks_validation, 0.0001), 0.9999)

  # ---- Calibration plot for the left-out center ----
  PredictionTools::val.surv.mi(p=predicted_risks_validation, y=S_10_validation,
                               time=horizon, g=2,
                               CI.metrics=TRUE, main=left_out_center,
                               show.metrics=c(rep(TRUE, 6), FALSE))
}

## 
## Wald Statistic Information
## 
## Variance Inflation Factors Due to Imputation:
## 
##         age        age'        bili       bili'     albumin  sex=Female 
##        1.01        1.00        1.01        1.01        1.01        1.01 
## ascites=Yes  hepato=Yes     stage=2     stage=3     stage=4 
##        1.08        1.11        1.11        1.14        1.16 
## 
## Fraction of Missing Information:
## 
##         age        age'        bili       bili'     albumin  sex=Female 
##        0.01        0.00        0.01        0.01        0.01        0.01 
## ascites=Yes  hepato=Yes     stage=2     stage=3     stage=4 
##        0.07        0.10        0.10        0.12        0.13 
## 
## d.f. for t-distribution for Tests of Single Coefficients:
## 
##         age        age'        bili       bili'     albumin  sex=Female 
##    66572.61   332639.56    54254.44    50821.03    62116.70    63436.06 
## ascites=Yes  hepato=Yes     stage=2     stage=3     stage=4 
##      714.21      379.60      393.47      278.47      219.94 
## 
## The following fit components were averaged over the 5 model fits:
## 
##   linear.predictors means stats center 
## 
## 
## Wald Statistic Information
## 
## Variance Inflation Factors Due to Imputation:
## 
##     age    age'    bili   bili' albumin 
##       1       1       1       1       1 
## 
## Fraction of Missing Information:
## 
##     age    age'    bili   bili' albumin 
##       0       0       0       0       0 
## 
## d.f. for t-distribution for Tests of Single Coefficients:
## 
##     age    age'    bili   bili' albumin 
##     Inf     Inf     Inf     Inf     Inf 
## 
## The following fit components were averaged over the 5 model fits:
## 
##   linear.predictors means stats center 
## 
## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  sex     0.09   1    0.7673 0.09     1    0.7673 -1.91
##  stage   1.64   3    0.6512 1.72     4    0.7864 -6.28
##  ascites 0.96   1    0.3265 2.69     5    0.7482 -7.31
## 
## Approximate Estimates after Deleting Factors
## 
##                Coef    S.E.  Wald Z         P
## age         0.03819 0.03166  1.2064 2.277e-01
## age'        0.01178 0.03101  0.3798 7.041e-01
## bili        0.77535 0.15184  5.1065 3.282e-07
## bili'      -1.85404 0.42622 -4.3500 1.361e-05
## albumin    -1.48511 0.28401 -5.2290 1.704e-07
## hepato=Yes  0.56938 0.26562  2.1436 3.206e-02
## 
## Factors in Final Model
## 
## [1] age     bili    albumin hepato

## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  sex     0.06   1    0.8005 0.06     1    0.8005 -1.94
##  stage   2.91   3    0.4050 2.98     4    0.5615 -5.02
##  ascites 1.88   1    0.1699 4.86     5    0.4330 -5.14
##  hepato  3.77   1    0.0520 8.64     6    0.1951 -3.36
## 
## Approximate Estimates after Deleting Factors
## 
##             Coef    S.E.  Wald Z         P
## age      0.03802 0.03159  1.2038 2.287e-01
## age'     0.01198 0.03113  0.3849 7.003e-01
## bili     0.83894 0.15028  5.5827 2.368e-08
## bili'   -2.02884 0.42276 -4.7990 1.594e-06
## albumin -1.57702 0.28483 -5.5366 3.083e-08
## 
## Factors in Final Model
## 
## [1] age     bili    albumin
## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  sex     0.13   1    0.7188 0.13     1    0.7188 -1.87
##  ascites 0.26   1    0.6085 0.39     2    0.8220 -3.61
##  hepato  0.63   1    0.4276 1.02     3    0.7961 -4.98
##  stage   5.16   3    0.1605 6.18     6    0.4033 -5.82
## 
## Approximate Estimates after Deleting Factors
## 
##             Coef    S.E.  Wald Z         P
## age      0.03737 0.03170  1.1788 2.385e-01
## age'     0.01256 0.03118  0.4028 6.871e-01
## bili     0.83756 0.15049  5.5656 2.612e-08
## bili'   -2.03034 0.42394 -4.7892 1.674e-06
## albumin -1.57049 0.28673 -5.4771 4.322e-08
## 
## Factors in Final Model
## 
## [1] age     bili    albumin
## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  stage   1.62   3    0.6556 1.62     3    0.6556 -4.38
##  sex     0.24   1    0.6255 1.85     4    0.7624 -6.15
##  ascites 0.70   1    0.4030 2.55     5    0.7683 -7.45
## 
## Approximate Estimates after Deleting Factors
## 
##                Coef    S.E.  Wald Z         P
## age         0.04137 0.03180  1.3007 1.934e-01
## age'        0.01011 0.03117  0.3244 7.456e-01
## bili        0.78976 0.15021  5.2578 1.458e-07
## bili'      -1.89340 0.42132 -4.4940 6.991e-06
## albumin    -1.44804 0.29169 -4.9644 6.893e-07
## hepato=Yes  0.55255 0.26009  2.1244 3.364e-02
## 
## Factors in Final Model
## 
## [1] age     bili    albumin hepato

## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  sex     0.13   1    0.7212 0.13     1    0.7212 -1.87
##  ascites 0.22   1    0.6364 0.35     2    0.8391 -3.65
##  hepato  0.84   1    0.3580 1.20     3    0.7540 -4.80
##  stage   5.32   3    0.1498 6.52     6    0.3679 -5.48
## 
## Approximate Estimates after Deleting Factors
## 
##             Coef    S.E.  Wald Z         P
## age      0.03770 0.03169  1.1898 2.341e-01
## age'     0.01194 0.03115  0.3833 7.015e-01
## bili     0.83437 0.15027  5.5524 2.818e-08
## bili'   -2.02149 0.42337 -4.7747 1.800e-06
## albumin -1.56362 0.28645 -5.4586 4.798e-08
## 
## Factors in Final Model
## 
## [1] age     bili    albumin

## 
## Wald Statistic Information
## 
## Variance Inflation Factors Due to Imputation:
## 
##         age        age'        bili       bili'     albumin  sex=Female 
##        1.02        1.04        1.03        1.03        1.04        1.03 
## ascites=Yes  hepato=Yes     stage=2     stage=3     stage=4 
##        1.23        1.06        1.01        1.00        1.01 
## 
## Fraction of Missing Information:
## 
##         age        age'        bili       bili'     albumin  sex=Female 
##        0.02        0.04        0.03        0.03        0.04        0.03 
## ascites=Yes  hepato=Yes     stage=2     stage=3     stage=4 
##        0.19        0.06        0.01        0.00        0.01 
## 
## d.f. for t-distribution for Tests of Single Coefficients:
## 
##         age        age'        bili       bili'     albumin  sex=Female 
##     7782.08     2437.34     4020.35     4711.53     2839.42     4813.61 
## ascites=Yes  hepato=Yes     stage=2     stage=3     stage=4 
##      116.07     1090.99   124286.36   863231.11    98640.49 
## 
## The following fit components were averaged over the 5 model fits:
## 
##   linear.predictors means stats center 
## 
## 
## Wald Statistic Information
## 
## Variance Inflation Factors Due to Imputation:
## 
##         age        age'        bili       bili' ascites=Yes 
##        1.02        1.04        1.03        1.02        1.20 
## 
## Fraction of Missing Information:
## 
##         age        age'        bili       bili' ascites=Yes 
##        0.02        0.04        0.03        0.02        0.17 
## 
## d.f. for t-distribution for Tests of Single Coefficients:
## 
##         age        age'        bili       bili' ascites=Yes 
##    11956.63     2561.00     4814.62     8212.24      142.89 
## 
## The following fit components were averaged over the 5 model fits:
## 
##   linear.predictors means stats center 
## 
## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  hepato  0.01   1    0.9357 0.01     1    0.9357 -1.99
##  albumin 0.52   1    0.4729 0.52     2    0.7704 -3.48
##  sex     2.15   1    0.1430 2.67     3    0.4458 -3.33
##  stage   7.14   3    0.0674 9.81     6    0.1328 -2.19
## 
## Approximate Estimates after Deleting Factors
## 
##                Coef    S.E. Wald Z         P
## age          0.1262 0.03983  3.167 1.540e-03
## age'        -0.1083 0.04203 -2.577 9.960e-03
## bili         1.0661 0.17377  6.135 8.519e-10
## bili'       -3.7125 0.68721 -5.402 6.579e-08
## ascites=Yes  0.8852 0.36108  2.452 1.422e-02
## 
## Factors in Final Model
## 
## [1] age     bili    ascites

## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  albumin 0.33   1    0.5680  0.33    1    0.5680 -1.67
##  hepato  0.25   1    0.6197  0.57    2    0.7512 -3.43
##  sex     2.81   1    0.0940  3.38    3    0.3370 -2.62
##  stage   7.78   3    0.0508 11.16    6    0.0837 -0.84
## 
## Approximate Estimates after Deleting Factors
## 
##                Coef    S.E. Wald Z         P
## age          0.1235 0.04043  3.054 2.259e-03
## age'        -0.1114 0.04238 -2.627 8.610e-03
## bili         1.0664 0.17562  6.072 1.262e-09
## bili'       -3.7327 0.69420 -5.377 7.572e-08
## ascites=Yes  1.0448 0.34171  3.058 2.231e-03
## 
## Factors in Final Model
## 
## [1] age     bili    ascites

## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  hepato  0.41   1    0.5219 0.41     1    0.5219 -1.59
##  albumin 0.94   1    0.3310 1.36     2    0.5079 -2.64
##  sex     1.83   1    0.1762 3.18     3    0.3641 -2.82
##  stage   5.34   3    0.1486 8.52     6    0.2022 -3.48
## 
## Approximate Estimates after Deleting Factors
## 
##                Coef    S.E. Wald Z         P
## age          0.1331 0.04186  3.180 1.474e-03
## age'        -0.1186 0.04323 -2.744 6.062e-03
## bili         1.0362 0.17433  5.944 2.784e-09
## bili'       -3.6308 0.68553 -5.296 1.182e-07
## ascites=Yes  1.2113 0.35425  3.419 6.278e-04
## 
## Factors in Final Model
## 
## [1] age     bili    ascites

## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  hepato  0.19   1    0.6598 0.19     1    0.6598 -1.81
##  albumin 0.62   1    0.4309 0.81     2    0.6656 -3.19
##  sex     2.68   1    0.1015 3.50     3    0.3213 -2.50
## 
## Approximate Estimates after Deleting Factors
## 
##                Coef    S.E.  Wald Z         P
## age          0.1083 0.04082  2.6525 7.990e-03
## age'        -0.1047 0.04279 -2.4467 1.442e-02
## bili         1.0086 0.18158  5.5548 2.780e-08
## bili'       -3.5320 0.71667 -4.9283 8.294e-07
## ascites=Yes  0.9394 0.35670  2.6336 8.448e-03
## stage=2     -0.3880 1.07814 -0.3599 7.189e-01
## stage=3      0.5584 1.04055  0.5366 5.915e-01
## stage=4      0.9056 1.04241  0.8688 3.850e-01
## 
## Factors in Final Model
## 
## [1] age     bili    ascites stage

## 
##      Backwards Step-down - Original Model
## 
##  Deleted Chi-Sq d.f. P      Residual d.f. P      AIC  
##  hepato  0.25   1    0.6156 0.25     1    0.6156 -1.75
##  albumin 0.87   1    0.3507 1.12     2    0.5703 -2.88
##  sex     1.40   1    0.2366 2.52     3    0.4710 -3.48
## 
## Approximate Estimates after Deleting Factors
## 
##                 Coef    S.E.  Wald Z         P
## age          0.10432 0.04031  2.5878 9.658e-03
## age'        -0.09357 0.04140 -2.2600 2.382e-02
## bili         0.94415 0.18016  5.2405 1.601e-07
## bili'       -3.28659 0.71022 -4.6275 3.700e-06
## ascites=Yes  0.86770 0.34026  2.5501 1.077e-02
## stage=2     -0.23685 1.07737 -0.2198 8.260e-01
## stage=3      0.48136 1.04051  0.4626 6.436e-01
## stage=4      0.97350 1.04011  0.9360 3.493e-01
## 
## Factors in Final Model
## 
## [1] age     bili    ascites stage

dev.off()

## png 
##   2

3.3 Exercises: Part 3

Try each exercise yourself first. Then open the field below it to see the answer, the code and its results.

Exercise 1. Which variables are selected when center 1 or center 2 is left out? Are they the same as in the final model? (Tip: see the objects selected_terms_per_model and shrinkage_per_model in the code above.)

Show answer, code and results

The objects selected_terms_per_model and shrinkage_per_model hold exactly this information. Here we put them in one table:

data.frame(development_data=c("All data (final model)",
                              paste("Without", unique(as.character(data$center)))),
           selected_terms=selected_terms_per_model,
           shrinkage_factor=round(shrinkage_per_model, 3))

Compare the selected terms in each row with the final model on all data. Different selections in different development sets show that variable selection is unstable.

Exercise 2. Compare the calibration plots of the two centers. Would you expect big differences in real data with truly different centers? What could cause them?

Show answer

There is no code for this one. In this tutorial the centers are artificial, so the two plots should look alike, and differences are due to chance. In real data with truly different centers you would expect larger differences. Possible causes are differences in case mix (the patients differ), in how predictors and outcomes are measured, in treatment practice, and in follow-up.

Exercise 3. Why do we develop the model on the other centers and not on a random half of the data?

Show answer

There is no code for this one. A random split gives validation data that are very similar to the development data, so it mainly shows internal validity and wastes half of the data. Leaving out a whole center tests the model in a different setting from the development settings, which is closer to how the model will be used in practice, while still using most of the data for development.

Exercise 4. How would you change the code if you had 5 centers instead of 2?

Show answer and code (not run)

The loop already works for any number of centers, because it uses unique(data$center). Only the figure needs more room. This code is not run, because the tutorial data have 2 centers:

# one panel per center in the figure
n_centers <- length(unique(data$center))
png(file="Internal_external_cross_validation.png",
    width=2000*n_centers, height=2000, res=300)
layout(matrix(1:n_centers, ncol=n_centers, byrow=TRUE))

# ... the same loop as above ...

Also think about the number of patients per center. With 5 centers each center has fewer patients, so you may want fewer risk groups in the calibration plot (the argument g=2 of val.surv.mi). Each development set now contains 4 of the 5 centers, so more of the data is used for development than with 2 centers.