0
$\begingroup$

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 plotforecasts

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?

$\endgroup$

1 Answer 1

2
$\begingroup$

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.
$\endgroup$
3
  • $\begingroup$ Using plot(sfit, include.y = TRUE, conf.int = TRUE) works. But, both the mean and median are still shown with the $95\%$ PI. Is there an option to remove the mean? Is there an option to remove the median? $\endgroup$
    – JLee
    Commented Sep 18, 2019 at 12:25
  • 1
    $\begingroup$ You can use estimator = “mean” or “median” $\endgroup$ Commented Sep 18, 2019 at 12:28
  • $\begingroup$ One thing I forgot to ask: You mention different parameterizations which specify the type of association structure between the longitudinal and survival processes. The option "param" in jointModelBayes allows you to change the parameterization. However, the cumulative effects parametrization is not an option. Is there an option I am missing that allows this? $\endgroup$
    – JLee
    Commented Sep 18, 2019 at 18:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.