Multivariate Joint Models for Multiple Longitudinal Outcomes and a Time-to-Event

A new, development version, of package JMbayes has been rolled-out on the dedicated GitHub repo. The major addition in this version is a set of new functions that can fit multivariate joint models for multiple longitudinal outcomes and a time-to-event. As with other GitHub packages, it can be easily installed using the install_github() function from package devtools. The new function that fits multivariate joint models is called mvJointModelBayes() and has a very similar syntax as the jointModelBayes() function described in a previous post. However, contrary to jointModelBayes() that is entirely written in R, the main bulk of computations of mvJointModelBayes() are based on C++ code building upon the excellent Rcpp and RcppArmadillo packages.

To make the connection between the two functions from a practical viewpoint more clear, let’s fit first a univariate joint model using mvJointModelBayes(). As for jointModelBayes(), we first need to separately fit a mixed-effects model for the longitudinal outcome, and a Cox model for the survival outcome. However, contrary to jointModelBayes(), which required the mixed model to be fitted with function lme() from the nlme package, for mvJointModelBayes() the mixed model needs to be fitted with the new function mvlgmer(). This follows the syntax of lmer() from package lme4. For example, the code below fits a linear mixed model for the longitudinal biomarker serum bilirubin from the Mayo Clinic PBC data set:



The main arguments of this function are: formulas a list of lme4-like formulas (a formula per outcome), data a data.frame that contains all the variables specified in formulas (NAs allowed), and families a list of family objects specifying the type of each outcome (currently the gaussian, binomial and poisson families are allowed). Following, we fit a Cox model using function coxph() from the survival package; argument model of coxph() needs to be set to TRUE:


Then the univariate joint model is fitted by the following simple call to mvJointModelBayes():


To fit a multivariate joint model for multiple longitudinal outcomes, the only thing that needs to be adapted from the procedure described above, is the call to mvglmer(). For example, to fit a bivariate joint model for serum bilirubin and presence of spiders, we use the syntax:


In addition, similarly to jointModelBayes(), mvJointModelBayes() allows to specify different functional forms for the longitudinal outcomes that are included in the Cox model.  As an example, we extend model JMFit2 by including the value and slope term for bilirubin, and the area term for spiders (in the log-odds scale). To include these terms into the model we specify the Formulas argument. This is specified in a similar manner as the derivForms argument of jointModelBayes(). In particular, it should be a list of lists. Each component of the outer list should have as name the name of the corresponding outcome variable. Then in the inner list we need to specify four components, namely, the fixed & random R formulas specifying the fixed and random effects part of the term to be included, and indFixed & indRandom integer indices specifying which of the original fixed and random effects are involved in the calculation of the new term. In the inner list you can also optionally specify a name for the term you want to include. A couple of notes:

1. For terms not specified in the Formulas list, the default value functional form is used.
2. If for a particular outcome you want to include both the value functional form and an extra term, then you need to specify that in the Formulas using two entries. To include the value functional form you only need to set the corresponding to “value”, and for the second term to specify the inner list. See example below on how to include the value and slope for serum bilirubin (for example, if the list below the entry “log(serBilir)” = “value” was not given, then only the slope term would have been included in the survival submodel).

We further extend the previous model by including the interactions terms between the terms specified in Formulas and the randomized treatment indicator drug. The names specified in the list that defined the interaction factors should match with the names of the output from JMFit3. Because of the many association parameters we have, we place a shrinkage prior on the association coefficients ‘alpha’. In particular, if we have K association parameters, we assume that alpha_k ~ N(0, tau * phi_k), k = 1, …, K. The precision parameters tau and phi_k are given Gamma priors. Parameter tau is global shrinkage parameter, and phi_k a specific shrinkage parameter per alpha coefficient:



Future updates of the package will include additional features, including among others

• a method for the generic functions survfitJM() and predict() to compute dynamic predictions for multivariate longitudinal data setting for both the longitudinal and survival outcomes.
• allow for competing risks.
• make the C++ code (even) faster.