Data Engineering and BI courses are free this week!
jointModel(lmeObject, survObject, timeVar,
parameterization = c("value", "slope", "both"),
method = c("weibull-PH-GH", "weibull-PH-aGH", "weibull-AFT-GH",
"weibull-AFT-aGH", "piecewise-PH-GH", "piecewise-PH-aGH",
"Cox-PH-GH", "Cox-PH-aGH", "spline-PH-GH", "spline-PH-aGH",
"ch-Laplace"),
interFact = NULL, derivForm = NULL, lag = 0, scaleWB = NULL,
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
value a formula for the interaction terms corresponding to the
value parameterization, slope a formula for the interaction terms corresponding to the
slope parameterizatifixed a formula representing the derivative of the fixed-effects part of the
liner mixed model with respect to time, indFixed a numeric vector indicating which fixed effects of lmeObject method = "weibull-AFT-GH" or method = "weibull-PH-GH". The default NULL means that the scale
parametercontrol argument.jointModelObject for the components of the fit.jointModel 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 let $w_i$ denote 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$),
$\alpha$ the association parameter for $m_i(t)$, $m_i'(t)$ the derivative of $m_i(t)$ with respect to $t$, and
$\alpha_d$ the association parameter for $m_i'(t)$. 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 these options the linear predictor for the
survival submodel is written as parameterization = "value", parameterization = "slope", and parameterization = "both", where in all the above the value
of $k$ is specified by the lag argument and $m_i'(t) = d m_i(t) / dt$. If interFact is specified, then
$m_i{min(t-k, 0)}$ and/or $m_i'{min(t-k, 0)}$ are multiplied with the design matrices derived from the formulas
supplied as the first two arguments of interFact, respectively. In this case $\alpha$ and/or $\alpha_s$ become vectors of
association parameters.
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 standard Gauss-Hermite quadrature rule when the chosen option for the method
argument contains the string "GH", and the (pseudo) adaptive Gauss-Hermite rule when the chosen option for the method
argument contains the string "aGH". For method = "ch-Laplace" the fully exponential Laplace approximation described in
Rizopoulos et al. (2009) is used. The (pseudo) adaptive Gauss-Hermite and the Laplace approximation are particularly useful 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,
rocJM,
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)
# we also include an interaction term of log(serBilir) with drug
fitJOINT <- jointModel(fitLME, fitSURV, timeVar = "year",
interFact = list(value = ~ drug, data = pbc2.id))
fitJOINT
summary(fitJOINT)
# a joint model in which the risk for and event depends both on the true value of
# marker and the true value of the slope of the longitudinal trajectory
lmeFit <- lme(sqrt(CD4) ~ obstime * drug, random = ~ obstime | patient, data = aids)
coxFit <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)
# to fit this model we need to specify the 'derivForm' argument, which is a list
# with first component the derivative of the fixed-effects formula of 'lmeFit' with
# respect to 'obstime', second component the indicator of which fixed-effects
# coefficients correspond to the previous defined formula, third component the
# derivative of the random-effects formula of 'lmeFit' with respect to 'obstime',
# and fourth component the indicator of which random-effects correspond to the
# previous defined formula
dForm <- list(fixed = ~ 1 + drug, indFixed = c(2, 4), random = ~ 1, indRandom = 2)
jointModel(lmeFit, coxFit, timeVar = timeVar, method = "spline-PH-GH",
parameterization = "both", derivForm = dForm)
# 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 with a spline approximate baseline hazard
fitJOINT <- jointModel(fitLME, fitCOX,
timeVar = "obstime", method = "spline-PH-GH")
fitJOINT
summary(fitJOINT)Run the code above in your browser using DataLab