Multivariate Joint Models - Theory

Joint Model Specification

Joint models for longitudinal and survival data have gained a lot of attention the recent years. There have been extended to handle among others multivariate longitudinal data, competing risks and recurrent events, and nowadays there also exist several freely available software packages for their implementation. From the aforementioned extensions, the one that is most practically relevant is the multivariate longitudinal data one. Even though this extension is mathematically straightforward, from a computational viewpoint joint models with multiple longitudinal outcomes remain difficult to fit in practice due to the high number of random effects they require. This difficulty has also hampered to a degree their practical application. The current version of JMbayes implements a novel approach that enables fitting such joint models in realistic computing times.

We start with a general definition of the framework of multivariate joint models for multiple longitudinal outcomes and an event time. Let \(\mathcal D_n = \{T_i, T_i^U, \delta_i, y_i; i = 1, \ldots, n\}\) denote a sample from the target population, where we let \(T_i^*\) denote the true event time for the \(i\)-th subject, \(T_i\) and \(T_i^U\) the observed event times, and \(\delta_i \in \{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_i^U\)). Assuming \(K\) longitudinal outcomes we let \(y_{ki}\) denote the \(n_{ki} \times 1\) longitudinal response vector for the \(k\)-th outcome (\(k = 1, \ldots, K\)) and the \(i\)-th subject, with elements \(y_{kij}\) denoting the value of the \(k\)-th longitudinal outcome taken at time point \(t_{kij}\), \(j = 1, \ldots, 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 \[g_k \bigl [ E \{ y_{ki}(t) \mid b_{ki} \} \bigr ] = \eta_{ki}(t) = x_{ki}^\top(t) \beta_k + z_{ki}^\top(t) b_{ki},\] where \(g_k(\cdot)\) denotes a known one-to-one monotonic link function, and \(y_{ki}(t)\) denotes the value of the \(k\)-th longitudinal outcome for the \(i\)-th subject at time point \(t\), \(x_{ki}(t)\) and \(z_{ki}(t)\) denote the design vectors for the fixed-effects \(\beta_k\) and for 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 mix 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 effects \(b_i = (b_{1i}, b_{2i}, \ldots, b_{Ki})^\top\) 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 \(\eta_i(t)\) and/or the random effects. More specifically, we have \[\begin{eqnarray} \nonumber h_i (t \mid \mathcal H_i(t), w_i(t)) & = & \lim_{\Delta t \rightarrow 0} \Pr \{ t \leq T_i^* < t + \Delta t \mid T_i^* \geq t, \mathcal H_i(t), w_i(t) \} \big / \Delta t, \quad t > 0\\ & = & h_0(t) \exp \biggl [\gamma^\top w_i(t) + \sum \limits_{k = 1}^K \sum \limits_{l = 1}^{L_k} f_{kl} \{\mathcal H_{ki}(t), w_i(t), b_{ki}, \alpha_{kl} \} \biggr], \end{eqnarray}\] where \(\mathcal H_{ki}(t) = \{ \eta_{ki}(s), 0 \leq s < t \}\) denotes the history of the underlying longitudinal process up to \(t\), \(h_0(\cdot)\) denotes the baseline hazard function, \(w_i(t)\) is a vector of exogenous, possibly time-varying, covariates with corresponding regression coefficients \(\gamma\). Functions \(f_{kl}(\cdot)\), parameterized by vector \(\alpha_{kl}\), specifies which components/features of each longitudinal outcome are included in the linear predictor of the relative risk model. Some examples, motivated by the literature, are (subscripts \(kl\) have been dropped in the following expressions but are assumed): \[\begin{eqnarray} f \{\mathcal H_i(t), w_i(t), b_i, \alpha \} & = & \alpha \eta_i(t),\\ f \{\mathcal H_i(t), w_i(t), b_i, \alpha \} & = & \alpha_1 \eta_i(t) + \alpha_2 \eta_i'(t), \mbox{ with } \eta_i'(t) = \frac{d\eta_i(t)}{dt},\\ f \{\mathcal H_i(t), w_i(t), b_i, \alpha \} & = & \alpha \int_0^t \eta_i(s) \, ds. \end{eqnarray}\] These formulations of \(f(\cdot)\) 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(\cdot)\) is modeled flexibly using a B-splines approach, i.e., \[\log h_0(t) = \sum \limits_{q = 1}^Q \gamma_{h_0,q} B_q(t, v),\] where \(B_q(t, v)\) denotes the \(q\)-th basis function of a B-spline with knots \(v_1, \ldots, v_Q\) and \(\gamma_{h_0}\) the vector of spline coefficients. To avoid the task of choosing the appropriate number and position of the knots, we include a relatively high number of knots (e.g., 15 to 20) and appropriately penalize the B-spline regression coefficients \(\gamma_{h_0}\) for smoothness using the differences penalty.

