Description Usage Arguments Details Value Competing risks Separate models Note Author(s) References See Also Examples
This generic function fits a joint model with random latent association, building on the formulation described in Wulfsohn and Tsiatis (1997) while allowing for the presence of longitudinal and survival covariates, and three choices for the latent process. The link between the longitudinal and survival processes can be proportional or separate. When failure is attributable to 2 separate causes, a competing risks joint model is fitted as per Williamson et al. (2008).
1 2 3 4 5 6 7 8 9 10 11 12 13 14 
data 
an object of class 
long.formula 
a formula object with the response variable, and the covariates to include in the longitudinal submodel. 
surv.formula 
a formula object with the survival time, censoring
indicator and the covariates to include in the survival submodel. The
response must be a survival object as returned by the

model 
a character string specifying the type of latent association.
This defaults to the intercept and slope version as seen in Wulfsohn and
Tsiatis (1997). For association via the random intercept only, choose

sepassoc 
logical value: if 
longsep 
logical value: if 
survsep 
if 
gpt 
the number of quadrature points across which the integration with
respect to the random effects will be performed. Defaults to 
lgpt 
the number of quadrature points which the loglikelihood is
evaluated over following a model fit. This defaults to 
max.it 
the maximum number of iterations of the EM algorithm that the
function will perform. Defaults to 
tol 
the tolerance level before convergence of the algorithm is deemed
to have occurred. Default value is 
verbose 
if 
The joint
function fits a joint model to survival and
longitudinal data. The formulation is similar to Wulfsohn and Tsiatis
(1997). A linear mixed effects model is assumed for the longitudinal data,
namely
Y_i = X_{i1}(t_i)^Tβ_1 + D_i(t_i)^T U_i + ε_i,
where U_i is a vector of random effects, (U_{0i}, … U_{qi}) whose length depends on the model chosen, i.e. q = 1 for the random intercept model. D_i is the random effects covariate matrix, which will be timedependent for all but the random intercept model. X_{i1} is the longitudinal design matrix for unit i, and t_i is the vector of measurement times for subject i. Measurement error is represented by ε_i.
The Cox proportional hazards model is adopted for the survival data, namely
λ(t) = λ_0(t) \exp\{{X_{i2}(t)^T β_2 + D_i(t)(γ^T U_i)}\}.
The parameter γ determines the level of association between the two processes. For the intercept and slope model with separate association we have
D_i(t) (γ^T U_i) = γ_0 U_{0i} + γ_1 U_{1i} t,
whereas under proportional association
D_i(t)(γ^T U_i) = γ (U_{0i} + U_{1i} t).
X_{i2} is the vector of survival covariates for unit i. The baseline hazard function is λ_0(t).
The function uses an EM algorithm to estimate parameters in the joint
model. Starting values are provided by calls to standard R functions
lme
and coxph
for the
longitudinal and survival components, respectively.
A list containing the parameter estimates from the joint model and, if required, from either or both of the separate analyses. The combined loglikelihood from a separate analysis and the loglikelihood from the joint model are also produced as part of the fit.
If failure can be attributed to 2 causes, i.e. socalled competing risks events data, then a causespecific hazards model is adopted, namely
λ_g(t) = λ_{0g}(t) \exp\{{X_{i2}(t)^T β_2^{(g)} + D_i(t)(γ^T U_i)}\},
where g = 1,2 denotes the failure type, β_2^{(g)} (g =
1,2) are causespecific hazard parameters corresponding to the same
covariates, and λ_{0g}(t) are causespecific baseline hazard
functions. For this data, a proportional association structure is assumed
(i.e. sepassoc = FALSE
) and a randomintercepts and randomslopes
model must be used (i.e. model = "intslope"
). Note that the function
only permits 2 failure types. The model is specified in full by Williamson
et al. (2008). The function joint()
automatically detects whether
competing risks are present by counting the number of unique components in
the event column on the event time data.
Both longsep
and survsep
ignore any
latent association (i.e. γ = 0) between the longitudinal and
survival processes but their output can be used to compare with the results
from the joint model. If interest is solely in the individual processes
then the user should instead make use of the functions
lme
and coxph
mentioned above.
Furthermore, if interest is in the separate effect of each random effect
(this is for intercept and slope or quadratic models only) upon the
survival data, the user should set sepassoc = TRUE
.
Since numerical integration is required, it is advisable to check the
stability of the maximum likelihood estimates with an increasing number of
GaussHermite quadrature points. joint()
uses gpt = 3
by
default, as this has been adequate for many datasets. However, for certain
datasets and models, this might be too small.
Pete Philipson
Wulfsohn MS, Tsiatis AA. A joint model for survival and longitudinal data measured with error. Biometrics. 1997; 53(1): 330339.
Henderson R, Diggle PJ, Dobson A. Joint modelling of longitudinal measurements and event time data. Biostatistics. 2000; 1(4): 465480.
Williamson PR, KolamunnageDona R, Philipson P, Marson AG. Joint modelling of longitudinal and competing risks data. Stat Med. 2008; 27: 64266438.
lme
, coxph
,
jointdata
, jointplot
.
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  ## Standard joint model
data(heart.valve)
heart.surv < UniqueVariables(heart.valve,
var.col = c("fuyrs", "status"),
id.col = "num")
heart.long < heart.valve[, c("num", "time", "log.lvmi")]
heart.cov < UniqueVariables(heart.valve,
c("age", "hs", "sex"),
id.col = "num")
heart.valve.jd < jointdata(longitudinal = heart.long,
baseline = heart.cov,
survival = heart.surv,
id.col = "num",
time.col = "time")
fit < joint(data = heart.valve.jd,
long.formula = log.lvmi ~ 1 + time + hs,
surv.formula = Surv(fuyrs, status) ~ hs,
model = "intslope")
## Competing risks joint model (same data as Williamson et al. 2008)
## Not run:
data(epileptic)
epileptic$interaction < with(epileptic, time * (treat == "LTG"))
longitudinal < epileptic[, c(1:3, 13)]
survival < UniqueVariables(epileptic, c(4, 6), "id")
baseline < UniqueVariables(epileptic, "treat", "id")
data < jointdata(longitudinal = longitudinal,
survival = survival,
baseline = baseline,
id.col = "id",
time.col = "time")
fit2 < joint(data = data,
long.formula = dose ~ time + treat + interaction,
surv.formula = Surv(with.time, with.status2) ~ treat,
longsep = FALSE, survsep = FALSE,
gpt = 3)
summary(fit2)
## End(Not run)

Loading required package: survival
Competing risks detected
Call:
joint(data = data, long.formula = dose ~ time + treat + interaction,
surv.formula = Surv(with.time, with.status2) ~ treat, longsep = FALSE,
survsep = FALSE, gpt = 3)
Random effects joint model
Data: data
Loglikelihood: 4056.936
Longitudinal submodel fixed effects: dose ~ time + treat + interaction
(Intercept) 1.9730179636
time 0.0003545439
treatLTG 0.1380959688
interaction 0.0006185222
Competing risks submodel effects: Surv(with.time, with.status2) ~ treat
Failure cause 1:
treatLTG 0.02721022
Failure cause 2:
treatLTG 0.6601433
Latent association:
gamma_1 0.5894143
gamma_2 0.9258597
Variance components:
U_0 U_1 Residual
7.321131e01 1.211086e06 1.964523e01
Convergence at iteration: 19
Number of observations: 2797
Number of groups: 605
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.