Joint Models with Multiple Longitudinal Outcomes and a Time-to-Event Outcome: a Corrected Two-Stage Approach

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.

Joint models for longitudinal and survival data have become a valuable asset in the toolbox of modern data scientists. After the seminal papers of Faucett and Thomas [1996] and Wulfsohn and Tsiatis [1997], several extensions of these models have been proposed in the literature. These include, amongst others, flexible specification of the longitudinal model [Brown et al., 2005a], consideration of competing risks [Elashoff et al., 2008, Andrinopoulou et al., 2014 and multi-state models [Ferrer et al., 2016], and the calculation of dynamic predictions [Proust-Lima and Taylor, 2009, Rizopoulos, 2011, Andrinopoulou and Rizopoulos, 2016, Rizopoulos et al., 2017, Andrinopoulou et al., 2018. A particularly useful and practical extension is that which allows for the inclusion of multiple longitudinal outcomes [Rizopoulos and Ghosh, 2011a, Chi and Ibrahim, 2006, Brown et al., 2005b, Lin et al., 2002. In medical settings in particular, data collection is likely to be complex: while the standard joint model allows us to determine the association between a survival outcome and a single longitudinal outcome (biomarker), there are more often than not multiple biomarkers relevant to the event of interest. Extending the univariate joint model to accommodate these multiple longitudinal outcomes allows us to incorporate more information, improving prognostication and enabling us to better make sense of the complex underlying nature of the disease dynamics. A motivating example of this is the Bio-SHiFT cohort study; a prospective observational study conducted in the Netherlands on chronic heart failure (CHF) patients. The primary focus of the study was to determine whether or not disease progression in individual CHF patients can be assessed using longitudinal measurements of several blood biomarkers [van Boven et al., 2018]. Previous work on this data has focused mainly on the association between each individual biomarker and a single composite event, but it is likely that the predictive value of the biomarkers will be more accurately determined when they are assessed in concert.
Extension to the multivariate case is mathematically straightforward, and may be easily combined with other extensions, allowing for longitudinal outcomes of varying types; left, right and interval censoring; and the inclusion of competing risks, amongst others. There are also now a number of excellent software packages available, which make for easier implementation of the more complex models. There are however technical challenges which hamper the widespread use of these models. As the number of longitudinal outcomes increases, and thus the number of random effects, standard methods become computationally prohibitive: under a Bayesian approach, the number of parameters to sample becomes unreasonably large, and in the case of maximum likelihood, we are required to numerically approximate the integrals over the random effects, which is challenging in high dimensions. The practical solution most commonly used in such settings is that of the two-stage approach, wherein a multivariate mixed model is first used for the longitudinal outcomes, following which, the output of this model is used to fit a survival submodel. Unfortunately, substantial research on this topic indicates that this approach results in biased estimates [Tsiatis and Davidian, 2004, Rizopoulos, 2012, Ye et al., 2008. In this paper, we propose an adaptation of the simple two-stage approach which eliminates the bias and substantially reduces computational time. We propose the use of a correction factor, based on importance sampling theory [Press et al., 2007, Section 7.9]. This correction factor allows us to re-weight each realization of the MCMC sample obtained from the Bayesian estimation of the twostage approach, such that the resulting estimates more closely approximate those obtained via the full multivariate joint model. The weights are given by the target distribution (the full posterior distribution of the multivariate joint model), divided by the product of the posterior distributions for each of the two stages, evaluated for each iteration of the MCMC sample. The use of this correction factor alone is not enough to eliminate the bias, but, prior to its application, the two-stage approach is itself modified: where before, in the second stage, only the parameters of the survival submodel were updated, we now also update the random effects. These adaptations combined, achieve unbiased estimates in a fraction of the time required to compute the full multivariate model. The rest of the paper is organised as follows: Section 2 introduces the full multivariate joint model, and Section 3 discusses the estimation of the model under the Bayesian paradigm. Section 4 introduces the importance-sampling corrected two-stage approach, and presents the results of a simple simulation, and Section 5 the importance-sampling corrected two-stage approach with updated random effects. Section 6 presents the results of a more complex simulation, and finally in Section 7 we look at an analysis of the Bio-SHiFT data.

Joint Model Specification
We start with a general definition of the framework of multivariate joint models for multiple longitudinal outcomes and an event time. Let D n = {T i , T U i , δ i , y i ; i = 1, . . . , n} denote a sample from the target population, where T * i denotes the true event time for the i-th subject, T i and T U i the observed event times. Then δ i ∈ {0, 1, 2, 3} denotes the event indicator, with 0 corresponding to right censoring (T * i > T i ), 1 to a true event (T * i = T i ), 2 to left censoring (T * i < T i ), and 3 to interval censoring (T i < T * i < T U i ). Assuming K longitudinal outcomes we let y ki denote the n ki × 1 longitudinal response vector for the k-th outcome (k = 1, . . . , K) and the i-th subject, with elements y kij denoting the value of the k-th longitudinal outcome for the i-th subject, taken at time point t kij , j = 1, . . . , n ki . To accommodate multivariate longitudinal responses of different types in a unified framework, we postulate a generalized linear mixed effects model.
In particular, the conditional distribution of y ki given a vector of random effects b ki is assumed to be a member of the exponential family, with linear predictor given by where g k (·) denotes a known one-to-one monotonic link function, y ki (t) denotes the value of the k-th longitudinal outcome for the i-th subject at time point t, and x ki (t) and z ki (t) denote the design vectors for the fixed-effects β k and the random effects b ki , respectively. The dimensionality and composition of these design vectors is allowed to differ between the multiple outcomes, and they may also contain a combination of baseline and time-varying covariates. To account for the association between the multiple longitudinal outcomes we link their corresponding random effects. More specifically, the complete vector of random ef- is assumed to follow a multivariate normal distribution with mean zero and variance-covariance matrix D.
For the survival process, we assume that the risk for an event depends on a function of the subject-specific linear predictor η i (t) and/or the random effects. More specifically, we have where H ki (t) = {η ki (s), 0 ≤ s < t} denotes the history of the underlying longitudinal process up to t, h 0 (·) denotes the baseline hazard function, and w i (t) is a vector of exogenous, possibly time-varying, covariates with corresponding regression coefficients γ. Functions f kl (·), parameterized by vector α kl , specify which components/features of each longitudinal outcome are included in the linear predictor of the relative risk model Brown [2009], Rizopoulos and Ghosh [2011b], Rizopoulos [2012], . Some examples, motivated by the literature, are (subscripts kl have been dropped in the following expressions but are assumed): These formulations of f (·) postulate that the hazard of an event at time t may be associated with the underlying level of the biomarker at the same time point, the slope of the longitudinal profile at t or the accumulated longitudinal process up to t. In addition, the specified terms from the longitudinal outcomes may also interact with some covariates in the w i (t). Furthermore, note, that we allow a combination of L k functional forms per longitudinal outcome. Finally, the baseline hazard function h 0 (·) is modeled flexibly using a B-splines approach, i.e., where B q (t, v) denotes the q-th basis function of a B-spline with knots v 1 , . . . , v Q and γ h0 the vector of spline coefficients; typically Q = 15 or 20.

Likelihood and Priors
As explained in Section 1, we use a Bayesian approach for the estimation of the joint model's parameters.
The posterior distribution of the model parameters given the observed data is derived under the assumptions that given the random effects, the longitudinal outcomes are independent from the event times, the multiple longitudinal outcomes are independent of each other, and the longitudinal responses of each subject in each outcome are independent. Under these assumptions the posterior distribution is analogous to: where θ denotes the full parameter vector, and with ψ kij (b ki ) and ϕ denoting the natural and dispersion parameters in the exponential family, respectively, and c k (·), a k (·), and d k (·) are known functions specifying the member of the exponential family. For the survival part accordingly we have where I(·) denotes the indicator function. The integral in the definition of the cumulative hazard function does not have a closed-form solution, and thus a numerical method is employed for its evaluation. Standard options are the Gauss-Kronrod and Gauss-Legendre quadrature rules. For the parameters of the longitudinal outcomes we use standard default priors. The covariance matrix of the random effects is parameterized in terms of a correlation matrix Ω and a vector of σ d . For the correlation matrix Ω we use the LKJ-Correlation prior proposed by Lewandowski et al. [2009] with parameter ζ = 1.5. For each element of σ d we use a half-Student's t prior with 3 degrees of freedom. For the regression coefficients γ of the relative risk model we assume independent normal priors with zero mean and variance 1000. The same prior is also assumed for the vector of association parameters α. However, when α becomes high dimensional (e.g., when several functional forms are considered per longitudinal outcome), we opt for a global-local ridge-type shrinkage prior. More specifically, for the s-th element of α we assume αs ∼ N (0, τ ψs), The global smoothing parameter τ has sufficient mass near zero to ensure shrinkage, while the local smoothing parameter ψ s allows individual coefficients to attain large values. The motivation for using this type of prior distribution in this case is that we expect the different terms behind the specification of f (·) to be correlated, and many of the corresponding coefficients to be non-zero. Nonetheless, other options of shrinkage or variableselection priors could also be used Andrinopoulou and Rizopoulos [2016]. Finally, the penalized version of the B-spline approximation to the baseline hazard is specified using the following hierarchical prior for γ h0 Lang and Brezger [2004]: where τ h is the smoothing parameter that takes a Gamma(1, τ hδ ) prior distribution, with a hyper-prior This limitation has hampered the use of multivariate joint models in practice. The two-stage approach, which entails fitting the longitudinal and survival outcomes separately, is the solution most often used to overcome this computational deadlock. Using this approach, under the Bayesian framework, we would have the following two stages: S-I: We fit a multivariate mixed model for the longitudinal outcomes using either MCMC or HMC, and we obtain a sample {θ (m) y , b (m) ; m = 1, . . . , M } of size M from the posterior, where θ y denotes the subset of the parameters that are included in the definition of the longitudinal submodels (including the parameters in the random-effects distribution).
S-II: Utilizing the sample from Stage I, we obtain a sample for the parameters of the survival submodel where θ t denotes the subset of the parameters that are included in the definition of the survival submodel, and T = (T, T U ). This two-stage procedure essentially entails the same number of iterations as the full Bayesian estimation of the multivariate joint model. The computational benefits stem from the fact that we do not need to numerically integrate the survival submodel density function in Stage I. Even though this approach greatly reduces the computational burden, there exists a substantial body of work demonstrating that it results in biased estimates, even in the simpler case of univariate joint models [see Davidian, 2004, Rizopoulos, 2012, andreferences therein]. This bias is a result of not working with the full joint distribution, which would produce estimates of θ y and b that are appropriately corrected for informative dropout relating to the occurrence of an event.
To overcome this issue, we propose the correction of the estimates we obtain from the two-stage approach using importance sampling weights [Press et al., 2007, Section 7.9]. In particular, we consider that the . . , M } that we have obtain using the two-stage approach can be considered a weighted sample from the full posterior of the multivariate joint model with weights given by: The numerator in this expression is the posterior distribution of the multivariate joint model, and the denominator, the corresponding posterior distributions from each of the two stages. As previously stated, from (6) we observe that the difference between fitting the full joint model versus the two-stage approach comes from the second term in the numerator and denominator. By expanding these two terms we obtain The resulting weights involve a marginal likelihood calculation, which we perform using a Laplace approximation, namely and q denotes the dimensionality of the θ t vector. The extra computational burden of performing this Laplace approximation is minimal in practice, since good initial values can be provided from one iteration m to the next m + 1, which substantially reduces the number of required optimization iterations for finding is provided as an initial value to find θ (m+1) t ).

Performance
To evaluate whether the introduction of the importance sampling weights alleviates the bias observed with the simple two-stage approach (i.e., without the weights), we perform a 'proof-of-concept' simulation study.
In particular, we compare the proposed corrected two stage approach with the simple two-stage approach, Comparing the first two of these approaches, we see that the calculation of the importance-sampling weights in the corrected two-stage approach had minimal computational cost, with the full multivariate joint model taking substantially more time to fit. Figure D.2 shows boxplots of posterior means from the 500 datasets for the parameters of the two longitudinal submodels. We observe that all three approaches provide very similar results with minimal bias. Figure 1 shows the corresponding boxplots of posterior means for the parameters of the survival submodel. As expected, the full multivariate joint model returns unbiased results. Similarly, as has previously been reported in the literature, the simple two-stage approach exhibits considerable bias.
We see that this bias persists for the corrected two-stage approach, although theoretically the use of the importance sampling weights should alleviate it (by adjusting the posterior means obtained via the simple two-stage approach such that they more closely resemble those from the full multivariate model).
5 Corrected Two-Stage Approach with Random Effects 5.1 Importance sampling correction with random effects The above result is unexpected, since (as per Figure D.2), the corrected two-stage (and indeed the simple twostage) approach unbiasedly estimates both the fixed effects and the variance components of the longitudinal submodels. However, further investigation shows that there is a considerable difference between the corrected two-stage approach and the multivariate joint model with regards to the posterior of the random effects.
This is depicted in Figure 2 for one of the longitudinal outcomes we have simulated.
[ Figure 2 about here.] The data have been simulated such that higher values for longitudinal outcome y 1 are associated with a higher hazard of the event. From Figure 2 we observe that the random effect estimates for the multivariate mixed model, and especially the random slope estimates for subjects with and without an event differ from those for the multivariate joint model. In particular, we observe that the random slope estimates from the joint model are larger for subjects with an event compared to the linear mixed model, and vice versa for subjects without an event. This observation suggests that we could improve the weights given in (6) by updating (in the second stage) not only the parameters of the survival submodel θ t but also the random effects b. That is, we obtain a sample for the parameters of the survival submodel {θ from the corresponding joint posterior distribution, , the corresponding second stage presented in Section 4, since we now also need to calculate the densities of the mixed-effect models for the K longitudinal outcomes. Nonetheless, the computational gains compared to fitting the full joint model remain significant. Under this second stage (8) the importance sampling weights now take the form: Similarly to (6), the new weights have been formulated such that the difference lies in the second term in both the numerator and denominator. By doing an expansion of these two terms similar to that used in the previous section, we obtain: and the self-normalized weights arew The integrals in the numerator are once again approximated using the Laplace method, namely, we let denote the Hessian matrix for the random effects, and analogously, denote the Hessian matrix for the θ t parameters. Then, we approximate the inner integral by where κ denotes the number of random effects for each subject i. Similarly, the outer integral is approximated Given the requirement for a double Laplace approximation, and the fact that the denominator does not simplify, the calculation of the ̟ (m) weights given by (10) is more computationally intensive than the ones presented in Section 4. Nevertheless, these required computations still remain many orders of magnitude faster than fitting the full joint model.

Performance
To assess whether updating the random effects in the importance sampling weights alleviates the bias we observed in Section 4.2, we have re-analyzed the same simulated datasets. The details are again given in Appendix A.1. The results from 500 simulated datasets are presented in Figures 3, D.1, and D.3. As anticipated, the corrected two-stage approach with updated random effects added only a small computational cost, with the full multivariate joint model still taking considerably more time to fit than either of the corrected two-stage approaches ( Figure D.1). The boxplots depicting the posterior means from the 500 datasets for the parameters of the longitudinal submodels once again demonstrate similar results for all three approaches (Figure D.3). Figure 3 shows the posterior means for the parameters of the survival submodel. We observe that the bias seen for the corrected two-stage approach is now eliminated, with the posterior means from the approach with updated random effects closely approximating those from the full multivariate joint model.

Extra Simulations
Further simulations were performed in order to assess the performance of the importance-sampling-corrected two-stage approach with the updated random effects, in different scenarios. Details of these simulations are given in Appendices A.2 and A.3.

Scenario II
Scenario II included 6 continuous longitudinal outcomes. Owing to the increased number of outcomes, the full multivariate joint model was not run. Table 1 shows the bias for the parameters of the survival submodel, together with the RMSE and coverage (based on the 2.5% and 97.5% credibility intervals for each parameter). Table C.1 in the appendix shows the same information for the parameters of the 6 longitudinal outcomes. [

Scenario III
Scenario III again included 6 longitudinal outcomes, now of varying types: 3 continuous and 3 binary. Table 2 demonstrates yet again the alleviation of the bias achieved by updating the random effects.
[ hospitalisation for heart failure, cardiac death, LVAD placement and heart transplantation). Biomarkers were measured at inclusion and subsequently every 3 months until the end of follow-up. We focus on 6 biomarkers: the glomerular marker Cystatin C (CysC), two tubular markers; urinary N-acetyl-beta-D-glucosaminidase (NAG) and kidney-injury-molecule (KIM)-1, and the markers N-terminal propBNP (NT-proBNP), cardiac troponin T (HsTNT), and C-reactive protein (CRP). The latter three markers are known to be related to poor outcomes in CHF patients, and measure various aspects of heart failure pathophysiology (wall stress, myocyte damage and inflammation respectively). All biomarkers were logarithmically transformed for further analysis (log base 2) due to skewness. For each of NT-probnp, HsTNT and CRP, we included natural cubic splines in both the fixed and random effects parts of their longitudinal models, with differing numbers of knots per outcome ( Figure D.4). Simple linear models with random intercept and slope were used for CysC, NAG and KIM-1. Thus, for each of CysC, NAG and KIM-1 (k = 1, 2, 3), we fit: For the remaining outcomes (k = 4, 5, 6), we have: where B kn (t; λ p ) denotes the B-spline basis matrix for a natural cubic spline of time with two internal knots placed at the 25th and 75th percentiles of the follow up times for NT-probnp (p = 1, 2, 3), and one internal knot placed at the 50th percentile of the follow up times for each of HsTNT and CRP (p = 1, 2).
Boundary knots were set at the 5th and 95th percentiles. We assume a multivariate normal distribution for the random effects, b We fit 3 joint models, using the global-local ridge-type shrinkage prior previously described in each case. Model 1 included only the current underlying value of the longitudinal marker for each of the 6 markers. Model 2 included the current value and slope for each marker, and model 3 included the integrated longitudinal profile for each marker (AUC). We thus have: The parameter estimates and 95% credibility intervals for the event process are presented in Table 3. Hazard ratios are presented per doubling of level, slope or AUC at any point in time.

Discussion
In this paper, we presented a novel approach for fitting joint models which allows for the inclusion of multivariate longitudinal outcomes with realistic computing times. We demonstrated once again, the bias of the estimated parameters for the survival process characteristic of the standard two-stage approach, and proposed the use of an importance-sampling corrected two-stage approach, with updated random effects, in its place. Our approach was shown to be successful, producing satisfactory results in a number of simulation scenarios: both survival and longitudinal estimates were unbiased, and computing times were reduced by several orders of magnitude, compared to the full multivariate joint model. We were easily able to incorporate multiple outcomes in the analysis of the Bio-SHiFT data, obtaining very similar results to those previously noted for the CHF-related biomarkers (CRP, HsTNT and NT-proBNP). We did not find any significant associations between any of the renal markers (CysC, NAG and KIM-1) and the risk of the composite event in the multivariate analysis, indicating that their predictive value may not be independent of the CHF-related markers. While the simulations included up to 6 multiple outcomes of varying types, it would be interesting to confirm our results in even more complex settings, (perhaps incorporating competing risks such as those present in the Bio-SHiFT study), and to try determine the limits of the methodology. A further topic for research would be methods for increasing the speed of computation involved in fitting the multivariate mixed model itself, so as to extend the number of outcomes even further. The proposed importance-sampling corrected two-stage estimation approach is implemented in function mvJointModelBayes() in the freely available package JMbayes (version 0.8-0) for the R programming language (freely available from the Comprehensive R Archive Network). An example of how these functions should be used can be found in the appendices.     D). The variance-covariance matrix D has general form:

Posterior
and for Scenario I: Time was simulated from a uniform distribution between 0 and 25. For the survival outcome, adjusting for group allocation, we used: The baseline risk was simulated from a Weibull distribution h 0 (t) = φt φ−1 , with φ = 1.65. For the simulation of the censoring times, an exponential censoring distribution was selected, with mean µ = 15, such that the censoring rate was between 60% and 70%. More details are presented in Table C.3.

A.2 Scenario II
Scenario II is an extension of Scenario I, such that we now have K = 6 continuous longitudinal outcomes. We again simulate 500 patients with a maximum of 15 repeated measurements per patient. The k longitudinal outcomes each had form: for group allocation as in Scenario I, we used: The baseline risk was simulated using B-splines with knots specified apriori. An exponential censoring distribution was used for the simulation of the censoring times, with mean µ = 15, such that the censoring rate was between 60% and 70%. Further details are again available in Table C.3.

A.3 Scenario III
In Scenario III, we simulate 500 patients with a maximum of 15 repeated measurements per patient, including 3 continuous and 3 binary longitudinal outcomes, such that: where g k (·) denotes the canonical link function appropriate to the response type (identity and logit for the gaussian and binomial outcomes respectively), and b Scenario III maintains the use of the uniform distribution between 0 and 25 for time, and the use of B-splines for the simulation of the baseline hazard. The censoring times were simulated using an exponential censoring distribution as before, with mean µ = 15. Table C.3 provides additional information.

B Example R Code
The below code fits a multivariate joint model for K = 3 longitudinal outcomes: y 1 , y 2 and y 3 , where y 1 is binary and both y 2 and y 3 are continuous. We fit a linear mixed model for y 1 with random intercept and slope (time is f utime), and use natural cubic splines with two knots, (at f utime = 6 and f utime = 15 respectively) in both the fixed and random parts of the models for y 2 and y 3 . The survival submodel adjusts for continuous baseline predictors x 1 and x 2 .