Likelihood and Priors

For the estimation of joint model’s parameters we will use a Bayesian approach. 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: \[\begin{eqnarray} p(\theta, b) \propto \prod \limits_{i = 1}^n \prod \limits_{k = 1}^K \prod \limits_{j = 1}^{n_{ki}} p (y_{kij} \mid b_{ki}, \theta) \; p(T_i, T_i^U, \delta_i \mid b_{ki}, \theta) \; p(b_{ki} \mid \theta) \; p(\theta), \end{eqnarray}\] where \(\theta\) denotes the full parameter vector, and \[p (y_{kij} \mid b_{ki}, \theta) = \exp \biggl \{ \Bigl [y_{kij} \psi_{kij}(b_{ki}) - c_k\{\psi_{kij}(b_{ki})\} \Bigr ] \Big / a_k(\varphi) - d_k(y_{kij}, \varphi) \biggr \},\] with \(\psi_{kij}(b_{ki})\) and \(\varphi\) denoting the natural and dispersion parameters in the exponential family, respectively, \(c_k(\cdot)\), \(a_k(\cdot)\), and \(d_k(\cdot)\) are known functions specifying the member of the exponential family. For the survival part accordingly we have \[\begin{eqnarray} p(T_i, T_i^U, \delta_i \mid b_i, \theta) & = &\Bigr \{h_i (T_i \mid \mathcal H_i(T_i), w_i(T_i)) \Bigr \}^{I(\delta_i = 1)} \exp \biggl \{- \int_0^{T_i} h_i (s \mid \mathcal H_i(s), w_i(s)) \; ds \biggr \}^{I(\delta_i \in \{0, 1\})}\\ \nonumber & \times & \Biggl \{1 - \exp \biggl \{- \int_0^{T_i} h_i (s \mid \mathcal H_i(s), w_i(s)) \; ds \biggr \}\Biggr\}^{I(\delta_i = 2)}\\ & \times & \Biggl \{\exp \biggl \{- \int_0^{T_i} h_i (s \mid \mathcal H_i(s), w_i(s)) \; ds \biggr \} - \exp \biggl \{- \int_0^{T_i^U} h_i (s \mid \mathcal H_i(s), w_i(s)) \; ds \biggr \} \Biggr\}^{I(\delta_i = 3)}, \end{eqnarray}\] where \(I(\cdot)\) 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 must be 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. More specifically, independent normal priors with zero mean and variance 1000 for the fixed effects and inverse Gamma priors for scale parameters. The covariance matrix of the random effects is parameterized in terms of a correlation matrix \(\Omega\) and a vector of \(\sigma_{d}\). For the correlation matrix \(\Omega\) we use LKJ-Correlation prior proposed by with parameter \(\zeta = 1.5\). For each element of \(\sigma_{d}\) we use a half-Student’s t prior with 3 degrees of freedom. For the regression coefficients \(\gamma\) 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 \(\alpha\). However, when \(\alpha\) becomes high dimensional (e.g., when several functional forms are considered per longitudinal outcome), we opt for a global-local ridge-type shrinkage prior, i.e., for the \(s\)-th element of \(\alpha\) we assume \[\alpha_s \sim \mathcal N (0, \tau \psi_s), \quad \tau^{-1} \sim \mbox{Gamma}(0.1, 0.1), \quad \psi_s^{-1} \sim \mbox{Gamma}(1, 0.01).\] The global smoothing parameter \(\tau\) has sufficiently mass near zero to ensure shrinkage, while the local smoothing parameter \(\psi_s\) allows individual coefficients to attain large values. Other options of shrinkage or variable-selection priors could be used as well. Finally, the penalized version of the B-spline approximation to the baseline hazard is specified using the following hierarchical prior for \(\gamma_{h_0}\): \[p(\gamma_{h_0} \mid \tau_h) \propto \tau_h^{\rho(K)/2}\exp \Bigl (-\frac{\tau_{h}}{2} \gamma_{h_0}^\top K \gamma_{h_0} \Bigr ),\] where \(\tau_h\) is the smoothing parameter that takes a \(\mbox{Gamma}(1, \tau_{h\delta} )\) prior distribution, with a hyper-prior \(\tau_{h\delta} \sim \mbox{Gamma}(10^{-3}, 10^{-3})\), which ensures a proper posterior distribution for \(\gamma_{h_0}\), \(K = \Delta_r^\top \Delta_r + 10^{-6}I\), with \(\Delta_r\) denoting the \(r\)-th difference penalty matrix, and \(\rho(K)\) denotes the rank of \(K\).

