Predictive performance of joint models versus standard survival models

Question

I am trying to show that predictions based on repeated measures of markers (using joint modelling of repeated markers and time to event models: JMbayes package) are better than those based on only one single measure of these markers (using standard survival model). This is in fact the main motivation of using joint modelling approach. I did this for the pbc2 data and was surprised to find out that the AUC from the standard model was similar or even better than the one from the joint model! Any explanation? Please see the code below.

library(JMbayes)
library(timeROC)
library(dplyr)
library(purrr)

#percentage of training and prediction time horizon
train = 0.7
timeHorizon = 3
#Read the pbc2 data#
data("pbc2")
pbc2$event <- as.numeric(pbc2$status != "alive")
pbc2$id    <- as.numeric(as.character(pbc2$id))
#Use complete case analysis: remove missing values for dependent and independent covariates#
pbc2 <- pbc2 %>% filter(!is.na(serBilir), !is.na(spiders), !is.na(age), !is.na(drug) )
pbc2 <- pbc2 %>% group_by(id) %>% 
                 mutate(Flag=if_else((years - year) < timeHorizon & event==0 , NA_real_, abs(years - year))) %>%  
                filter(!is.na(Flag)) %>% ungroup() 

#Generate training and validation data sets#
set.seed(123)
train_ids      <- sample(unique(pbc2$id), size = ceiling(n_distinct(pbc2$id)*train), replace = FALSE)
pbc2_longtrain <- filter(pbc2, id %in% train_ids)
pbc2_ep1train  <- pbc2_longtrain %>% filter(year==0) 
pbc2_longval   <- filter(pbc2, !(id %in% train_ids))

#######################################
#A. Deriving the predictive tools   ###
#######################################

#A1. Standard survival: Based on baseline observations
fit_stand <- coxph(Surv(years, event) ~ drug + age + log(serBilir) + spiders, data = pbc2_ep1train, model = TRUE)
#A2. Joint survival Model: Based on repeated measures of serBilir and spiders
MixedModelFit <- mvglmer(list(log(serBilir) ~ year + (year | id),
                              spiders ~ year + (1 | id)), 
                         data = pbc2_longtrain,
                         families = list(gaussian, binomial))

CoxFit        <- coxph(Surv(years, event) ~ drug + age, data = pbc2_ep1train, model = TRUE)
fit_joint         <- mvJointModelBayes(MixedModelFit, CoxFit, timeVar = "year")

#########################################
#B. Validation of the predictive tools ##
########################################

#B.1. Standard survival method#
    #Validation Data set to be used#
df_val_stand <- pbc2_longval %>% group_by(id) %>% 
                           mutate(age=age+year, years_orig=years, years=timeHorizon, Y.s=years_orig - year) %>% 
                           slice(n()) %>% ungroup()

    #Prediction using current values#
df_val_stand$S_stand <- exp(-predict(fit_stand, newdata = df_val_stand, type="expected")) 
AUC_stand      <- timeROC(T=df_val_stand$Y.s, delta=df_val_stand$event, marker=I(1 - df_val_stand$S_stand), cause=1, times=timeHorizon)

#B2. Joint Model
   #Validation Data set to be used#
df_val_joint <- pbc2_longval %>% group_by(id) %>% 
                                mutate(s = last(year), t_plus_s = s + timeHorizon, Y.s = years-s) %>%
                                ungroup()

   #Prediction using repeated markers#
d.sJointE   <- df_val_joint %>% group_split(id) %>% map_dfr(function(sgr) {
  predict <- JMbayes::survfitJM(object = fit_joint, newdata = sgr, last.time = unique(sgr$s) , survTimes = unique(sgr$t_plus_s), idVar = "id")
  res <- data.frame(id=sgr$id[1], years=sgr$years[1], s=sgr$s[1], Y.s=sgr$Y.s[1], t_plus_s=sgr$t_plus_s[1], event=sgr$event[1], predict[['summaries']][[1]])
  res
})

AUC_joint <- timeROC(T=d.sJointE$Y.s, delta=d.sJointE$event, marker=I(1 - d.sJointE$Mean), cause=1, times=timeHorizon)

#################
#C. Results #####
#################
AUC_stand
# Time-dependent-Roc curve estimated using IPCW  (n=93, without competing risks). 
# Cases Survivors Censored AUC (%)
# t=0     0        93        0      NA
# t=3    45        48        0   90.28
# 
# Method used for estimating IPCW:marginal 
# 
# Total computation time : 0.01  secs.

AUC_joint
# Time-dependent-Roc curve estimated using IPCW  (n=93, without competing risks). 
# Cases Survivors Censored AUC (%)
# t=0     0        93        0      NA
# t=3    45        48        0   89.86
# 
# Method used for estimating IPCW:marginal 
# 
# Total computation time : 0  secs.

Answer 1

Echoing the comment of Dimitris, aucJM() is the JMbayes function for computing time-varying AUC or AUC(t). Because the conventional Cox model cannot be estimated in JMbayes, it is unclear if your goal can be achieved. Using timeROC() with the Cox model and aucJM() with the joint model might be an apples-to-oranges comparison. The manual page for timeROC() is a bit vague (to me), but the function possibly computes the IPCW estimate of the c-index, which is not the same as the AUC(t) of aucJM(). The IPCW c-index forces censoring at time t and uses the observed event times, whereas IPCW AUC(t) uses the binary event status at t. As Blanche et al. (2019) show in their Figure 1, the IPCW c-index can be larger or smaller than IPCW AUC(t) depending on t. Perhaps your question of the importance of repeated measures can be addressed by using the extended joint models or alternative association structures available in JMbayes (see this vignette). For example, you can specify a model to examine the effect of the current value and the slope at time t. You will have to think it through, but it seems that if the repeated measurements are important, then the slope will show an effect. Alternatively, you can specify a random-intercepts and random-slopes model to examine the predictive effect of baseline level and rate of change.

Answer 2

Dear Dimitris and Jeffrey, Thank you all for your comments. I just realized now that the comparison between standard and joint survival models cannot be done in the way I described above. The main reason for this is that for each value s, the joint model provides a prediction rule P(s) = (T > t + s | T > s) which depends on the history of the markers up to s. In my example I was, unintentionally, creating a composite prediction rule P = (P(s_1), P(s_2), ...., P(s_n)) where n is the number of subjects in the validation set and I was comparing the predictive ability of this composite prediction rule P to the one of the standard survival model. Saying that, I still believe that it is very important to show, from a real validation data sets, that predictions made from joint models are better than those made from standard survival models. I came across a paper published in statistics in medicine (2018) mentioning another paper published in the same journal in 2016 where apparently this issue was addressed. They said: "In terms of survival predictions, Yang et al. have recently found increased predictive performance for a variety of joint models compared with baseline‐only and time‐dependent Cox models". I looked at the original paper of Yang et al (2016) but I didn't quite understand how they did it (it wasn't explained). The paper can be accessed here (see page 10, paragraph starting with For JM1 and JM2 models).

All the best

Abderrahim