The longitudinal outcomes $y_i(t_{ik})$ ($k=1,\ldots,n_i$, $i=1,\ldots,N$) for $N$ subjects are described by a linear mixed model and the risk of the terminal event is represented by a proportional hazard risk model. The joint model is constructed assuming that the processes are linked via a latent structure (Wulfsohn and Tsiatsis 1997):
$$\left{ \begin{array}{ll} y_{i}(t_{ik})=\bold{X}_{Li}(t_{ik})^\top \bold{\beta}_L +\bold{ Z}_i(t_{ik})^\top \bold{b}_i + \epsilon_i(t_{ik}) & \mbox{(Longitudinal)} \ \lambda_i(t|\bold{b}_i)=\lambda_0(t)\exp(\bold{X}_{Ti}(t)\bold{\beta}_T+h(\bold{b}_i,\bold{\beta}_L,\bold{Z}_i(t),\bold{X}_{Li}(t))^\top \bold{\eta}_T ) & \mbox{(Terminal)} \ \end{array} \right.$$
where $\bold{X}_{Li}(t)$ and $\bold{X}_{Ti}$ are vectors of fixed effects covariates and $\bold{\beta}_L$ and $\bold{\beta}_T$ are the associated coefficients. Measurements errors $\epsilon_i(t_{ik})$ are iid normally distributed with mean 0 and variance $\sigma_{\epsilon}^2$. The random effects $\bold{b}_i = (b_{0i},\ldots, b_{qi})^\top\sim \mathcal{N}(0,\bold{B}_1)$ are associated to covariates $\bold{Z}_i(t)$ and independent from the measurement error. The relationship between the two processes is explained via $h(\bold{b}_i,\bold{\beta}_L,\bold{Z}_i(t),\bold{X}_{Li}(t))$ with coefficients $\bold{\eta}_T$. Two forms of the function $h(\cdot)$ are available: the random effects $\bold{b}_i$ and the current biomarker level $m_i(t)=\bold{X}_{Li}(t_{ik})^\top \bold{\beta}_L +\bold{ Z}_i(t_{ik})^\top \bold{b}_i$.
We consider that the longitudinal outcome can be a subject to a quantification limit, i.e. some observations, below a level of detection $s$ cannot be quantified (left-censoring).
longiPenal(formula, formula.LongitudinalData, data, data.Longi,
random, id, intercept = TRUE, link = "Random-effects",
left.censoring = FALSE, n.knots, kappa, maxit = 350,
hazard = "Splines", init.B, init.Random, init.Eta,
method.GH = "Standard", n.nodes, LIMparam = 1e-3,
LIMlogl = 1e-3, LIMderiv = 1e-3, print.times = TRUE)formula.formula.LongitudinalData."1".TRUE."Random-effects" for the association directly via the random effects of the biomarker, "Current-level" for the association via the true current level of thFALSE, ie. no left-censoring. In case of a left-censored biomarker, this argument must be equal to the threshold $s$."Splines" for semiparametric hazard functions using equidistant intervals or "Splines-per" using percentile with the penalized likelihood estimation,
"Weibull" for the parametric "Standard" for the standard non-adaptive Gaussian quadrature, "Pseudo-adaptive" for the pseudo-adaptive Gaussian quadrature and "HRMSYM" for the algorithm for the multivariatefrailtyPenal function), $10^{-3}$ by default.frailtyPenal function), $10^{-3}$ by default.frailtyPenal function), $10^{-3}$ by default.The method of the Gauss-Hermite quadrature for approximations of the multidimensional integrals, i.e. length of random is 2, can be chosen among the standard, non-adaptive, pseudo-adaptive in which the quadrature points are transformed using the information from the fitted mixed-effects model for the biomarker (Rizopoulos 2012) or multivariate non-adaptive procedure proposed by Genz et al. 1996 and implemented in FORTRAN subroutine HRMSYM.
The choice of the method is important for estimations. The standard non-adaptive Gauss-Hermite quadrature ("Standard") with a specific number of points gives accurate results but can be time consuming. The non-adaptive procedure ("HRMSYM") offers advantageous computational time but in case of datasets in which some individuals have few repeated observations (biomarker measures or recurrent events), this method may be moderately unstable.
The pseudo-adaptive quadrature uses transformed quadrature points to center and scale the integrand by utilizing estimates of the random effects from an appropriate linear mixed-effects model. This method enables using less quadrature
points while preserving the estimation accuracy and thus lead to a better computational time.
NOTE. Data frames data and data.Longi must be consistent. Names and types of corresponding covariates must be the same, as well as the number and identification of individuals.
A. Krol, L. Ferrer, JP. Pignon, C. Proust-Lima, M. Ducreux, O. Bouche, S. Michiels, V. Rondeau (2015). Joint Model for Left-Censored Longitudinal Data, Recurrent Events and Terminal Event: Predictive Abilities of Tumor Burden for Cancer Evolution with Application to the FFCD 2000-05 Trial. Submitted.
D. Rizopoulos (2012). Fast fitting of joint models for longitudinal and event time data using a pseudo-adaptive Gaussian quadrature rule. Computational Statistics and Data Analysis 56, 491-501.
M.S. Wulfsohn, A.A. and Tsiatis, A. A. (1997). A joint model for survival and longitudinal data measured with error. Biometrics 53, 330-9.
plot.longiPenal,print.longiPenal,summary.longiPenal###--- Joint model for longitudinal data and a terminal event ---###
data(colorectal)
data(colorectalLongi)
# Survival data preparation - only terminal events
colorectalSurv <- subset(colorectal, new.lesions == 0)
# Baseline hazard function approximated with splines
# Random effects as the link function
model.spli.RE <- longiPenal(Surv(time1, state) ~ age + treatment + who.PS
+ prev.resection, tumor.size ~ year * treatment + age + who.PS ,
colorectalSurv, data.Longi = colorectalLongi, random = c("1", "year"),
id = "id", link = "Random-effects", left.censoring = -3.33,
n.knots = 7, kappa = 2)
# Weibull baseline hazard function
# Current level of the biomarker as the link function
model.weib.CL <- longiPenal(Surv(time1, state) ~ age + treatment + who.PS
+ prev.resection, tumor.size ~ year * treatment + age + who.PS ,
colorectalSurv, data.Longi = colorectalLongi, random = c("1", "year"),
id = "id", link = "Current-level", left.censoring = -3.33, hazard = "Weibull")Run the code above in your browser using DataLab