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 risks of the recurrent and terminal events are represented by proportional hazard risk models. The joint model is constructed assuming that the processes are linked via a latent structure (Krol et al. 2015):
$$\left{ \begin{array}{lc} 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)} \ r_{ij}(t|\bold{b}_i)=r_0(t)\exp(v_i+\bold{X}_{Rij}(t)\bold{\beta}_R+g(\bold{b}_i,\bold{\beta}_L,\bold{Z}_i(t),\bold{X}_{Li}(t))^\top \bold{\eta}_R ) & \mbox{(Recurrent)} \ \lambda_i(t|\bold{b}_i)=\lambda_0(t)\exp(\alpha v_i+\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)$, $\bold{X}_{Rij}(t)$ and $\bold{X}_{Ti}$ are vectors of fixed effects covariates and $\bold{\beta}_L$, $\bold{\beta}_R$ 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 biomarker and recurrent events is explained via $g(\bold{b}_i,\bold{\beta}_L,\bold{Z}_i(t),\bold{X}_{Li}(t))$ with coefficients $\bold{\eta}_R$ and between the biomarker and terminal event 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 functions $g(\cdot)$ and $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$. The frailty term $v_i$ is gaussian with mean 0 and variance $\sigma_v$. Together with $\bold{b}_i$ constitutes the random effects of the model: $$\bold{u}_i=\left(\begin{array}{c} \bold{b}_{i}\v_i \ \end{array}\right) \sim \mathcal{N}\left(\bold{0}, \left(\begin{array} {cc} \bold{B}_1&\bold{0} \ \bold{0} & \sigma_v^{2}\\end{array}\right)\right),$$
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).
trivPenal(formula, formula.terminalEvent, formula.LongitudinalData,
data, data.Longi, random, id, intercept = TRUE,
link = "Random-effects", left.censoring = FALSE,
recurrentAG = FALSE, n.knots, kappa, maxit = 350,
hazard = "Splines", init.B, init.Random, init.Eta, init.Alpha,
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" FALSE, otherwise the value of the threshold must be given."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 multivariateThe 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 transformation does not include the frailty in the trivariate model, for which the standard method is used). 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.
plot.trivPenal,print.trivPenal,summary.trivPenal###--- Trivariate joint model for longitudinal data, ---###
###--- recurrent events and a terminal event ---###
data(colorectal)
data(colorectalLongi)
# Parameter initialisation for covariates - longitudinal and terminal part
# Survival data preparation - only terminal events
colorectalSurv <- subset(colorectal, new.lesions == 0)
initial.longi <- 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 = 6, kappa = 2, method.GH="Pseudo-adaptive",
maxit=40, n.nodes=7)
# Parameter initialisation for covariates - recurrent part
initial.frailty <- frailtyPenal(Surv(time0, time1, new.lesions) ~ cluster(id)
+ age + treatment + who.PS, data = colorectal,
recurrentAG = TRUE, RandDist = "LogN", n.knots = 6, kappa =2)
# Baseline hazard function approximated with splines
# Random effects as the link function, Calendar timescale
# (computation takes around 40 minutes)
model.spli.RE.cal <-trivPenal(Surv(time0, time1, new.lesions) ~ cluster(id)
+ age + treatment + who.PS + terminal(state),
formula.terminalEvent =~ age + treatment + who.PS + prev.resection,
tumor.size ~ year * treatment + age + who.PS, data = colorectal,
data.Longi = colorectalLongi, random = c("1", "year"), id = "id",
link = "Random-effects", left.censoring = -3.33, recurrentAG = TRUE,
n.knots = 6, kappa=c(0.01, 2), method.GH="Standard", n.nodes = 7,
init.B = c(-0.07, -0.13, -0.16, -0.17, 0.42, #recurrent events covariates
-0.16, -0.14, -0.14, 0.08, 0.86, -0.24, #terminal event covariates
2.93, -0.28, -0.13, 0.17, -0.41, 0.23, 0.97, -0.61)) #biomarker covariates
# Weibull baseline hazard function
# Random effects as the link function, Gap timescale
# (computation takes around 30 minutes)
model.weib.RE.gap <-trivPenal(Surv(gap.time, new.lesions) ~ cluster(id)
+ age + treatment + who.PS + prev.resection + terminal(state),
formula.terminalEvent =~ age + treatment + who.PS + prev.resection,
tumor.size ~ year * treatment + age + who.PS, data = colorectal,
data.Longi = colorectalLongi, random = c("1", "year"), id = "id",
link = "Random-effects", left.censoring = -3.33, recurrentAG = FALSE,
hazard = "Weibull", method.GH="Pseudo-adaptive",n.nodes=7)Run the code above in your browser using DataLab