jointModelBayes output

Question

Consider the following example from a workshop by Dimitris Rizopoulos workshop.

We have the following joint model:

\begin{equation} \begin{aligned} y_i(t) &= m_i(t) + \epsilon_i(t) \\ &= \beta_0 + \beta_1t + \beta_2(t \times \text{ddI}_i) + b_{i0} + b_{i1}t + \epsilon_i(t), \hspace{1cm} \epsilon_i(t) \sim N(0, \sigma^2) \\ \\ h_i(t)& = h_0(t)\exp(\gamma \text{ddI}_i + \alpha m_i(t)). \end{aligned} \end{equation}

Here $b = (b_0, b_1)' \sim N(0, D)$, and

\begin{equation} log(h_0(t)) = \gamma_{h_0,0} + \sum_{q=1}^Q\gamma_{h_0,q}B_q(t,v). \end{equation}

Using code from slide 96:

library(JM)
library(JMbayes)
#Fitting linear mixed effects model
lmeFit <- lme(CD4 ~ obstime + obstime:drug,
              random = ~ obstime | patient, data = aids)

#Fitting Cox PH model
coxFit <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)

#Fitting joint model with Bayesian approach
jointFitBayes <- jointModelBayes(lmeFit, coxFit, timeVar = "obstime")
summary(jointFitBayes)

The summary displays:

Call:
jointModelBayes(lmeObject = lmeFit, survObject = coxFit, timeVar = "obstime")

Data Descriptives:
Longitudinal Process        Event Process
Number of Observations: 1405    Number of Events: 188 (40.3%)
Number of subjects: 467

Joint Model Summary:
Longitudinal Process: Linear mixed-effects model
Event Process: Relative risk model with penalized-spline-approximated 
        baseline risk function
Parameterization: Time-dependent value 

 LPML      DIC       pD
 -Inf 12743.78 1708.641

Variance Components:
             StdDev    Corr
(Intercept)  4.5930  (Intr)
obstime      0.5745 -0.0268
Residual     1.7174        

Coefficients:
Longitudinal Process
                  Value Std.Err Std.Dev    2.5%   97.5%      P
(Intercept)      7.1967  0.0057  0.2217  6.7563  7.6278 <0.001
obstime         -0.2301  0.0011  0.0427 -0.3134 -0.1449 <0.001
obstime:drugddI  0.0072  0.0015  0.0600 -0.1109  0.1245  0.897

Event Process
           Value Std.Err  Std.Dev    2.5%    97.5%      P
drugddI   0.3417  0.0081   0.1908 -0.0550   0.7022   0.08
Assoct   -0.2966  0.0021   0.0438 -0.3858  -0.2143 <0.001
tauBs   186.5855 14.9435 146.5236 21.7815 565.1833     NA

MCMC summary:
iterations: 20000 
adapt: 3000 
burn-in: 3000 
thinning: 10 
time: 1.9 min

Now, if you look into the jointFitBayes object and look at the mcmc list, we find 8 samples (corresponding to 8 parameter vectors, where a vector could be of length 1).

The list contains samples for (among others) D and tauBs. The posterior mean for tauBs is 186.5855, and the posterior mean for D is

\begin{bmatrix} 21.0959 & -0.0708 \\ -0.0708 & 0.3300 \end{bmatrix}.

In the help section it states:

?jointModelBayes
tau
 the precision parameter from the linear mixed effects model (i.e., τ = 1/σ^2 with σ denoting the error terms standard deviation).
 invD
 the inverse variance-covariance matrix of the random effects.

But, there is no description about tauBs or D. If D is the variance-covariance matrix of the random effects, then why does it not match the variance component part of the summary?:

\begin{bmatrix} 4.5930 & -0.0268 \\ -0.0268 & 0.5745 \end{bmatrix}.

In addition, from the summary, the posterior mean for $\sigma$ is 1.7174. This would mean that the posterior mean for $\tau$ should be 0.3390448. But, the posterior mean for tauBs is 186.5855.

What can tauBs be?

Example extended: we can use the joint models to predict survival probabilities in the future. This can be done as follows (see slide 157):

sfit <- survfitJM(jointFitBayes, newdata = aids[aids$patient == "2", ], idVar = "patient")

However, if I do the following, cf slides 127 and 157 in the workshop, I receive an error:

#Fitting linear mixed effects model
lmeFit <- lme(CD4 ~ obstime + obstime:drug,
              random = ~ obstime | patient, data = aids)

#Fitting Cox PH model
coxFit <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)

#Fitting joint model using MLE
jointFit <- jointModel(lmeFit, coxFit, timeVar = "obstime",
                       method = "piecewise-PH-aGH")
sfit <- survfitJM(jointFit, newdata = aids[aids$patient == "2", ], idVar = "patient")

Error in UseMethod("survfitJM") : 
  no applicable method for 'survfitJM' applied to an object of class "jointModel"

Finally, using

sfit <- survfitJM(jointFitBayes, newdata = aids[aids$patient == "3", ], idVar = "patient")

plot(sfit, include.y = TRUE)

I obtain the following plot

However, in the workshop on slide 157 (with corresponding plots on slides 156), the same code plots the median survival probability with corresponding $95\%$ prediction intervals.

Why is this?

Answer 1

There are several questions here. Some answers:

• tauBs is the smoothing parameter for the penalized B-spline approximation of the log baseline hazard function.
• D is indeed the variance-covariance matrix of the random effects. Note that as the name suggests, this includes variances and covariances. In the output, you obtain standard deviations and correlations.
• You should not have both packages JM and JMbayes loaded.
• To include a confidence interval for the dynamic predictions from a joint model you need to use the arguments include.y = TRUE and conf.int = TRUE in the call to the plot() method.