jointModel(lmeObject, survObject, timeVar,
method = c("weibull-AFT-GH", "weibull-PH-GH", "piecewise-PH-GH", "ch-GH", "ph-GH", "ch-Laplace"),
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 = "ph-GH" for which only the EM algorithm is available.method = "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" or when method = "ch-GH" or method = "ch-Laplace".
For method = "piecewise-PH-GH" default is to place equally-spaced lng.in.kn knots (see below) in the
quantiles of the observed event times. For method = "ch-GH" and method = "ch-Laplace" default is
to place lng.in.kn knots (see below) at the quantiles of the uncensored log survival times.method = "piecewise-PH-GH" where it
denotes the number of internal knots for the piecewise constant baseline risk function or when method = "ch-GH"
or method = "ch-Laplace" where it denotes the number of internal knots for B-splines approximation of the log cumulative
baseline hazard. Default is 6 when method = "piecewise-PH-GH" and 3 otherwise.method = "ch-GH" or method = "ch-Laplace".method = "weibull-AFT-GH}, code{method = "weibull-PH-GH and method = "ch-GH" default is 21 for one- or
two-dimensional integration and 11 otherwise. For method = "ph-GH" default is 21 for one-dimensional
integration, 11 for two-dimensional integration and 9 otherwise.TRUE the parameter estimates and the log-likelihood value are printed during
the optimization procedure. Default is FALSE.controljointModel 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 = "ch-GH" and 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 = "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 = "ch-GH" and method = "ph-GH",
whereas for method = "ch-Laplace" the fully exponential Laplace approximation described in Rizopoulos et al. (2008) 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# 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, under the additive log cumulative hazard model
fitJOINT <- jointModel(fitLME, fitCOX,
timeVar = "obstime", method = "ch-GH")
fitJOINT
summary(fitJOINT)Run the code above in your browser using DataLab