Description Usage Arguments Details Value Note Author(s) References See Also Examples
This function applies a maximum likelihood approach to fit the semiparametric joint models of survival and normal longitudinal data. The survival model is assumed to come from a class of transformation models, including the Cox proportional hazards model and the proportional odds model as special cases. The longitudinal process is modeled by nonparametric multiplicative random effects (NMRE) model.
1 2 
fitLME 
an object inheriting from class 
fitCOX 
an object inheriting from class 
data 
a data.frame containing all the variables included in the joint modeling. See Note. 
model 
an indicator specifying the dependency between the survival and longitudinal outcomes. Default is 1. See Details. 
rho 
a nonnegative real number specifying the transformation model you would like to fit. Default is 0, i.e. the Cox proportional hazards model. See Details. 
timeVarY 
a character string indicating the time variable in the NMRE model. See Examples. 
timeVarT 
a character string indicating the time variable in the 
control 
a list of control values for the estimation algorithm with components:

... 
additional options to be passed to the 
The jmodelMult
function fits joint models for survival and longitudinal data. Nonparametric multiplicative random effects models (NMRE) are assumed for the longitudinal processes. With the Cox proportional hazards model and the proportional odds model as special cases, a general class of transformation models are assumed for the survival processes. The baseline hazard functions are left unspecified, i.e. no parametric forms are assumed, thus leading to semiparametric models. For detailed model formulation, please refer to Xu, Hadjipantelis and Wang (2017).
The longitudinal model (NMRE) is written as
Y_i(t)=μ_i(t) + \varepsilon_i(t) = b_i\times\mathbf{B}^\top(t)\boldsymbolγ + \varepsilon_i(t),
where \mathbf{B}(t)=(B_1(t), \cdots, B_L(t)) is a vector of Bspline basis functions and b_i is a random effect \sim\mathcal{N}(1, σ_b^2). Note that we also allow the inclusion of baseline covariates as columns of \mathbf{B}(t). If model = 1
, then the linear predictor for the survival model is expressed as
η(t) = \mathbf{W}_i^\top(t)\boldsymbolφ + αμ_i(t),
indicating that the entire longitudinal process (free of error) enters the survival model as a covariate. If other values are assigned to the model
argument, the linear predictor for the surival model is then expressed as
η(t) = \mathbf{W}_i^\top(t)\boldsymbolφ + α b_i,
suggesting that the survival and longitudinal models only share the same random effect.
The survival model is written as
Λ(tη(t)) = G≤ft[\int_0^t\exp\{η(s)\}dΛ_0(s)\right],
where G(x) = \log(1 + ρ x) / ρ with ρ ≥q 0 is the class of logarithmic transfomrations. If rho = 0
, then G(x) = x, yielding the Cox proportional hazards model. If rho = 1
, then G(x) = \log(1 + x), yielding the proportional odds model. Users could assign any nonnegative real value to rho
.
An expectationmaximization (EM) algorithm is implemented to obtain parameter estimates. The convergence criterion is either of (i) \max \{  \boldsymbolθ^{(t)}  \boldsymbolθ^{(t  1)}  / (  \boldsymbolθ^{(t  1)}  + .Machine\$double.eps \times 2 ) \} < tol.P, or (ii)  L(\boldsymbolθ^{(t)})  L(\boldsymbolθ^{(t  1)}) / (  L(θ^{(t  1)})  + .Machine\$double.eps \times 2 ) < tol.L, is satisfied. Here \boldsymbolθ^{(t)} and \boldsymbolθ^{(t1)} are the vector of parameter estimates at the tth and (t1)th EM iterations, respectively; L(\boldsymbolθ) is the value of the loglikelihood function evaluated at \boldsymbolθ. Users could specify the tolerance values tol.P
and tol.L
through the control
argument.
For standard error estimation for the parameter estimates, three methods are provided, namely "PRES"
, "PFDS"
and "PLFD"
(detailed information are referred to Xu, Baines and Wang (2014)). In the control
argument, if SE.method = "PRES"
, numerically differentiating the profile Fisher score vector with Richardson extrapolation is applied; if SE.method = "PFDS"
, numerically differentiating the profile Fisher score vector with forward difference is applied; if SE.method = "PLFD"
, numerially (second) differentiating the profile likelihood with forward difference is applied. Generally, numerically differentiating a function f(x) (an arbitrary function) with forward difference is expressed as
f^\prime(x) = \frac{f(x + δ)  f(x)}{δ},
and that with Richardson extrapolation is expressed as
f^\prime(x) = \frac{f(x  2δ)  8f(x  δ) + 8f(x + δ)  f(x + 2δ)}{12δ}.
Users could specify the value of δ through the delta
item in the control
argument.
See jmodelMultObject
for the components of the fit.
1. To fit a nonparametric multiplicative random effects model, the fixed effect in the fitLME
object should be a matrix of Bspline basis functions (an object from the bs
function) with the possibility of including baseline covariates and the random effect should only include a random intercept. In the bs
function, it is a good practice to specify the boundary knots through the Boundary.knots
argument, where the upper boundary knot is typically the longest followup time among all subjects. See Examples.
2. Currently, jmodelMult()
could only handle the fitLME
object with a simple randomeffects structure (only the pdDiag()
class). Moreover, the withingroup correlation and heteroscedasticity structures in the fitLME
object (i.e. the correlation
and weights
argument of lme()
) are ignored.
3. The data
argument in jmodelMult()
, lme()
and coxph()
should be the same data frame.
4. For the fitCOX
object, only the W_i(t) in the linear predictor η(t) for the survial model (see Details) should be involved in the formula
argument of coxph{}
. Since coxph()
uses the same data frame as lme()
does, a timedependent Cox model must be fitted by coxph()
although W_i(t) may only contain timeindependent covariates. See Examples.
5. If W_i(t) in the linear predictor η(t) for the survial model (see Details) does involve timedependent covariate, then timeVarT
must specify the name of the time variable involved (see Examples).
6. The standard error estimates are obtained by numerical approximations which is naturally subject to numerical errors. Therefore, in extreme cases, there may be NA
values for part of the standard error estimates.
Cong Xu [email protected] Pantelis Z. Hadjipantelis [email protected]
Dabrowska, D. M. and Doksun K. A. (1988) Partial Likelihood in Transformation Models with Censored Data. Scandinavian Journal of Statistics 15, 1–23.
Ding, J. and Wang, J. L. (2008) Modeling longitudinal data with nonparametric multiplicative random effects jointly with survival data. Biometrics 64, 546–556.
Tsiatis, A. A. and Davidian, M. (2004) Joint modeling of longitudinal and timetoevent data: an overview. Statistica Sinica 14, 809–834.
Xu, C., Baines, P. D. and Wang, J. L. (2014) Standard error estimation using the EM algorithm for the joint modeling of survival and longitudinal data. Biostatistics 15, 731–744
Xu, C., Hadjipantelis, P. Z. and Wang, J. L. (2017) Semiparametric joint modeling of survival and longitudinal data: the R package JSM. Journal of Statistical Software, Under review.
Zeng, D. and Lin, D. (2007) Maximum likelihood estimation in semiparametric regression models with censored data. Journal of the Royal Statistical Society: Series B 69, 507–564.
jmodelMultObject
,
lme
,
coxph
,
Surv
,
bs
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  # linear mixedeffects model fit where the fixed effect is modeled by
# quadratic Bsplie basis with no internal knots
fitLME < lme(log(serBilir) ~ bs(obstime, degree = 2, Boundary.knots = c(0, 15)),
random = ~ 1  ID, data = pbc)
# Cox proportional hazards model fit with a single timeindependent covariate
fitCOX < coxph(Surv(start, stop, event) ~ drug, data = pbc, x = TRUE)
# joint model fit which assumes the Cox proportional hazards model for the survival process
# and NMRE for the longitudinal process. Use 'max.iter = 25', 'nknot = 3' and
# the 'PFDS' method to calculate standard error estimates as a quick toy example
fitJTMult.ph < jmodelMult(fitLME, fitCOX, pbc, timeVarY = "obstime",
control = list(SE.method = 'PFDS', max.iter = 25, nknot = 3))
summary(fitJTMult.ph)
## Not run:
# joint model fit with the default control
fitJTMult.ph2 < jmodelMult(fitLME, fitCOX, pbc, timeVarY = "obstime")
summary(fitJTMult.ph2)
# joint model fit where the survival and longitudinal processes only share
# the same random effect
fitJTMult.ph3 < jmodelMult(fitLME, fitCOX, pbc, model = 2, timeVarY = "obstime")
summary(fitJTMult.ph3)
# joint model fit which assumes the proportional odds model for the survival process
# and NMRE model for the longitudinal process
fitJTMult.po < jmodelMult(fitLME, fitCOX, pbc, rho = 1, timeVarY = "obstime")
summary(fitJTMult.po)
# joint model fit where the survival and longitudinal processes only share
# the same random effect
fitJTMult.po2 < jmodelMult(fitLME, fitCOX, pbc, model = 2, rho = 1, timeVarY = "obstime")
summary(fitJTMult.po2)
# allow baseline covariates in the NMRE model for the longitudinal process
fitLME2 < lme(log(serBilir) ~ drug + bs(obstime, degree = 2, Boundary.knots = c(0, 15)),
random = ~1  ID, data = pbc)
fitJTMult.ph4 < jmodelMult(fitLME2, fitCOX, pbc, timeVarY = "obstime")
summary(fitJTMult.ph4)
# Cox proportional hazards model fit with a timedependent covariate
fitCOX2 < coxph(Surv(start, stop, event) ~ drug + as.numeric(drug) : obstime,
data = pbc, x = TRUE)
# joint model fit in which \code{timeVarT} must be specified
fitJTMult.ph5 < jmodelMult(fitLME, fitCOX2, pbc, timeVarY = "obstime", timeVarT = 'obstime',
control = list(max.iter = 300))
summary(fitJTMult.ph5)
## End(Not run)

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