Multivariate Joint Models - Practice

The basic multivariate joint model

The function that fits multivariate joint models in JMbayes is called mvJointModelBayes() and has a very similar syntax as the jointModelBayes() function. 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 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:

MixedModelFit1 <- mvglmer(list(log(serBilir) ~ year + (year | id)), data = pbc2,
                          families = list(gaussian))

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:

pbc2.id$Time <- pbc2.id$years
pbc2.id$event <- as.numeric(pbc2.id$status != "alive")
CoxFit <- coxph(Surv(Time, event) ~ drug + age, data = pbc2.id, model = TRUE)

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

JMFit1 <- mvJointModelBayes(MixedModelFit1, CoxFit, timeVar = "year")
summary(JMFit1)

## 
## Call:
## mvJointModelBayes(mvglmerObject = MixedModelFit1, survObject = CoxFit, 
##     timeVar = "year")
## 
## Data Descriptives:
## Number of Groups: 312        Number of events: 169 (54.2%)
## Number of Observations:
##   log(serBilir): 1945
## 
## Random-effects covariance matrix:
##       StdDev    Corr
## (I)1  1.0081  (Int)1
## yer1  0.1735  0.4016
## 
## Survival Outcome:
##                     PostMean  StDev  StErr    2.5%  97.5%     P
## drugD-penicil        -0.0580 0.1578 0.0039 -0.3632 0.2446 0.728
## age                   0.0430 0.0078 0.0003  0.0284 0.0578 0.000
## log(serBilir)_value   1.3054 0.0829 0.0022  1.1492 1.4789 0.000
## 
## Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
##             PostMean  StDev  StErr   2.5%  97.5% P
## (Intercept)   0.4956 0.0596 0.0019 0.3767 0.6096 0
## year          0.1764 0.0130 0.0004 0.1522 0.2026 0
## sigma         0.3489 0.0067 0.0002 0.3366 0.3627 0
## 
## MCMC summary:
## iterations: 900 
## burn-in: 600 
## thinning: 300 
## time: 2.6 min

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:

MixedModelFit2 <- mvglmer(list(log(serBilir) ~ year + (year | id),
                               spiders ~ year + (1 | id)), data = pbc2,
                          families = list(gaussian, binomial))

JMFit2 <- mvJointModelBayes(MixedModelFit2, CoxFit, timeVar = "year")
summary(JMFit2)

## 
## Call:
## mvJointModelBayes(mvglmerObject = MixedModelFit2, survObject = CoxFit, 
##     timeVar = "year")
## 
## Data Descriptives:
## Number of Groups: 312        Number of events: 169 (54.2%)
## Number of Observations:
##   log(serBilir): 1945
##   spiders: 1887
## 
## Random-effects covariance matrix:
##       StdDev    Corr        
## (I)1  1.0066  (Int)1   year1
## yer1  0.1723  0.4135        
## (I)2  3.0810  0.4736  0.4288
## 
## Survival Outcome:
##                     PostMean  StDev  StErr    2.5%  97.5%     P
## drugD-penicil        -0.0985 0.1536 0.0035 -0.3728 0.2178 0.516
## age                   0.0422 0.0079 0.0002  0.0267 0.0577 0.000
## log(serBilir)_value   1.2678 0.0957 0.0024  1.0879 1.4629 0.000
## spiders_value         0.0562 0.0435 0.0014 -0.0280 0.1376 0.188
## 
## Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
##             PostMean  StDev  StErr   2.5%  97.5% P
## (Intercept)   0.4967 0.0592 0.0019 0.3795 0.6094 0
## year          0.1786 0.0129 0.0004 0.1533 0.2037 0
## sigma         0.3495 0.0068 0.0002 0.3360 0.3649 0
## 
## Longitudinal Outcome: spiders (family = binomial, link = logit)
##             PostMean  StDev  StErr    2.5%   97.5% P
## (Intercept)  -1.6716 0.2212 0.0150 -2.1535 -1.2896 0
## year          0.1830 0.0303 0.0032  0.1234  0.2408 0
## 
## MCMC summary:
## iterations: 900 
## burn-in: 600 
## thinning: 300 
## time: 3.9 min

