LOOIC for joint models?

[danielinteractive]

Hi @gowerc and @gravesti ,

I am not sure if I am doing something wrong, but while I am able to calculate the LOOIC criterion with the $loo() method for CmdStanMCMC objects when I either have a longitudinal or a survival model from jmpost, I cannot do so when I have a true joint model with both. I get then as error:

> os_looic <- mcmc_joint_results$loo(r_eff = FALSE)
Error: Can't find the following variable(s) in the output: log_lik

Is this expected?

Thanks,
cheers
Daniel

[danielinteractive]

@mercifr1 fyi

[gowerc]

Hey @danielinteractive

Yes this is expected. I'm struggling to remember the exact argument / issue that led to this but the core of it was that "an observation" isn't well defined for a joint model as loo on the OS model means leaving out a whole subject but this corresponds to multiple observations on the TGI model

[gowerc]

Related issues (though it doesn't look like I explained much more than what I've already put here):

• Compare loglik values with brms #391
• Update loo implementation to be for individual sub models only #402

[danielinteractive]

Thanks @gowerc for the quick reply,

I see - I mean in principle it would be well defined if we defined as one observation being a whole subject (including its OS and TGI data points). I guess the question is then whether this works correctly using the Stan mechanics. In the past, I remember having calculated the LOOIC based on the subject level log-likelihood contributions, and it did work ok in the sense that the LOOIC was decomposed into the TGI and OS LOOIC contributions in an additive way.

[gravesti]

@danielinteractive To that end, I think it's possible to do this with after the sampling with generated quantities:

generated quantities {
    vector[n_subjects] log_lik;

    for (i in 1:n_subjects) {
        // survival data, one per patient
        log_lik[i] = os_subj_log_lik[i];
    }
    for (i in 1:n_tumour_all) {
        // longitudinal data, multiple per patient
        log_lik[subject_tumour_index[i]] += long_obvs_log_lik[i];
    }
}

and then add this to the model

library(jmpost)
library(cmdstanr)
library(loo)
mp <- sampleStanModel(...)
mc <- as.CmdStanMCMC.JointModelSamples(mp)
tmp <- tempfile(fileext = ".stan") # file extension for syntax highlighting
write_stan(jm, destination = tmp)
file.edit(tmp) # opens the Stan file in the default editor

# add generated quantities block and save

mod_gq <- cmdstan_model(tmp)
fit_gq <- mod_gq$generate_quantities(mc, data = mc$data_file(), seed = 123)

loo(fit_gq$draws("log_lik"), r_eff=relative_eff(exp(fit_gq$draws("log_lik"))))

[danielinteractive]

Thanks @gravesti , that could indeed work as a temporary workaround!

I think if we could add that to jmpost that would be even better 😉

[gowerc]

Thanks @gowerc for the quick reply,

I see - I mean in principle it would be well defined if we defined as one observation being a whole subject (including its OS and TGI data points). I guess the question is then whether this works correctly using the Stan mechanics. In the past, I remember having calculated the LOOIC based on the subject level log-likelihood contributions, and it did work ok in the sense that the LOOIC was decomposed into the TGI and OS LOOIC contributions in an additive way.

I guess this is where my limitation of the maths/methods comes in; As in, given that what is being calculated is an estimate of the LOO values using Pareto Smoothed Importance Sampling rather than exact values I wasn't sure if the methods generalised / were valid for the case of removing multiple values from the input dataset(s). For example I don't think I've seen any articles saying this method can be used for K-fold cross-validation which is essentially what we'd be doing. I stress though that this is just me not understanding the method and it may all be fine.

[gravesti]

Ok, thanks for bringing that up @gowerc . From my reading of the cross validation FAQ, it wouldn't be a problem per se to leave out multiple observations, but what is a problem is leaving out all observations that contribute to a parameter. With patient specific random effects this will be a problem for us.

In that case they propose two options:

• if it's only 1 or 2 parameters integrate them out. <-- we have much more than 2 patient specific parameters
• do the usual K-fold cross-validation.

[danielinteractive]

Makes sense, however also when we go for K-fold CV we would still benefit from having the log_lik available in the model samples already. So I think that is still very helpful