# Predictive performance of joint models versus standard survival models

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
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.

• @Dimitris Rizopoulos Dear Dimitris, any thought? Thanks – Abderrahim Nov 24 '19 at 12:46
• You should try the aucJM() function. – Dimitris Rizopoulos Nov 24 '19 at 17:47
• @DimitrisRizopoulos Thanks Dimitris. I believe that aucJM() cannot compare a joint to a standard survival model. It can be used to compare different joints models. The reason is that aucJM() uses a fixed starting time point (s) for all subjects. In my case, s is subject-specific and represent the time of the subject’s last observation. I also think that aucJM() is model-dependent because it uses the fitted joint model in order to deal with observations censored within the prediction period (i.e. Thoriz). Please excuse my long text but I am a big fun of JMbayes! Best wishes. Abderrahim – Abderrahim Nov 25 '19 at 13:16

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.