Different Functional Forms

Function 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:

For terms not specified in the Formulas list, the default value functional form is used.
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).

Forms <- list("log(serBilir)" = "value",
              "log(serBilir)" = list(fixed = ~ 1, random = ~ 1,
                                     indFixed = 2, indRandom = 2, name = "slope"),
              "spiders" = list(fixed = ~ 0 + year + I(year^2/2), random = ~ 0 + year,
                               indFixed = 1:2, indRandom = 1, name = "area"))

JMFit3 <- update(JMFit2, Formulas = Forms)
summary(JMFit3)

## 
## Call:
## mvJointModelBayes(mvglmerObject = MixedModelFit2, survObject = CoxFit, 
##     timeVar = "year", Formulas = Forms)
## 
## Data Descriptives:
## Number of Groups: 312        Number of events: 169 (54.2%)
## Number of Observations:
##   log(serBilir): 1945
##   spiders: 1887
## 
## Random-effects covariance matrix:
##       StdDev    Corr        
## (I)1  1.0066  (Int)1   year1
## yer1  0.1723  0.4135        
## (I)2  3.0810  0.4736  0.4288
## 
## Survival Outcome:
##                     PostMean  StDev  StErr    2.5%  97.5%     P
## drugD-penicil        -0.1139 0.1677 0.0041 -0.4455 0.2137 0.490
## age                   0.0430 0.0079 0.0002  0.0267 0.0583 0.000
## log(serBilir)_value   1.2288 0.1106 0.0023  1.0125 1.4473 0.000
## log(serBilir)_slope   2.0192 0.9277 0.0293  0.2282 3.8321 0.034
## spiders_area         -0.0077 0.0078 0.0002 -0.0237 0.0072 0.306
## 
## Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
##             PostMean  StDev  StErr   2.5%  97.5% P
## (Intercept)   0.4967 0.0592 0.0019 0.3795 0.6094 0
## year          0.1786 0.0129 0.0004 0.1533 0.2037 0
## sigma         0.3495 0.0068 0.0002 0.3360 0.3649 0
## 
## Longitudinal Outcome: spiders (family = binomial, link = logit)
##             PostMean  StDev  StErr    2.5%   97.5% P
## (Intercept)  -1.6716 0.2212 0.0150 -2.1535 -1.2896 0
## year          0.1830 0.0303 0.0032  0.1234  0.2408 0
## 
## MCMC summary:
## iterations: 900 
## burn-in: 600 
## thinning: 300 
## time: 4.2 min

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 place thre shrinkage prior on the association coefficients described above. In particular:

Ints <- list("log(serBilir)" = ~ drug, "log(serBilir)_slope" = ~ drug,
             "spiders_area" = ~ drug)

JMFit4 <- update(JMFit3, Interactions = Ints, priors = list(shrink_alphas = TRUE))
summary(JMFit4)

