Research into joint modelling methods has grown substantially over recent years. Previous research has predominantly concentrated on the joint modelling of a single longitudinal outcome and a single time-to-event outcome. In clinical practice, the data collected will often be more complex, featuring multiple longitudinal outcomes and/or multiple, recurrent or competing event times. Harnessing all available measurements in a single model is advantageous and should lead to improved and more specific model predictions.
Notwithstanding the increased flexibility and better predictive capabilities, the extension of the classical univariate joint modelling framework to a multivariate setting introduces a number of technical and computational challenges. These include the high-dimensional numerical integrations required, modelling of multivariate unbalanced data, and proper estimation of standard errors. Consequently, software capable of fitting joint models to multivariate data is lacking. Building on recent methodological developments, we extend the classical joint model to multiple continuous longitudinal outcomes, and describe how to fit it using a Monte Carlo Expectation-Maximization algorithm with antithetic simulation for variance reduction. The development of a new R package, joineRML (https://CRAN.R-project.org/package=joineRML), will be demonstrated and an application to a clinical trial dataset will be presented.