jointModel(lmeObject, survObject, timeVar,
method = c("weibull-AFT-GH", "weibull-PH-GH",
"piecewise-PH-GH", "Cox-PH-GH", "spline-PH-GH", "ch-Laplace"),
lag = 0, init = NULL, control = list(), ...)lme (see also Note).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
jointModelObject for the components of the fit.TRUE only the EM algorithm is used in the optimization, otherwise if
convergence has not been achieved a quasi-Newton algorithm is initiated. Default is FALSE except for
method = "Cox-PH-GH" for which only the EM algorithm is available.method = "Cox-PH-GH" for which
the default is 200.sqrt(.Machine$double.eps).numeriDeriv = "cd" a larger value (e.g., 1e-04) is suggested.parscale control argument for optim(), or the scale argument for
nlminb(). It should be a numeric vector of length equal to the number of parameters. Default is 0.01
for all parameters.method = "ch-Laplace". Default is 0.1.method = "ch-Laplace".method = "piecewise-PH-GH", method = "spline-PH-GH" or method = "ch-Laplace".
The default is to place equally-spaced lng.in.kn knots in the quantiles of the observed event times. For stratified models
fitted with method = "spline-PH-GH" this should be a list with elements numeric vectors of knots positions for each strata.method = "piecewise-PH-GH" where it
denotes the number of internal knots for the piecewise constant baseline risk function or when method = "spline-PH-GH"
or method = "ch-Laplace" where it denotes the number of internal knots for B-splines approximation of the log
baseline hazard. Default is 6 when method = "piecewise-PH-GH" and 5 otherwise.TRUE (the default), then the same knots are used in the approximation of the
baseline risk function in different strata when method = "spline-PH-GH".method = "spline-PH-GH" or method = "ch-Laplace".method = "weibull-PH-GH" and method = "weibull-AFT-GH" 15 are used,
whereas for method = "piecewise-PH-GH" 7.TRUE, the parameter estimates and the log-likelihood value are printed during
the optimization procedure. Default is FALSE.controljointModel fits joint models for longitudinal and survival data (more detailed information about the formulation of these
models can be found in Rizopoulos (2010)). 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,
$w_i$ denotes the vector of baseline covariates in survObject, with associated parameter vector $\gamma$, $m_i(t)$ the value
of the longitudinal outcome at time point $t$ as approximated by the linear mixed model (i.e., $m_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 method = "spline-PH-GH" it is also allowed to include stratification factors. These should be included in the specification of
the survObject using function strata(). Note that in this case survObject must only be a 'coxph' object.
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 method = "Cox-PH-GH" 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).jointModelObject,
anova.jointModel,
coef.jointModel,
fixef.jointModel,
ranef.jointModel,
fitted.jointModel,
residuals.jointModel,
plot.jointModel,
survfitJM,
dynC# 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