## 
## Call:
## mvJointModelBayes(mvglmerObject = MixedModelFit2, survObject = CoxFit, 
##     timeVar = "year", Formulas = Forms, Interactions = Ints, 
##     priors = list(shrink_alphas = TRUE))
## 
## Data Descriptives:
## Number of Groups: 312        Number of events: 169 (54.2%)
## Number of Observations:
##   log(serBilir): 1945
##   spiders: 1887
## 
## Random-effects covariance matrix:
##       StdDev    Corr        
## (I)1  1.0066  (Int)1   year1
## yer1  0.1723  0.4135        
## (I)2  3.0810  0.4736  0.4288
## 
## Survival Outcome:
##                                   PostMean  StDev  StErr    2.5%  97.5%
## drugD-penicil                      -0.1388 0.3267 0.0109 -0.8010 0.4682
## age                                 0.0412 0.0065 0.0002  0.0281 0.0535
## log(serBilir)_value                 1.2508 0.1168 0.0053  1.0251 1.4756
## log(serBilir)_value:drugD-penicil  -0.0072 0.1352 0.0045 -0.2803 0.2643
## log(serBilir)_slope                 0.7845 0.8128 0.0578 -0.1999 2.6384
## log(serBilir)_slope:drugD-penicil   0.3660 0.6454 0.0204 -0.4409 2.2402
## spiders_area                       -0.0126 0.0102 0.0003 -0.0323 0.0079
## spiders_area:drugD-penicil          0.0110 0.0131 0.0004 -0.0159 0.0382
##                                       P
## drugD-penicil                     0.662
## age                               0.000
## log(serBilir)_value               0.000
## log(serBilir)_value:drugD-penicil 0.970
## log(serBilir)_slope               0.312
## log(serBilir)_slope:drugD-penicil 0.540
## spiders_area                      0.232
## spiders_area:drugD-penicil        0.374
## 
## Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
##             PostMean  StDev  StErr   2.5%  97.5% P
## (Intercept)   0.4967 0.0592 0.0019 0.3795 0.6094 0
## year          0.1786 0.0129 0.0004 0.1533 0.2037 0
## sigma         0.3495 0.0068 0.0002 0.3360 0.3649 0
## 
## Longitudinal Outcome: spiders (family = binomial, link = logit)
##             PostMean  StDev  StErr    2.5%   97.5% P
## (Intercept)  -1.6716 0.2212 0.0150 -2.1535 -1.2896 0
## year          0.1830 0.0303 0.0032  0.1234  0.2408 0
## 
## MCMC summary:
## iterations: 900 
## burn-in: 600 
## thinning: 300 
## time: 4.6 min

Time-varying Effects

A common feature of all aforementioned models is that the regression parameters \(\alpha\) that measure the strength of the association between the longitudinal and survival outcomes were assumed to be constant in time. Nonetheless, in some settings it may be plausible that the associationg parameters \(\alpha\) change during follow-up. To accomodate such settings mvJointModelBayes() allows modeling the functions \(\alpha(t)\) using P-splines. In particular, the definition of the risk model becomes: \[\begin{eqnarray} h_i (t \mid \mathcal H_i(t), w_i(t)) & = & h_0(t) \exp \biggl [\gamma^\top w_i(t) + \sum \limits_{k = 1}^K \sum \limits_{l = 1}^{L_k} f_{kl} \{\mathcal H_{ki}(t), w_i(t), b_{ki}, \alpha_{kl}(t) \} \biggr], \end{eqnarray}\] and with regard to the estimation we postulate the following set of hierarchical priors: \[\begin{eqnarray} \lambda &\sim& \mathcal N(0, \tau_{\lambda} M)\\&&\\ \tau_{\lambda} &\sim& \mbox{inv-Gamma}(1, 0.005), \end{eqnarray}\] where \(K = \Delta_r^\top \Delta_r + 10^{-6} \mathbf I\), and \(\Delta_r\) denotes the \(r\)-th order differences matrix as in the specification of the baseline hazard presented previously.

The time-varying coefficients are specified via the Interactions argument of mvJointModelBayes(). Just as an illustration we extend the univariate model JMFit1 to allow for a time-varying association parameter:

Ints_tveffect <-  list("log(serBilir)_value" = ~ 0 + tve(Time, df = 8))
JMFit1_tveffect <- mvJointModelBayes(MixedModelFit1, CoxFit, timeVar = "year",
                                    Interactions = Ints_tveffect)

The plot() method can be used to depict the time-varying effect

plot(JMFit1_tveffect, which = "tv_effect")

The following function can be used to include in the same plot the time-constant and time-varying effects

