jointModel: Joint Models for Longitudinal and Survival Data

Description

This function fits shared parameter models for the joint modelling of normal longitudinal responses and time-to-event data under a maximum likelihood approach. Various options for the survival model are available.

Usage

jointModel(lmeObject, survObject, timeVar, 
    method = c("weibull-AFT-GH", "weibull-PH-GH", 
    "piecewise-PH-GH", "Cox-PH-GH", "spline-PH-GH", "ch-Laplace"), 
    init = NULL, control = list(), ...)

Arguments

lmeObject

an object inheriting from class lme (see also Note).

survObject

an object inheriting from class coxph or class survreg. In the call to coxph() or survreg(), you need to specify the argument x = TRUE such that the design matrix is contained in

timeVar

a character string indicating the time variable in the linear mixed effects model.

method

a character string specifying the type of joint model to fit. See Details.

init

a list of user-specified initial values. The initial values list must have the following components: [object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object] If the user-specified list of in

control

a list of control values with components: [object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],

...

options passed to the control argument.

Value

See jointModelObject for the components of the fit.

Details

The jointModel function fits joint models for longitudinal and survival data. For the longitudinal responses the linear mixed effects model represented by the lmeObject is assumed. For the survival times three options are available. In particular, let $t_i$ denote the time-to-event for the $i$th sample unit, $x_i$ denote the vector of baseline covariates in survObject, with associated parameter vector $\gamma$, $W_i(t)$ the value of the longitudinal outcome at time point $t$ (i.e., $W_i(t)$ equals the fixed effects part + random effects part of the linear mixed effects model for sample unit $i$), and $\alpha$ the association parameter. Then, for method = "weibull-AFT-GH" a time-dependent Weibull model under the accelerated failure time formulation is assumed. For method = "weibull-PH-GH" a time-dependent relative risk model is postulated with a Weibull baseline risk function. For method = "piecewise-PH-GH" a time-dependent relative risk model is postulated with a piecewise constant baseline risk function. For method = "spline-PH-GH" a time-dependent relative risk model is assumed in which the log baseline risk function is approximated using B-splines. For method = "ch-Laplace" an additive model on the log cumulative hazard scale is assumed (see Rizopoulos et al., 2009 for more info). Finally, for method = "Cox-PH-GH" a time-dependent relative risk model is assumed where the baseline risk function is left unspecified (Wulfsohn and Tsiatis, 1997). For all survival models except for the time-dependent proportional hazards model, the optimization algorithm starts with EM iterations, and if convergence is not achieved, it switches to quasi-Newton iterations (i.e., BFGS in optim() or nlminb(), depending on the value of the optimizer control argument). For the time-dependent proportional hazards model only the EM algorithm is used. During the EM iterations, convergence is declared if either of the following two conditions is satisfied: (i) $L(\theta^{it}) - L(\theta^{it - 1}) < tol_3 { | L(\theta^{it - 1}) | + tol_3 }$, or (ii) $\max { | \theta^{it} - \theta^{it - 1} | / ( | \theta^{it - 1} | + tol_1) } < tol_2$, where $\theta^{it}$ and $\theta^{it - 1}$ is the vector of parameter values at the current and previous iterations, respectively, and $L(.)$ is the log-likelihood function. The values for $tol_1$, $tol_2$ and $tol_3$ are specified via the control argument. During the quasi-Newton iterations, the default convergence criteria of either optim() or nlminb() are used. The required integrals are approximated using the Gauss-Hermite quadrature rule for method = "weibull-AFT-GH", method = "weibull-PH-GH", method = "piecewise-PH-GH", method = "spline-PH-GH" and method = "Cox-PH-GH", whereas for method = "ch-Laplace" the fully exponential Laplace approximation described in Rizopoulos et al. (2009) is used. This last option is more suitable when high-dimensional random effects vectors are considered (e.g., when modelling nonlinear subject-specific trajectories with splines or high-order polynomials).

References

Henderson, R., Diggle, P. and Dobson, A. (2000) Joint modelling of longitudinal measurements and event time data. Biostatistics 1, 465--480. Hsieh, F., Tseng, Y.-K. and Wang, J.-L. (2006) Joint modeling of survival and longitudinal data: Likelihood approach revisited. Biometrics 62, 1037--1043. Rizopoulos, D., Verbeke, G. and Lesaffre, E. (2009) Fully exponential Laplace approximation for the joint modelling of survival and longitudinal data. Journal of the Royal Statistical Society, Series B 71, 637--654. Rizopoulos, D., Verbeke, G. and Molenberghs, G. (2009) Multiple-imputation-based residuals and diagnostic plots for joint models of longitudinal and survival outcomes. Biometrics, to appear (doi: 10.1111/j.1541-0420.2009.01273.x). Tsiatis, A. and Davidian, M. (2004) Joint modeling of longitudinal and time-to-event data: an overview. Statistica Sinica 14, 809--834. Wulfsohn, M. and Tsiatis, A. (1997) A joint model for survival and longitudinal data measured with error. Biometrics 53, 330--339.

Examples

Run this code

# linear mixed model fit (random intercepts)
fitLME <- lme(log(serBilir) ~ drug * year, random = ~ 1 | id, data = pbc2)
# survival regression fit
fitSURV <- survreg(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)
# joint model fit, under the (default) Weibull model
fitJOINT <- jointModel(fitLME, fitSURV, timeVar = "year")
fitJOINT
summary(fitJOINT)

# linear mixed model fit (random intercepts + random slopes)
fitLME <- lme(log(serBilir) ~ drug * year, random = ~ year | id, data = pbc2)
# survival regression fit
fitSURV <- survreg(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)
# joint model fit, under the (default) Weibull model
fitJOINT <- jointModel(fitLME, fitSURV, timeVar = "year")
fitJOINT
summary(fitJOINT)

# linear mixed model fit
fitLME <- lme(sqrt(CD4) ~ obstime * drug - drug, 
    random = ~ 1 | patient, data = aids)
# cox model fit
fitCOX <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)
# joint model fit
fitJOINT <- jointModel(fitLME, fitCOX, 
    timeVar = "obstime", method = "spline-PH-GH")
fitJOINT
summary(fitJOINT)

Run the code above in your browser using DataLab