Description Usage Arguments Details Value Note Author(s) References See Also Examples
This function fits shared parameter models for the joint modelling of normal longitudinal responses and timetoevent data under a maximum likelihood approach. Various options for the survival model are available.
1 2 3 4 5 6 7 8  jointModel(lmeObject, survObject, timeVar,
parameterization = c("value", "slope", "both"),
method = c("weibullPHaGH", "weibullPHGH", "weibullAFTaGH",
"weibullAFTGH", "piecewisePHaGH", "piecewisePHGH",
"CoxPHaGH", "CoxPHGH", "splinePHaGH", "splinePHGH",
"chLaplace"),
interFact = NULL, derivForm = NULL, lag = 0, scaleWB = NULL,
CompRisk = FALSE, init = NULL, control = list(), ...)

lmeObject 
an object inheriting from class 
survObject 
an object inheriting from class 
timeVar 
a character string indicating the time variable in the linear mixed effects model. 
parameterization 
a character string indicating the type of parameterization. See Details 
method 
a character string specifying the type of joint model to fit. See Details. 
interFact 
a list with components 
derivForm 
a list with components 
lag 
a numeric scalar denoting a lag effect in the timedependent covariate represented by the mixed model; default is 0. 
scaleWB 
a numeric scalar denoting a fixed value for the scale parameter of the Weibull hazard; used only when

CompRisk 
logical; should a competing risks joint model be fitted. 
init 
a named list of userspecified initial values:
When this list of initial values does not contain some of these components or contains components not of the appropriate length, then the default initial values are used instead. 
control 
a list of control values with components:

