  Joint models for longitudinal and survival data have gained a lot of
attention in recent years, with the development of myriad extensions to the
basic model, including those which allow for multivariate longitudinal data,
competing risks and recurrent events. Several software packages are now also
available for their implementation. Although mathematically straightforward,
the inclusion of multiple longitudinal outcomes in the joint model remains
computationally difficult due to the large number of random effects required,
which hampers the practical application of this extension. We present a novel
approach that enables the fitting of such models with more realistic
computational times. The idea behind the approach is to split the estimation of
the joint model in two steps; estimating a multivariate mixed model for the
longitudinal outcomes, and then using the output from this model to fit the
survival submodel. So called two-stage approaches have previously been
proposed, and shown to be biased. Our approach differs from the standard
version, in that we additionally propose the application of a correction
factor, adjusting the estimates obtained such that they more closely resemble
those we would expect to find with the multivariate joint model. This
correction is based on importance sampling ideas. Simulation studies show that
this corrected-two-stage approach works satisfactorily, eliminating the bias
while maintaining substantial improvement in computational time, even in more
difficult settings.