plot_tveffects <- function (objTVeff, objCeff, cut_time = NULL, 
                            scale = c("logHR", "HR"), include = c("const", "tv", "both")) {
    scale <- match.arg(scale)
    include <- match.arg(include)
    tv_effect <- plot(objTVeff, "tv_effect", return_values = TRUE)[[1]]
    const_effect <- cbind(times = tv_effect[, 1],
                          est = objCeff$statistics$postMeans$alphas,
                          low = objCeff$statistics$CIs$alphas[1, ],
                          upp = objCeff$statistics$CIs$alphas[2, ])
    if (!is.null(cut_time)) {
        tv_effect <- tv_effect[tv_effect[, "times"] <= cut_time, ]
        const_effect <- const_effect[const_effect[, "times"] <= cut_time, ]
    }
    if (scale == "HR") {
        tv_effect[, -1] <- exp(tv_effect[, -1])
        const_effect[, -1] <- exp(const_effect[, -1])
    }
    matplot(tv_effect[, 1], cbind(tv_effect[, -1], const_effect[, -1]), type = "n",
            xlab = "Time (years)", 
            ylab = if (scale == "HR") "Hazard Ratio" else "log Hazard Ratio")
    if (include %in% c("const", "both")) {
        polygon(c(const_effect[, 1], rev(const_effect[, 1])), 
                c(const_effect[, 3], rev(const_effect[, 4])), 
                col = "lightgrey", border = NA)
    }
    if (include %in% c("tv", "both")) {
        polygon(c(tv_effect[, 1], rev(tv_effect[, 1])), 
                c(tv_effect[, 3], rev(tv_effect[, 4])), 
                col = "grey", border = NA)
    }
    if (include %in% c("const", "both")) {
        matlines(const_effect[, 1], const_effect[, -1], lty = c(1, 2, 2), 
                 col = 1, lwd = c(2, 1, 1))
    }
    if (include %in% c("tv", "both")) {
        matlines(tv_effect[, 1], tv_effect[, -1], lty = c(1, 2, 2), 
                 col = 2, lwd = c(2, 1, 1))
    }
}

plot_tveffects(JMFit1_tveffect, JMFit1, include = "both", scale = "HR", cut_time = 12)

Left, Right and Interval Censored Data

Function mvJointModelBayes() can also simultaneously handled, right, left and interval censored data. As an example, we first create artificial right, left and interval censoring in the PBC data set:

pbc2$status2 <- as.numeric(pbc2$status != "alive")
pbc2.id$status2 <- as.numeric(pbc2.id$status != "alive")
pbc2$status3 <- as.character(pbc2$status)
ff <- function (x) {
    out <- if (x[1L] %in% c('dead', 'transplanted')) 'dead' else 
        switch(sample(1:3, 1), '1' = "alive", '2' = "left", '3' = "interval")
    rep(out, length(x))
}
pbc2$status3 <- unlist(with(pbc2, lapply(split(status3, id), ff)), use.names = FALSE)
pbc2$status3 <- unname(with(pbc2, sapply(status3, function (x) 
    switch(x, 'dead' = 1, 'alive' = 0, 'left' = 2, 'interval' = 3))))
pbc2$yearsU <- as.numeric(NA)
pbc2$yearsU[pbc2$status3 == 3] <- pbc2$years[pbc2$status3 == 3] + 
    runif(sum(pbc2$status3 == 3), 0, 4)
pbc2.id <- pbc2[!duplicated(pbc2$id), ]

Following, we fit a Weibull accelerated failure time (AFT) model accounting for the different types of censoring:

survFit <- survreg(Surv(years, yearsU, status3, type = "interval") ~ drug + age, 
                   data = pbc2.id, model = TRUE)

Finally, we fit the joint model with the same call to mvJointModelBayes() as we have previously used.

JMFit1_intcens <- mvJointModelBayes(MixedModelFit1, survFit, timeVar = "year")
summary(JMFit1_intcens)

## 
## Call:
## mvJointModelBayes(mvglmerObject = MixedModelFit1, survObject = survFit, 
##     timeVar = "year")
## 
## Data Descriptives:
## Number of Groups: 312        Number of events: 169 (54.2%)
## Number of Observations:
##   log(serBilir): 1945
## 
## Random-effects covariance matrix:
##       StdDev    Corr
## (I)1  1.0081  (Int)1
## yer1  0.1735  0.4016
## 
## Survival Outcome:
##                     PostMean  StDev  StErr    2.5%  97.5%     P
## drugD-penicil        -0.0508 0.1221 0.0026 -0.2823 0.2002 0.702
## age                   0.0228 0.0063 0.0002  0.0109 0.0350 0.000
## log(serBilir)_value   0.7072 0.0583 0.0018  0.5957 0.8221 0.000
## 
## Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
##             PostMean  StDev  StErr   2.5%  97.5% P
## (Intercept)   0.4956 0.0596 0.0019 0.3767 0.6096 0
## year          0.1764 0.0130 0.0004 0.1522 0.2026 0
## sigma         0.3489 0.0067 0.0002 0.3366 0.3627 0
## 
## MCMC summary:
## iterations: 900 
## burn-in: 600 
## thinning: 300 
## time: 5.3 min

As previously formulated in the Theory Section, even if the survFit provided to mvJointModelBayes() is an AFT model, the function internally specifies a relative risk model.

Multivariate Joint Models

Dimitris Rizopoulos

2018-04-17