In this case study we show how to implement an estimator for the leave-one-out (LOO) version of the expected Brier score (denoted by LOO-BS) for Bayesian survival models, with potentially right-censored event times. We combine the “inverse of probability of censoring weighted estimator” by Gerds et al. with Pareto-smoothed importance sampling (PSIS), by Vehtari et al., in order to efficiently approximate the LOO-BS. We compare the PSIS approximation to the exact LOO-BS.

As a by-product we show how to implement proportional hazard models with a time-dependent baseline hazard by using monotone splines (M-splines).

Useful References:

For PSIS we recommend:

Vehtari, Aki, Andrew Gelman, and Jonah Gabry. “Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC.” Statistics and Computing 27.5 (2017): 1413-1432. APA

For the inverse of probability of censoring weighted estimator see:

Gerds, Thomas A., and Martin Schumacher. “Consistent estimation of the expected Brier score in general survival models with right‐censored event times.” Biometrical Journal 48.6 (2006): 1029-1040.

Monotone (M-splines) were introduced here:

Ramsay, James O. “Monotone regression splines in action.” Statistical science 3.4 (1988): 425-441.

See als this Stan discourse discussion.

knitr::opts_chunk$set(echo = TRUE)
library(rstan);rstan_options(auto_write = TRUE)
library(loo)
library(survival)
library(tibble)
library(dplyr)
library(purrr)
library(bayesplot)
library(splines2)
library(pec)
library(cowplot)
library(MCMCpack)
library(RColorBrewer)
qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))

sm <- stan_model("~/Desktop/Stan/Loo_Pec/survival_parametric_baseline_hazard_simplex_loo.stan")
expose_stan_functions("~/Desktop/Stan/Loo_Pec/survival_parametric_baseline_hazard_simplex_loo.stan")
df <- read.delim("~/Desktop/Stan/Loo_Pec/ncog.txt",sep=" ")

df <- df  %>% mutate(event=d,time=t*12/365) 
df <- df%>% mutate(isGroupB=as.double(arm=="B")) %>%
  dplyr::select(isGroupB, event, time)
N <- nrow(df)
times <- pull(df,time)
is_censored <- pull(df,event)==0
msk_censored <- is_censored == 1
time_range <- range(times)
time_min <- time_range[1]
time_max <- time_range[2]
X <- as.matrix(pull(df, isGroupB))

The dataset

The table below reports survival data from a randomized control trial run by the Nothern California Oncology Group, that compares two treatments for head and neck cancer: isGroupB==0 in case of chemotherapy and isGroup==1 if chemotherapy and radiation was applied. The response variable time is the survival time in months (actually days divided by 12/365) and the variable event is equal to 0 iff the survival time of that patient is only known to exceed the corresponding survival time. For details see (starting) chapter 9.2 in “Computer Age Statistical Inference” by B. Efron and T. Hastie.

Prior

We build a proportional hazard model with a baseline-hazard that is parametrized by M-Splines. Below we show prior predictive (checks) draws for the time-dependent baseline1 hazard and survival functions.

For details we refer to accompanying Stan Code in survival_parametric_baseline_hazard_simplex_loo.stan.

The functions surv_t and hazard_t are defined in the above file and are exposed above as R functions through expose_stan_functions. They allow the evaluation of the survival and hazard function, respectively, at fixed time, for different (latent) paremters and covariates. Both functions return a matrix with as many rows as covariates and as many columns as parameters. For their respective signatures, refer to the definition in the Stan file.