... 
options passed to the 
Function 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
γ, 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 fixedeffects part +
randomeffects part of the linear mixed effects model for sample unit i),
α the association parameter for m_i(t), m_i'(t) the derivative of m_i(t) with respect to t, and
α_d the association parameter for m_i'(t). Then, for method = "weibullAFTGH"
a timedependent Weibull model under
the accelerated failure time formulation is assumed. For method = "weibullPHGH"
a timedependent relative risk model is postulated
with a Weibull baseline risk function. For method = "piecewisePHGH"
a timedependent relative risk model is postulated with a
piecewise constant baseline risk function. For method = "splinePHGH"
a timedependent relative risk model is assumed in which the
log baseline risk function is approximated using Bsplines. For method = "chLaplace"
an additive model on the log cumulative hazard
scale is assumed (see Rizopoulos et al., 2009 for more info). Finally, for method = "CoxPHGH"
a timedependent 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
eta = gamma' w_i + alpha m_i{max(tk,0)},
when
parameterization = "value"
,
eta = gamma' w_i + alpha m_i'{max(tk, 0)}
when parameterization = "slope"
, and
eta = gamma' w_i + alpha m_i{max(tk, 0)} + α_s m_i'{max(tk, 0)},
when 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\{max(tk, 0)\} and/or m_i'\{max(tk, 0)\} are multiplied with the design matrices derived from the formulas
supplied as the first two arguments of interFact
, respectively. In this case α and/or α_s become vectors of
association parameters.
For method = "splinePHGH"
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 timedependent proportional hazards model, the optimization algorithm starts
with EM iterations, and if convergence is not achieved, it switches to quasiNewton iterations (i.e., BFGS in
optim()
or nlminb()
, depending on the value of the optimizer
control argument). For method = "CoxPHGH"
only the
EM algorithm is used. During the EM iterations, convergence is declared if either of the following two conditions is satisfied: (i)
L(θ^{it})  L(θ^{it  1}) < tol_3 \{  L(θ^{it  1})  + tol_3 \} , or (ii)
\max \{  θ^{it}  θ^{it  1}  / (  θ^{it  1}  + tol_1) \} < tol_2, where θ^{it} and
θ^{it  1} is the vector of parameter values at the current and previous iterations, respectively, and L(.) is the
loglikelihood function. The values for tol_1, tol_2 and tol_3 are specified via the control
argument. During the
quasiNewton iterations, the default convergence criteria of either optim()
or nlminb()
are used.
The required integrals are approximated using the standard GaussHermite quadrature rule when the chosen option for the method
argument contains the string "GH", and the (pseudo) adaptive GaussHermite rule when the chosen option for the method
argument contains the string "aGH". For method = "chLaplace"
the fully exponential Laplace approximation described in
Rizopoulos et al. (2009) is used. The (pseudo) adaptive GaussHermite and the Laplace approximation are particularly useful when
highdimensional random effects vectors are considered (e.g., when modelling nonlinear subjectspecific trajectories with splines
or highorder polynomials).
See jointModelObject
for the components of the fit.
1. The lmeObject
argument should represent a linear mixed model object with a simple randomeffects
structure, i.e., only the pdDiag()
class is currently allowed.
2. The lmeObject
object should not contain any withingroup correlation structure (i.e., correlation
argument of lme()
) or withingroup heteroscedasticity structure (i.e., weights
argument of lme()
).
3. It is assumed that the linear mixed effects model lmeObject
and the survival model survObject
have been
fitted to the same subjects. Moreover, it is assumed that the ordering of the subjects is the same for both
lmeObject
and survObject
, i.e., that the first line in the data frame containing the event times
corresponds to the first set of lines identified by the grouping variable in the data frame containing the repeated
measurements, and so on.
4. In the print
and summary
generic functions for class jointModel
, the estimated coefficients (and
standard errors for the summary
generic) for the event process are augmented with the element "Assoct" that
corresponds to the association parameter α and the element "Assoct.s" that corresponds to the parameter
α_s when parameterization
is "slope"
or "both"
(see Details).
5. The standard errors returned by the summary
generic function for class jointModel
when
method = "CoxPHGH"
are based on the profile score vector (i.e., given the NPMLE for the unspecified baseline
hazard). Hsieh et al. (2006) have noted that these standard errors are underestimated.
6. As it is the case for all types of mixed models that require numerical integration, it is advisable (especially in difficult datasets) to check the stability of the maximum likelihood estimates with an increasing number of GaussHermite quadrature points.
7. It is assumed that the scale of the time variable (e.g., days, months years) is the same in both lmeObject
and survObject
.
Dimitris Rizopoulos [email protected]
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. (2012a) Joint Models for Longitudinal and TimetoEvent Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.
Rizopoulos, D. (2012b) Fast fitting of joint models for longitudinal and event time data using a pseudoadaptive Gaussian quadrature rule. Computational Statistics and Data Analysis 56, 491–501.
Rizopoulos, D. (2011) Dynamic predictions and prospective accuracy in joint models for longitudinal and timetoevent data. Biometrics 67, 819–829.
Rizopoulos, D. (2010) JM: An R package for the joint modelling of longitudinal and timetoevent data. Journal of Statistical Software 35 (9), 1–33. http://www.jstatsoft.org/v35/i09/
Rizopoulos, D., Verbeke, G. and Lesaffre, E. (2009) Fully exponential Laplace approximations 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. (2010) Multipleimputationbased residuals and diagnostic plots for joint models of longitudinal and survival outcomes. Biometrics 66, 20–29.
Tsiatis, A. and Davidian, M. (2004) Joint modeling of longitudinal and timetoevent 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.
jointModelObject
,
anova.jointModel
,
coef.jointModel
,
fixef.jointModel
,
ranef.jointModel
,
fitted.jointModel
,
residuals.jointModel
,
plot.jointModel
,
survfitJM
,
rocJM
,
dynCJM
,
aucJM
,
prederrJM
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83  ## Not run:
# 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 fixedeffects formula of 'lmeFit' with
# respect to 'obstime', second component the indicator of which fixedeffects
# coefficients correspond to the previous defined formula, third component the
# derivative of the randomeffects formula of 'lmeFit' with respect to 'obstime',
# and fourth component the indicator of which randomeffects correspond to the
# previous defined formula
dForm < list(fixed = ~ 1 + drug, indFixed = c(2, 4), random = ~ 1, indRandom = 2)
jointModel(lmeFit, coxFit, timeVar = "obstime", method = "splinePHaGH",
parameterization = "both", derivForm = dForm)
# Competing Risks joint model
# we first expand the PBC dataset in the competing risks long format
# with two competing risks being death and transplantation
pbc2.idCR < crLong(pbc2.id, "status", "alive")
# we fit the linear mixed model as before
lmeFit.pbc < lme(log(serBilir) ~ drug * ns(year, 3),
random = list(id = pdDiag(form = ~ ns(year, 3))), data = pbc2)
# however, for the survival model we need to use the data in the long
# format, and include the competing risks indicator as a stratification
# factor. We also take interactions of the baseline covariates with the
# stratification factor in order to allow the effect of these covariates
# to be different for each competing risk
coxCRFit.pbc < coxph(Surv(years, status2) ~ (drug + sex)*strata + strata(strata),
data = pbc2.idCR, x = TRUE)
# the corresponding joint model is fitted simply by including the above
# two submodels as main arguments, setting argument CompRisk to TRUE,
# and choosing as method = "splinePHaGH". Similarly as above, we also
# include strata as an interaction factor to allow serum bilirubin to
# have a different effect for each of the two competing risks
jmCRFit.pbc < jointModel(lmeFit.pbc, coxCRFit.pbc, timeVar = "year",
method = "splinePHaGH",
interFact = list(value = ~ strata, data = pbc2.idCR),
CompRisk = TRUE)
summary(jmCRFit.pbc)
# 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 splineapproximated baseline hazard function
fitJOINT < jointModel(fitLME, fitCOX,
timeVar = "obstime", method = "splinePHaGH")
fitJOINT
summary(fitJOINT)
## End(Not run)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.