ISPED Agnieszka Król

Audrey Mauguen

Yassin Mazroui

Alexandre Laurent

Stefan Michiels

Virginie Rondeau

Tutorial in Joint Modeling and Prediction: A Statistical Software for Correlated Longitudinal Outcomes, Recurrent Events and a Terminal Event

Keywords: dynamic prediction, frailty, joint model, longitudinal data, predictive accuracy, R, survival analysis

Extensions in the field of joint modeling of correlated data and dynamic predictions improve the development of prognosis research. The R package frailtypack provides estimations of various joint models for longitudinal data and survival events. In particular, it fits models for recurrent events and a terminal event (frailtyPenal), models for two survival outcomes for clustered data (frailtyPenal), models for two types of recurrent events and a terminal event (multivPenal), models for a longitudinal biomarker and a terminal event (longiPenal) and models for a longitudinal biomarker, recurrent events and a terminal event (trivPenal). The estimators are obtained using a standard and penalized maximum likelihood approach, each model function allows to evaluate goodness-of-fit analyses and provides plots of baseline hazard functions. Finally, the package provides individual dynamic predictions of the terminal event and evaluation of predictive accuracy. This paper presents the theoretical models with estimation techniques, applies the methods for predictions and illustrates frailtypack functions details with examples.

Introduction

Joint models

Recent technologies allow the registration of greater and greater amount of data. In medical research different kinds of patient information are gathered over time together with clinical outcome data such as overall survival (OS). Joint models enable the analysis of correlated data of different types such as individual repeated data and clustered data together with OS. The repeated data may be recurrent events (e.g., relapses of a tumor, hospitalizations, appearance of new lesions) or a longitudinal outcome called biomarker (e.g., tumor size, prostate-specific antigen or CD4 cell counts). The correlated data that are not analyzed jointly with OS, are subjugated to a bias coming from ignoring the terminal event which is the competing event for the occurrence of repeated outcomes (not only does it preclude the outcomes from being observed but it also prevents them from occurring). On the other hand, the standard survival analysis for OS may lead to biased estimations, if the repeated data is considered as time-varying covariates or if it is completely ignored in the analysis.

For these correlated data one can consider joint models, e.g., a joint model for a longitudinal biomarker and a terminal event, which received most of the attention in the literature. This joint model estimates simultaneously the longitudinal and survival processes using the relationship via a latent structure of random effects [START_REF] Wulfsohn | A Joint Model for Survival and Longitudinal Data Measured with Error[END_REF]. Extensions of these include, among others, models for multiple longitudinal outcomes [START_REF] Hatfield | Multilevel Bayesian Models for Survival Times and Longitudinal Patient-Reported Outcomes with Many Zeros[END_REF], multiple failure times [START_REF] Elashoff | A Joint Model for Longitudinal Measurements and Survival Data in the Presence of Multiple Failure Types[END_REF] and both [START_REF] Chi | Joint Models for Multivariate Longitudinal and Multivariate Survival Data[END_REF]. A review of the joint modeling of longitudinal and survival data was already given elsewhere [START_REF] Mccrink | Advances in Joint Modelling: A Review of Recent Developments with Application to the Survival of End Stage Renal Disease Patients[END_REF][START_REF] Gould | Joint Modeling of Survival and Longitudinal Non-Survival Data: Current Methods and Issues. Report of the DIA Bayesian Joint Modeling Working Group[END_REF][START_REF] Asar | Joint Modelling of Repeated Measurement and Time-to-Event Data: An Introductory Tutorial[END_REF].

Another option for the analysis of correlated data are joint models for recurrent events and a terminal event [START_REF] Liu | Shared Frailty Models for Recurrent Events and a Terminal Event[END_REF][START_REF] Rondeau | Joint Frailty Models for Recurring Events and Death Using Maximum Penalized Likelihood Estimation: Application on Cancer Events[END_REF]. These models are usually called joint frailty models as the processes are linked via a random effect that represents the frailty of a subject (patient) to experience an event. They account for the dependence between two survival processes quantified by the frailty term, contrary to the alternative marginal approach [START_REF] Li | Use of the Wei-Lin-Weissfeld Method for the Analysis of a Recurring and a Terminating Event[END_REF], which does not model the dependence. Some extensions to joint frailty models include incorporation of a nonparametric covariate function [START_REF] Yu | A Joint Model of Recurrent Events and a Terminal Event with a Nonparametric Covariate Function[END_REF], inclusion of two frailty terms for the identification of the origin of the dependence between the processes [START_REF] Mazroui | General Joint Frailty Model for Recurrent Event Data with a Dependent Terminal Event: Application to Follicular Lymphoma Data[END_REF], consideration of the diseasespecific mortality process [START_REF] Belot | A Joint Frailty Model to Estimate the Recurrence Process and the Disease-Specific Mortality Process without Needing the Cause of Death[END_REF] and accommodation of time-varying coefficients [START_REF] Yu | Joint Model of Recurrent Events and a Terminal Event With Time-Varying Coefficients[END_REF][START_REF] Mazroui | Time-Varying Coefficients in a Multivariate Frailty Model: Application to Breast Cancer Recurrences of Several Types and Death[END_REF]. Finally, [START_REF] Mazroui | Multivariate Frailty Models for Two Types of Recurrent Events with a Dependent Terminal Event: Application to Breast Cancer Data[END_REF] proposed a model with two types of recurrent events following the approach of [START_REF] Zhao | Analysis of Multivariate Recurrent Event Data with Time-Dependent Covariates and Informative Censoring[END_REF]. A review of joint frailty models in the Bayesian context was given by [START_REF] Sinha | Current Methods for Recurrent Events Data with Dependent Termination: A Bayesian Perspective[END_REF].

Joint models for recurrent events and longitudinal data have received the least attention among the joint models so far. However, a marginal model based on the generalized linear mixed model was proposed by [START_REF] Efendi | A Joint Model for Longitudinal Continuous and Time-to-Event Outcomes with Direct Marginal Interpretation[END_REF]. This model allows several longitudinal and time-to-event outcomes and includes two sets of random ef-fects, one for the correlation within a process and between the processes, and the second to accommodate for overdispersion.

Finally, in some applications all the types of data, i.e., a longitudinal biomarker, recurrent events and a terminal event can be studied so that all of them are correlated to each other (in the following we call such models trivariate models). [START_REF] Liu | Analysis of Longitudinal Data in the Presence of Informative Observational Times and a Dependent Terminal Event, with Application to Medical Cost Data[END_REF] proposed a trivariate model for medical cost data. The longitudinal outcomes (amount of medical costs) were reported at informative times of recurrent events (hospitalizations) and were related to a terminal event (death). The dependence via random effects was introduced so that the process of longitudinal measurements and the process of recurrent events were related to the process of the terminal event. A relationship between longitudinal outcomes and recurrent events was not considered. This relationship was incorporated into a trivariate model proposed by [START_REF] Liu | Joint Analysis of Correlated Repeated Measures and Recurrent Events Processes in the Presence of Death, with Application to a Study on Acquired Immune Deficiency Syndrome[END_REF] for an application of an HIV study. In this parametric approach the focus was on the analysis of the associations present in the model and the effect of the repeated measures and recurrent events on the terminal event. [START_REF] Kim | Joint Models of Longitudinal Data and Recurrent Events with Informative Terminal Event[END_REF] analyzed the trivariate data using the transformation functions for the cumulative intensities for recurrent and terminal events. Finally, [START_REF] Król | 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[END_REF] proposed a trivariate model in which all the processes were related to each other via a latent structure of random effects and applied the model to the context of a cancer evolution and OS.

Prediction in joint models

In the framework of joint models that consider a terminal event, one can be interested in predictions of the event derived from the model. As joint models include information from repeated outcomes, these predictions are dynamic and they change with the update of the observations. Dynamic predictive tools were proposed for joint models for longitudinal data and a terminal event [START_REF] Proust-Lima | Development and Validation of a Dynamic Prognostic Tool for Prostate Cancer Recurrence Using Repeated Measures of Posttreatment PSA: A Joint Modeling Approach[END_REF][START_REF] Rizopoulos | Dynamic Predictions and Prospective Accuracy in Joint Models for Longitudinal and Time-to-Event Data[END_REF], joint models for recurrent events and a terminal event [START_REF] Mauguen | Dynamic Prediction of Risk of Death Using History of Cancer Recurrences in Joint Frailty Models[END_REF] and trivariate models [START_REF] Król | 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[END_REF].

For the evaluation of the predictive accuracy of a joint model few methods were proposed due to the complexity of the models in which the survival data are usually right-censored. The standard methods for the assessment of the predictive abilities derived from the survival analysis and adapted to the joint models context are the Brier score (Proust-Lima and Taylor 2009; [START_REF] Mauguen | Dynamic Prediction of Risk of Death Using History of Cancer Recurrences in Joint Frailty Models[END_REF][START_REF] Król | 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[END_REF]) and the receiver operating characteristic (ROC) methodology [START_REF] Rizopoulos | Dynamic Predictions and Prospective Accuracy in Joint Models for Longitudinal and Time-to-Event Data[END_REF][START_REF] Blanche | Quantifying and Comparing Dynamic Predictive Accuracy of Joint Models for Longitudinal Marker and Time-to-Event in Presence of Censoring and Competing Risks[END_REF]. Finally, the expected prognostic cross-entropy (EPOCE), a measure based on information theory, provides a useful tool for the evaluation of predictive abilities of a model [START_REF] Commenges | Choice of Prognostic Estimators in Joint Models by Estimating Differences of Expected Conditional Kullback-Leibler Risks[END_REF][START_REF] Proust-Lima | Joint Latent Class Models for Longitudinal and Time-to-Event Data: A Review[END_REF].

Software for joint models

Together with the theoretical development of the joint models, an increase of appropriate software is observed, mostly for standard frameworks. Among the available R (R Core Team 2017) packages, the joint analysis of a single longitudinal outcome and a terminal event can be performed using JM [START_REF] Rizopoulos | JM: An R Package for the Joint Modelling of Longitudinal and Timeto-Event Data[END_REF][START_REF] Rizopoulos | JM: Joint Modeling of Longitudinal and Survival Data. R package version 1.4[END_REF] based on the likelihood maximization approach using an expectation-maximization (EM) algorithm, JMbayes (Rizopoulos 2016b,a) built in the Bayesian context, and joineR [START_REF] Philipson | joineR: Joint Modelling of Repeated Measurements and Time-to-Event Data[END_REF], a frequentist approach that allows a flexible formulation of the dependence structure. For the approach based on joint latent class models, i.e., joint models that consider homogeneous latent subgroups of individuals sharing the same biomarker trajectory and risk of the terminal event, the extensive package lcmm (Proust-Lima, Philipps, and Liquet 2017b;Proust-Lima, Philipps, Diakite, and Liquet 2017a) provides estimations based on the frequentist approach. Apart from the R packages, a stjm package [START_REF] Crowther | stjm: Stata Module to Fit Shared Parameter Joint Models of Longitudinal and Survival Data[END_REF] in Stata (StataCorp. 2015) uses maximum likelihood estimation and provides flexible methods for modeling the longitudinal outcome based on polynomials or splines. Another possibility for the analysis of the joint models for a longitudinal outcome and a terminal event is the SAS (SAS Institute Inc. 2013) macro JMFit [START_REF] Zhang | JMFit: A SAS Macro for Joint Models of Longitudinal and Survival Data[END_REF]. Joint models using nonlinear mixed effects models can be estimated using the MONOLIX software [START_REF] Mbogning | Joint Modeling of Longitudinal and Repeated Time-to-Event Data Using Nnonlinear Mixed-Eeffects Models and the SAEM Aalgorithm[END_REF]. The packages JM, JMbayes and lcmm also provide predictive tools and predictive accuracy measures: the EPOCE estimator in lcmm, ROC and AUC (area under the curve) measures and prediction error in JM and JMbayes.

For joint models for correlated events (recurrent event or clustered event data) and a terminal event the available statistical software is limited. Among the available R packages, joint.Cox [START_REF] Emura | Cox: Penalized Likelihood Estimation under the Joint Cox Models Between TTP and OS for Meta-Analysis[END_REF]) provides the estimations using penalized likelihood under the copula-based joint Cox models for time to clustered events (progression) and a terminal event. Trivariate joint models proposed by [START_REF] Liu | Joint Analysis of Correlated Repeated Measures and Recurrent Events Processes in the Presence of Death, with Application to a Study on Acquired Immune Deficiency Syndrome[END_REF] were performed in SAS using the procedure NLMIXED.

The R package frailtypack [START_REF] Rondeau | frailtypack: General Frailty Models: Shared, Joint and Nested Frailty Models with Prediction[END_REF]) fits several types of joint models. Originally developed for the shared frailty models for correlated outcomes it has been extended into the direction of joint models. These extensions include the joint model for one or two types of recurrent events and a terminal event, the joint model for a longitudinal biomarker and a terminal event and the trivariate model for a longitudinal biomarker, recurrent events and a terminal event. Moreover, frailtypack includes prediction tools, such as concordance measures adapted to shared frailty models and predicted probability of events for the joint models. The package is available from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=frailtypack.

Characteristics of a previous version of the package, such as estimation of several shared and standard joint frailty models using penalized likelihood, were already described elsewhere [START_REF] Rondeau | frailtypack: An R Package for the Analysis of Correlated Survival Data with Frailty Models Using Penalized Likelihood Estimation or Parametrical Estimation[END_REF]. This work focuses on an overview and the developments of joint models included in the package (models for recurrent/clustered events and a terminal event, models for multivariate recurrent events and a terminal event, models for longitudinal data and a terminal event and models for longitudinal data, recurrent events and a terminal event) and the prediction tools accompanied by predictive accuracy measures. Finally, several options available for the models (e.g., two correlated frailties in the model for recurrent events and a terminal event, left-censoring for the longitudinal data) will be presented with the appropriate examples. A practical guide to different types of models included in the package along with available options is presented in a schematic table in Appendix A. This article firstly presents joint models with details on estimation methods (Section 2) and predictions of an event in the framework of these joint models (Section 3). In Section 4, the frailtypack package features are detailed. Section 5 contains some examples on real datasets and finally, Section 6 concludes the paper.

Models for correlated outcomes

Bivariate joint frailty models for two types of survival outcomes

Joint model for recurrent events and terminal event

We denote by i = 1, . . . , N the subject and by j = 1, . . . , n i the rank of the recurrent event.

For each subject i we observe the time of the terminal event T i = min(T * i , C i ) with T * i the true terminal event time and C i the censoring time. We denote the observed recurrent times

T ij = min(T * ij , C i , T * i )
with T * ij the real time of the recurrent event. Thus, for each rank j we can summarize the information with a triple {T ij , δ ij , δ i }, where

δ ij = I {T ij =T * ij } and δ i = I {T i =T * i } .
Let r 0 (.) and λ 0 (.) be baseline hazard functions related to the recurrent and terminal events, respectively. Let X Rij and X T i be two vectors of covariates associated with the risk of the related events. Let β R and β T be constant effects of the covariates whereas β R (t) and β T (t) denote a time-dependent effect. Finally, a frailty u i is a random effect that links the two risk functions and follows a given distribution D and α is a parameter allowing more flexibility. The hazard functions are defined by [START_REF] Liu | Shared Frailty Models for Recurrent Events and a Terminal Event[END_REF][START_REF] Rondeau | Joint Frailty Models for Recurring Events and Death Using Maximum Penalized Likelihood Estimation: Application on Cancer Events[END_REF]):

r ij (t|u i ) = u i r 0 (t) exp(X Rij β R ) = u i r ij (t) (recurrent event) λ i (t|u i ) = u α i λ 0 (t) exp(X T i β T ) = u α i λ i (t) (terminal event) (1)
where the frailty terms u i are iid (independent and identically distributed) and follow either the Gamma distribution with unit mean and variance θ (θ > 0), i.e., u i ∼ G( 1 θ , 1 θ ), or the lognormal distribution, i.e., u i = exp(v i ) ∼ lnN (0, σ 2 v ). The parameter α determines direction of the association (if found significant) between the processes.

For a given subject we can summarize all the information with

Θ i = {T (1) i , T i , δ (1) i , δ i }, where T (1) i = {T ij , j = 1, . . . , n i } and δ (1) i = {δ ij , j = 1, . . . , n i }. Finally, we are interested to estimate ξ = {r 0 (•), λ 0 (•), β R , β T , θ, α}.
In some cases, e.g., long follow-up, effects of certain prognostic factors can be varying in time. For this reason a joint model with time-dependent coefficients was proposed by [START_REF] Mazroui | Time-Varying Coefficients in a Multivariate Frailty Model: Application to Breast Cancer Recurrences of Several Types and Death[END_REF]. The coefficients β R (t) and β T (t) are modeled with regression B-splines of order q and m interior knots. The general form of an estimated time-dependent coefficient β(t) is

β(t) = m j=-q+1 ζ j B j,q (t),
where B j,q (t) is the basis of B-splines calculated using a recurring expression (De Boor 2001). Therefore, for every time-varying coefficient, m + q parameters ζ j , j = -q + 1, . . . , m need to be estimated. It has been shown that quadratic B-splines, i.e., q = 3, with a small number of interior knots (m ≤ 5) ensure stable estimation [START_REF] Mazroui | Time-Varying Coefficients in a Multivariate Frailty Model: Application to Breast Cancer Recurrences of Several Types and Death[END_REF]. The pointwise confidence intervals for β(t) can be obtained using the Hessian matrix.

The application of time-varying coefficients is helpful if there is a need to verify the proportional hazard (PH) assumption. Using a likelihood ratio (LR) test, the time-dependency of a covariate can be examined. A model with the time-dependent effect and a model with the constant effect are compared and the test hypotheses depend on whether the covariate is related to one of the events or both. If the time-dependency is tested for both events the null hypothesis is

H 0 : β R (t) = β R , β T (t) = β T and the alternative hypothesis is H 1 : β R (t) = β R , β T (t) = β T . The LR statistic has a χ 2 distribution of degree k(m + q -1),
where k is the number of events to which the covariate is related.

The LR test can also be used to verify whether a covariate with the time-dependent effect is significant for the events. In this case, a model with the time-varying effect covariate is compared to a model without this covariate. Again, the null hypotheses depend on the events considered. If the covariate is associated to both events, the null hypothesis

H 0 : β R (t) = 0, β T (t) = 0 is tested against H 1 : β R (t) = 0, β T (t) = 0. The LR statistic has a χ 2 distribution of degree k(m + q).

Joint model for two clustered survival events

An increasing number of studies favor the presence of clustered data, especially multi-center or multi-trial studies. In particular, meta-analyses are used to calculate surrogate endpoints in oncology. However, the clustering creates some heterogeneity between subjects, which needs to be accounted for. Using the above notations, the clustered joint model is similar to the model (1) but here, the index j (j = 1, . . . , n i ) represents a subject from cluster i (i = 1, . . . , N ) and the cluster-specific frailty term u i is shared by the subjects of a given group. Thus, the model can be written as [START_REF] Rondeau | A Joint Model for the Dependence between Clustered Times to Tumour Progression and Deaths: A Meta-Analysis of Chemotherapy in Head and Neck Cancer[END_REF]:

r ij (t|u i ) = u i r 0 (t) exp(X Rij β R ) = u i r ij (t) (time to event 1) λ ij (t|u i ) = u α i λ 0 (t) exp(X T ij β T ) = u α i λ ij (t) (time to event 2)
, and we assume the frailty terms u i to be iid Gamma distributed. The events can be chosen arbitrarily but it is assumed that the event 2 impedes the process of the event 1. Usually the event 2 is death of a patient and the other is an event of interest (e.g., surrogate endpoint for OS) such as time to tumor progression or progression-free survival.

The interest of using the joint model for the clustered data stems from the fact that it considers the dependency between the survival processes and respects that event 2 is a competitive event for event 1. The frailty term u i is common for a given group and represents the clustered association between the processes (at the cluster level) as well as the intra-cluster correlation.

The package frailtypack includes clustered joint models for two survival outcomes in presence of semi-competing risks (time-to-event, recurrent events are not allowed). To ensure identifiability of the models, it is assumed that α = 1. For these models, only a Gamma distribution of the frailty term is implemented in the package.

Joint model for recurrent events and a terminal event with two frailty terms

In model (1) the frailty term u i reflects the inter-and intra-subject correlation for the recurrent event as well as the association between the recurrent and the terminal events. In order to distinguish the origin of dependence, two independent frailty terms u i = (u i , v i ) can be considered [START_REF] Mazroui | General Joint Frailty Model for Recurrent Event Data with a Dependent Terminal Event: Application to Follicular Lymphoma Data[END_REF]):

r ij (t|u i ) = u i v i r 0 (t) exp(X Rij β R ) = u i v i r ij (t) (recurrent event) λ i (t|u i ) = u i λ 0 (t) exp(X T i β T ) = u i λ i (t) (terminal event) , ( 2
)
where

v i ∼ Γ( 1 η , 1 η ) (η > 0)
specific to the recurrent event rate and u i ∼ Γ( 1 θ , 1 θ ) (θ > 0) specific to the association between the processes. The variance of the frailty terms represents the heterogeneity in the data, associated with unobserved covariates. Moreover, a high value of the variance η indicates strong dependence between the recurrent events and a high value of the variance θ indicates that the recurrent and the terminal events are strongly dependent. The individual information is equivalent to Θ i from model (1) and the parameters to estimate are ξ = {r 0 (•), λ 0 (•), β R , β T , θ, η}.

In the package frailtypack for model (2), called also general joint frailty model, only the Gamma distribution is allowed for the random effects.

Multivariate joint frailty model

One of the extensions to the joint model for recurrent events and a terminal event (1) is the consideration of different types of recurrence processes and time-to-event data. Two types of recurrent events are taken into account in a multivariate frailty model proposed by [START_REF] Mazroui | Multivariate Frailty Models for Two Types of Recurrent Events with a Dependent Terminal Event: Application to Breast Cancer Data[END_REF]. The aim of the model is to analyze dependencies between all types of events. The recurrent event times are defined by

T (1) ij (j = 1, . . . , n (1) i ) and T (2) ij (j = 1, . . . , n
(2) i ) and both processes are censored by the terminal event T i . The joint model is expressed using the recurrent events and terminal event hazard functions:

       r (1) ij (t|u i , v i ) = r (1) 0 (t) exp(u i + X (1) Rij β (1) R ) = exp (u i )r (1) ij (t) (rec. event 1) r (2) ij (t|u i , v i ) = r (2) 0 (t) exp(v i + X (2) Rij β (2) R ) = exp(v i )r (2) ij (t) (rec. event 2) λ i (t|u i , v i ) = λ 0 (t) exp(α 1 u i + α 2 v i + X T i β T ) = exp(α 1 u i + α 2 v i )λ i (t) (terminal event) , (3) with vectors of regression coefficients β (1) R , β (2) R , β T and covariates X (1) Rij , X (2)
Rij , X T i for the first recurrent event, the second recurrent event and the terminal event, respectively. The frailty terms u i = (u i , v i ) explain the intra-correlation of the processes and the dependencies between them. For these correlated random effects the multivariate normal distribution of dimension 2 is considered:

u i = u i v i ∼ N 0, σ 2 u ρσ u σ v ρσ u σ v σ 2 v .
The variances σ 2 u and σ 2 v explain the within-subject correlation between occurrences of the recurrent event of type 1 and type 2, respectively. The dependency between the two recurrent events is explained by the correlation coefficient ρ and the dependency between the recurrent event 1 (event 2) and the terminal event is assessed by the term α 1 (α 2 ) in case of significant variance σ 2 u (σ 2 v ). Here, for each subject i we observe

Θ i = {T (1) i , T (2) i , T i , δ (1) i , δ (2) i , δ i }, where T (l) i = {T (l) ij , j = 1, . . . , n (l) i } and δ (l) i = {δ (l) ij , j = 1, . . . , n (l) i }, l = 1,
(•), r (2) 0 (•), λ 0 (•), β (1) R , β (2) R , β T , σ 2 u , σ 2 v , ρ, α 1 , α 2 }.

Bivariate joint model with longitudinal data

Here, instead of recurrent events we consider a longitudinal biomarker. For subject i we observe an l i -vector of longitudinal measurements y i = {y i (t ik ), k = 1, . . . , l i }. Again, the true terminal event time T * i impedes the longitudinal process and the censoring time does not stop it but values of the biomarker are no longer observed. For each subject i we observe

Θ i = {y i , T i , δ i }.
A random variable Y i (t) representing the biomarker is expressed using a linear mixed effects model and the terminal event time T * i using a proportional hazards model. The sub-models are linked by the random effects u i = b i :

Y i (t) = m i (t) + i (t) = X Li (t) β L + Z i (t) b i + i (t) (biomarker) λ i (t|b i ) = λ 0 (t) exp(X T i β T + h(b i , β L , Z i (t), X Li (t)) η T ) (death) , ( 4
)
where X Li (t) and X T i are vectors of fixed effects covariates. Coefficients β L and β T are constant in time. Measurements errors i (•) are iid normally distributed with mean 0 and variance σ 2 . We consider a q-vector of the random effects parameters b i = (b 0 , . . . , b q-1 ) ∼ N (0, B 1 ) associated to covariates Z i (t) and independent from the measurement error. In the frailtypack package, the maximum size of the matrix B 1 is 3. The relationship between the two processes is explained via a function h(b i , β L , Z i (t), X Li (t)) and quantified by coefficients η T . This possibly multivariate function represents the prognostic structure of the biomarker m i (t) (we assume that the measurement errors are not prognostic for the survival process) and in the frailtypack these are either the random effects b i or the current biomarker level m i (t). The structure of dependence is chosen a priori and should be designated with caution as it influences the model in terms of fit and predictive accuracy.

We consider that the longitudinal outcome y i (t ik ) can be subject to a quantification limit, i.e., some observations, below a level of detection s cannot be quantified (left-censoring). We introduce this in the vector y i which includes l o i observed values y o i and l c i censored values y c i

(l i = l o i + l c i ).
The aspect of the left-censored data is handled in the individual contribution from the longitudinal outcome to the likelihood. The individual observed outcomes are Θ i = {y o i , y c i , T i , δ i } (in case of no left-censoring y c i is empty and y o i is equivalent to y i ) and the parameters to estimate are ξ = {λ 0 (•), β L , β T , B 1 , σ 2 , η T }.

Trivariate joint model with longitudinal data

The trivariate joint model combines the joint model for recurrent events and a terminal event (model (1)) with the joint model for longitudinal data and a terminal event (4). We consider the longitudinal outcome Y observed in discrete time points y i = {y i (t ik ), k = 1, . . . , l i }, times of the recurrent event T ij (j = 1, . . . , n i ) and the terminal event time T i that is an informative censoring for the longitudinal data and recurrences. The respective processes are linked to each other via a latent structure. We define multivariate functions g(•) and h(•) for the associations between the biomarker and the recurrent events and between the biomarker and the terminal event, respectively. For the link between the recurrent events and the terminal event we use a frailty term v i from the normal distribution. This distribution was chosen here to facilitate the estimation procedure involving multiple numerical integration. Interpretation of v i should be performed with caution as the dependence between the recurrences and the terminal data is partially explained by the random effects of the biomarker trajectory. We define the model by [START_REF] Król | 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[END_REF]:

     Y i (t) = m i (t) + i (t) = X Li (t) β L + Z i (t) b i + i (t) (biomarker) r ij (t|v i , b i ) = r 0 (t) exp(v i + X Rij β R + g(b i , β L , Z i (t), X Li (t)) η R ) (recurrent event) , λ i (t|v i , b i ) = λ 0 (t) exp(αv i + X T i β T + h(b i , β L , Z i (t), X Li (t)) η T ) (terminal event)
(5) where the regression coefficients β L , β R , β T are associated to possibly time-varying covariates X Li (t), X Rij , X T i for the biomarker, recurrent events and the terminal event (baseline prognostic factors), respectively. The strength of the dependencies between the processes is quantified by α and vectors η R and η T . The vector of the random effects u i = (b i , v i ) of dimension q = q 1 + 1 follows the multivariate normal distribution:

u i = b i v i ∼ N 0, B 1 0 0 σ 2 v ,
where the dimension of the matrix B 1 in the frailtypack package is allowed to be maximum 3. It is convenient if one assume for the biomarker a random intercept and slope:

Y i (t) = m i (t) + i (t) = X Li (t) β L + b i0 + b i1 × t + i (t)
and then b i0 represents the heterogeneity of the baseline measures of the biomarker among the subjects and b i1 the heterogeneity of the slope of the biomarker's linear trajectory among the subjects. For non-linear trajectories, one can use parametric functions to assume two slopes of time, e.g., initial drop and long-term regrowth of tumor size under a treatment.

As for model (4) we allow the biomarker to be left-censored, again with y o i the observed outcomes and y c i the undetected ones. For individual i we observe then

Θ i = {y o i , y c i , T (1) i , T i , δ (1) i , δ i } and the interest is to estimate ξ = {r 0 (•), λ 0 (•), β L , β R , β T , B 1 , σ 2 v , σ 2 , α, η R , η T }.

Estimation

Maximum likelihood estimation

Estimation of a model's parameters ξ is based on the maximization of the marginal loglikelihood derived from the joint distribution of the observed outcomes that are assumed to be independent from each other given the random effects. Let u i represent random effects of a joint model (u i can be a vector, for models (2), (3), (4), (5) or a scalar, for model (1)) and f u i (u i ; ξ) the density function of the distribution of u i . The marginal individual likelihood is integrated over the random effects and is given by:

L i (Θ i ; ξ) = u i f y i |u i (y i |u i ; ξ) γ L f T (1) i |u i (T (1)
i , δ

(1)

i |u i ; ξ) γ R (1) f T (2) i |u i (T (2) i , δ (2) i |u i ; ξ) γ R (2) × f T i |u i (T i , δ i |u i ; ξ)f u i (u i ; ξ) du i , ( 6
)
where indicators γ • are introduced so that the likelihood is valid for all the joint models (1-5). Therefore, γ L = 1 in case of models (4) and ( 5) and 0 otherwise, γ R (1) = 1 in case of the models (1), (3) and ( 5) and 0 otherwise, γ R (2) = 1 in case of the model (3) and 0 otherwise. The conditional density of the longitudinal outcome f y i |u i is the density of the l i -dimensional normal distribution with mean m i = {m i (t ik ), k = 1, . . . , l i } and variance σ 2 I m i . In case of left-censored data, we observe l o i outcomes y o i and l c i outcomes y c i are left-censored (below the threshold s). Then f y i |u i can be written as a product of the density for the observed outcomes (normal with mean m i (t ik ) and variance σ 2 ) and the corresponding cumulative distribution function (cdf) F y i |u i for the censored outcomes:

f y i |u i (y i |u i ; ξ) = l o i k=1 f y o i (t ik )|u i (y o i (t ik )|u i ; ξ) l c i k=1 F y c i (t ik )|u i (s|u i ).
Modeling of the risk of an event can be performed either in a calendar timescale or in a gap timescale. This choice depends on the interest of the application, whether the focus is on time from a fixed date, e.g., beginning of a study, date of birth (calendar timescale) or on time between two consecutive events (gap timescale). The calendar timescale is often preferred in analyses of disease evolution, for instance, in a randomized clinical trial. In this approach, the time of entering the interval for the risk of experiencing the jth event is equal to the time of the (j -1)th event, an individual can be at risk of having the jth event only after the time of the (j -1)th event (delayed entry). On the other hand, in the gap timescale, after every recurrent event, the time of origin for the consecutive event is renewed to 0. In the medical context, this timescale is less natural than the calendar timescale, but it might be considered if the number of events is low or if the occurrence of the events does not importantly affect the subject's condition.

For the recurrent part of the model, the individual contribution from the recurrent event process of type l (l = 1, 2) is given by the contribution from the right-censored observations and the event times and in the calendar timescale it is:

f T (l) i |u i (T (l) i , δ (l) i |u i ; ξ) = n (l) i j=1 r (l) ij (T (l) ij |u i ; ξ) δ (l) ij exp   - T (l) ij T (l) i(j-1) r (l) ij (t|u i ; ξ)dt   .
For the gap timescale the lower limit of the integral is 0 and the upper limit is the gap between the time of the (j -1)th event and the jth event [START_REF] Duchateau | Evolution of Recurrent Asthma Event Rate over Time[END_REF]. In case of only one type of recurrent events, r

(l) ij is obviously r ij .
Similarly, the individual contribution from the terminal event process is:

f T i |u i (T i , δ i |u i ; ξ) = (λ i (T i |u i ; ξ)) δ i exp - T i 0 λ i (t|u i ; ξ)dt .
For all the analyzed joint models the marginal likelihood (Equation 6) does not have an analytic form and integration is performed using quadrature methods. If a model includes one random effect that follows the Gamma distribution, the Gauss-Laguerre quadrature is used for the integral. The integrals over normally distributed random effects can be approximated using the Gauss-Hermite quadrature.

For the multivariate joint frailty model (3), the bivariate joint model for longitudinal data and a terminal event (4) and the trivariate joint model ( 5) it is required (except for the bivariate joint model with a random intercept only) to approximate multidimensional integrals. In this case, the standard non-adaptive Gauss-Hermite quadrature, that uses a specific number of points, gives accurate results but often can be time consuming and thus alternatives have been proposed. The multivariate non-adaptive procedure using fully symmetric interpolatory integration rules proposed by [START_REF] Genz | Fully Symmetric Interpolatory Rules for Multiple Integrals over Infinite Regions with Gaussian Weight[END_REF] offers advantageous computational time but in case of datasets in which some individuals have few repeated measurements, this method may be moderately unstable. Another possibility is the pseudo-adaptive Gauss-Hermite quadrature that 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 [START_REF] Rizopoulos | Joint Models for Longitudinal and Time-to-Event Data: With Applications in R[END_REF].

Parametric approach

We propose to estimate the hazard functions using cubic M -splines, piecewise constant functions (PCF) or the Weibull distribution. In case of PCF and the Weibull distribution applied to the baseline hazard functions, the full likelihood is used for the estimation. Using the PCF for the baseline hazard function h 0 (t) (of recurrent events or a terminal event), an interval [0, τ ] (with τ the last observed time among N individuals) is divided into n int subintervals in one of two manners: equidistant so that all the subintervals are of the same length or percentile so that in each subinterval the same number of events is observed. Therefore, the baseline hazard function is expressed by

h 0 (t) = n int i=1 I {t∈(t i-1 ,t i )} c i , with c i ≥ 0.
In the approach of PCF the crucial point is to choose the appropriate number of the intervals so that the estimated hazard function could capture enough the flexibility of the true function. The other possibility in the parametric approach is to assume that the baseline hazard function h 0 (t) comes from the Weibull distribution with the shape parameter a > 0 and the scale parameter b > 0. Then the baseline hazard function is defined by h 0 (t) = (at a-1 )/b a . This approach is convenient given the small number of parameters to estimate (only two for each hazard function) but the resulting estimated functions are monotone and this constraint, in some cases, might be too limiting.

Semi-parametric approach based using penalized likelihood

In the semi-parametric approach, with regard to expected smooth baseline hazard functions, the likelihood of the model is penalized by terms depending on the roughness of the functions [START_REF] Joly | A Penalized Likelihood Approach for Arbitrarily Censored and Truncated Data: Application to Age-Specific Incidence of Dementia[END_REF]. Cubic M -splines, polynomial functions of 3rd order, are combined linearly to approximate the baseline hazard functions [START_REF] Ramsay | Monotone Regression Splines in Action[END_REF].

For the estimated parameters ξ the full log-likelihood, l joint (ξ) = N i=1 lnL i (Θ i , ξ) is penalized in the following way:

pl(ξ) = l joint (ξ) -κ 1 ∞ 0 r 0 (t) 2 dt -κ 2 ∞ 0 λ 0 (t) 2 dt,
for model (1), ( 2), ( 5)

pl(ξ) = l joint (ξ) -κ 1 ∞ 0 r (1) 0 (t) 2 dt -κ 2 ∞ 0 r (2) 0 (t) 2 dt -κ 3 ∞ 0 λ 0 (t) 2 dt, for model (3) pl(ξ) = l joint (ξ) -κ 1 ∞ 0 λ 0 (t) 2 dt,
for model ( 4)

.

The positive smoothing parameters (κ 1 , κ 2 and κ 3 ) provide a balance between the data fit and the smoothness of the functions. Both the full and penalized log-likelihood are maximized using the robust Marquardt algorithm [START_REF] Marquardt | An Algorithm for Least-Squares Estimation of Nonlinear Parameters[END_REF], a mixture between the Newton-Raphson and the steepest descent algorithm.

Goodness-of-fit

For verification of model assumptions, residuals are used as a standard statistical tool for visual assessment. In the context of the survival data (recurrent and terminal events) these are the martingale residuals and in the context of the longitudinal data, both, the residuals conditional on random effects and the marginal residuals are often used.

The martingale residuals model whether the number of observed events is correctly predicted by a model. Principally, it is based on the counting processes theory and for subject i and time t they are defined as the difference between the number of events of subject i until t and the Breslow estimator of the cumulative hazard function of t. Let N i (t) be the counting process of the event of type p (recurrent or terminal), u i represent the random effects and the process's intensity

ζ (p) i (t|u i ) = u i ζ (p) 0 (t) exp(X pi (t) β p ) = u i ζ (p) i (t) with ζ (p)
0 (t) the baseline risk and X pi (t) the prognostic factors. The martingale residuals can be expressed by (for details, see [START_REF] Commenges | Standardized Martingale Residuals Applied to Grouped Left Truncated Observations of Dementia Cases[END_REF]:

M i (t) = N i (t) -u i t 0 W i (s) ζ (p) i (s)ds, (7)
where W i (t) is equal to 1 if the individual is at risk of the event at time t and 0 otherwise. The martingale residuals are calculated at T i , that is at the end of the follow-up (terminal event).

The assessment of the model is performed visually, the mean of the martingale residuals at a given time point should be equal to 0.

For the longitudinal data estimated in the framework of linear mixed effects models, there exist marginal residuals averaged on the population level defined by R (m) i

= y i -X Li βL and conditional residuals that are subject-specific, R (c) i = y i -X Li βL -Z i bi . These raw residuals are recommended for checking homoscedasticity of the conditional and marginal variance. For verification of the normality assumption and detection of outlying observations, the Cholesky residuals are more adapted as they represent decorrelated residuals and are defined by:

R (m) * i = U (m) i R (m) i , R (c) * i = U (c) i R (c) i ,
where the raw residuals are multiplied by the upper-triangular matrices ( U

(m) i and U (c) i )
obtained by the Cholesky decomposition of the variance-covariance matrices:

V R (m) i = V i -X Li ( N i=1 X Li V i -1 X Li ) -1 X Li = U (m) i U (m) i , V R (c) i = σ 2 I n i V i -1 V R (m) i V i -1 σ 2 I n i = U (c) i U (c) i ,
where V i is the marginal variance-covariance matrix of the longitudinal outcome y i , equal to Z i B 1 Z i + σ 2 I n i and of dimension l i × l i . The marginal Cholesky residuals should be approximately normally distributed and thus, their normal Q-Q plot allows to verify the assumption of normality.

Both, for the calculation of the martingale residuals and the residuals of the longitudinal outcome, the values of the random effects are necessary. For this purpose, the empirical Bayes estimators of u i are calculated using the formula for the posterior probability function:

f (u i |Θ i ; ξ) = f (Θ i |u i ; ξ)f (u i ; ξ) f (Θ i ; ξ) ∝ f (Θ i |u i ; ξ)f (u i ; ξ).
For the joint models, this expression does not have an analytical solution and the numerical computation is applied that finds the u i that maximizes f (u i |Θ i ; ξ):

u i = arg max u i f (u i |Θ i ; ξ),
and this is obtained using the Marquardt algorithm.

Setting 1

Exactly 2 recurrent events (x) before t

Complete information

Setting 2

At least 2 recurrent events (x) before t

Incomplete information

Setting 3

Whatever the history of recurrent events before t Illustration with two recurrent events.

Prediction of risk of the terminal event

Specific predictions can be obtained in the framework of joint modeling. Prediction consists of estimating the probability of an event at a given time t + w knowing the available information at prediction time t. Using a joint model, it is possible to estimate the probability of having the terminal event at time t+w given the history of the individual (occurrences of the recurrent events and/or the measurements of the longitudinal biomarker) prior to t. For the joint model for recurrent events and a terminal event (1), three settings of prediction were developed by [START_REF] Mauguen | Dynamic Prediction of Risk of Death Using History of Cancer Recurrences in Joint Frailty Models[END_REF]. In the first one, all the available information is accounted for, and we consider that this information is complete. In the second setting, all the known information is accounted for, however we consider that this information may be incomplete. Finally, in the third setting, recurrence information is not accounted for and only covariates are considered.

All the settings are represented in Figure 1. Here, we focus on the first setting, for which we present the predictions and the measures of predictive accuracy but in the package all the three settings are implemented for the joint models for recurrent and terminal events.

For the joint models with a longitudinal outcome a complete history of the biomarker is considered [START_REF] Król | 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[END_REF]. Thus, for model (4) the individual's history is the whole observed trajectory of the biomarker and for model (5) it is the the whole observed trajectory of the biomarker and all the observed occurrences of the recurrent event (see Figure 2).

The proposed prediction can be performed for patients from the population used to develop the model, but also for "new patients" from other populations. This is possible as the probabilities calculated are marginal, i.e., the conditional probabilities are integrated over the distribution of the random effects. Thus, values of patients' frailties are not required to estimate the probabilities of the event. However, it should be noted that the predictions include individual deviations via the estimated parameters of the random effects' distribution in a joint model.

We denote by t the time at which the prediction is made and by w the window of prediction. Thus, we are interested in the probability of the event (death) at time t + w, knowing what happened before time t. The general formulation of predicted probability of the terminal event conditional on random effects and patient's history is: 4) and the bottom setting to the model ( 5).

P i (t, t + w; ξ) = P(T * i ≤ t + w|T * i > t, F i (t), X i ; ξ), ( 8
)
where X i are all the covariates included in the model, F i (t) corresponds to the complete repeated data of patient i observed until time t. We define the complete history of recurrences H

J(l) i (t) = {N R(l) i (t) = J (l) , T (l) * i1 < • • • < T (l) * iJ ≤ t} (N R i (t)
is the counting process of the recurrent events, l = 1, 2 for the two types of recurrent events, and in case of only one type of recurrent events in the models the index (l) is omitted) and the history of the biomarker Y i (t) = {y i (t iK ), t i1 < • • • < t iK < t}. Therefore, the individual's history F i (t) depends on the model considered and is equal to:

F i (t) = H J i (t), for model (1), (2) F i (t) = {H J(1) i (t), H J(2) i (t)}, for model (3) F i (t) = Y i (t), for model (4) F i (t) = {H J i (t), Y i (t)}, for model (5)
and for the recurrent events we assume T * i0 = 0 and T * i(J+1) > t. For the estimated probabilities, confidence intervals are obtained by the Monte Carlo (MC) method, using the 2.5th and 97.5th percentiles of the obtained distribution (percentile confidence interval).

Brier score

In order to validate the prediction ability of a given model, a prediction error is proposed using the weighted Brier score [START_REF] Gerds | Consistent Estimation of the Expected Brier Score in General Survival Models with Right-Censored Event Times[END_REF][START_REF] Mauguen | Dynamic Prediction of Risk of Death Using History of Cancer Recurrences in Joint Frailty Models[END_REF]. It measures the distance between the prediction (probability of event) and the actual outcome (dead or alive). Inverse probability weighting corrects the score for the presence of right censoring. At a given horizon of prediction t + w, the error of prediction is calculated by:

Err t+w = 1 N t Nt i=1 [I(T * i > t + w) -(1 -P i (t, t + w; ξ))] 2 w i (t + w, G N ), (9)
where the weights w i (t + w, G n ) are defined by:

w i (t + w, G n ) = I(T * i ≤ t + w)δ i G N (T * i )/ G N (t) + I(T * i > t + w) G N (t + w)/ G N (t) ,
and G N (t) is the Kaplan-Meier estimate of the population censoring distribution at time t, N t is the number of patients at risk of the event (alive and uncensored).

Direct calculations of the Brier score are not implemented in the package frailtypack but using predictions P i (t, t + w; ξ)), the Brier score can be obtained using weights w i (•) function pec from the package pec (Gerds 2016) (for details, see the colorectal example in Section 5.3).

For an internal validation (on a training dataset, i.e., used for estimation) of the model as a prediction tool for new patients, a k-fold cross-validation is used to correct for over-optimism.

In this procedure, the joint model estimations are performed k times on k -1 partitions from the random split and the predictions are calculated on the left partitions. Finally, the Brier score can be calculated using predictions for all individuals at several time points to obtain cross-validated prediction error curves.

EPOCE

Another method to evaluate a model's predictive accuracy is the EPOCE estimator that is derived using prognostic conditional log-likelihood [START_REF] Commenges | Choice of Prognostic Estimators in Joint Models by Estimating Differences of Expected Conditional Kullback-Leibler Risks[END_REF]). This measure is adapted both for external data and then the mean prognosis observed loss (MPOL) is computed, as well as for the training data using the approximated cross-validated prognosis observed loss (CVPOL a ).

This measure of predictive accuracy is the risk function of an estimator of the true prognostic density function f * T |F (t),T * ≥t , where F(t) denotes the history of repeated measurements and/or recurrent events until time t. Using information theory this risk is defined as the expectation of the loss function, that is the estimator derived from the joint model f T |F (t),T * ≥t conditioned on T * > t and this can be written as EPOCE(t) = E(-ln(f T |F (t),T * >t )|T * > t). In case where the model is evaluated on the training data the approximated leave-one-out CVPOL a is defined by:

CVPOL a (t) = - 1 N t Nt i=1 F i ( ξi , t) + N Trace(H -1 K t ), (10)
where H -1 is the inverted Hessian matrix of the joint log-likelihood (Equation 6),

K t = 1 Nt(N -1) N i=1 I {T i >t} vi d i with vi = ∂F i (ξ,t) ∂ξ | ξ and di = ∂l i (Θ,ξ) ∂ξ | ξ.
The individual contribution to the log-likelihood of a terminal event at t defined for the individuals that are still at risk of the event at t can be written as:

F i ( ξ, t) = ln        u i f Y i (t)|u i (Y i (t)|u i ; ξ) γ L l k=1 f H J(k) i (t)|u i (H J(k) i (t)|u i ; ξ) γ (k) R f T i |u i (T i , δ i |u i ; ξ)f u i (u i ) du i u i f Y i (t)|u i (Y i (t)|u i ; ξ) γ L l k=1 f H J(k) i (t)|u i (H J(k) i (t)|u i ; ξ) γ (k) R S t i (T i |u i ; ξ)f u i (u i ) du i       
S t i is the survival function of the terminal event for individual i. Again, l denotes the number of types of recurrent events included in the model.

In case of external data MPOL(t) is expressed by the first component of the sum in CVPOL a (t) (Equation 10). Indeed, the second component corresponds to the statistical risk introduced to CVPOL a in order to correct for over-optimism coming from the use of the approximated crossvalidation. The comparison of the models using EPOCE can be done visually by tracking the 95% intervals of difference in EPOCE.

The EPOCE estimator has some advantages compared to the prediction error, in particular for the training data the approximated cross-validation technique is less computationally demanding than the crude leave-one-out cross-validation method for the prediction error, which is important in case of joint models when time of calculation is often long [START_REF] Proust-Lima | Joint Latent Class Models for Longitudinal and Time-to-Event Data: A Review[END_REF].

Modeling and prediction using the R package frailtypack

Estimation of joint models

In the package frailtypack there are four different functions for the estimation of joint models, one for each model. The joint models for recurrent events and a terminal event (1) (and models for two clustered survival outcomes) as well as the general joint frailty models (2) are estimated with function frailtyPenal. The joint models for two recurrent events and a terminal event (3) are estimated with multivPenal. The estimation of joint models for a longitudinal outcome and a terminal event (4) is performed with longiPenal. Finally, function trivPenal estimates trivariate joint models (5). All these functions make calls to compiled Fortran codes programmed for computation and optimization of the log-likelihood. In the following, we detail each of the joint model functions.

frailtyPenal function frailtyPenal(formula, formula.terminalEvent, data, recurrentAG = FALSE, cross.validation = FALSE, jointGeneral, n.knots, kappa, maxit = 300, hazard = "Splines", nb.int, RandDist = "Gamma", betaknots = 1, betaorder = 3, initialize = TRUE, init.B, init.Theta, init.Alpha, Alpha, init.Ksi, Ksi, init.Eta, LIMparam = 1e-3, LIMlogl = 1e-3, LIMderiv = 1e-3, print.times = TRUE)

Argument formula is a two-sided formula for a survival object 'Surv' from the survival package [START_REF] Therneau | survival: Survival Analysis. R package version 2[END_REF]) and it represents the recurrent event process (the first survival outcome for the joint models in case of clustered data) with the combination of covariates on the right-hand side, the indication of a grouping variable (with term cluster(group)) and the indication of the variable for the terminal event (e.g., terminal(death)). It should be noted that the function cluster(x) is different from that included in the package survival. In both cases it is used for the identification of the correlated groups but in frailtypack it indicates the application of the frailty model and in survival, a GEE (generalized estimating equations) approach is used, without random effects. Argument formula.terminalEvent requires the combination of covariates related to the terminal event on the right-hand side. The name of the data.frame containing the variables used in the function should be specified in the argument data. The logical argument recurrentAG indicates whether the calendar timescale for recurrent events or clustered data with the counting process approach of [START_REF] Andersen | Cox's Regression Model for Counting Processes: A Large Sample Study[END_REF] (TRUE) or the gap timescale (FALSE; by default) is to be used. The argument for the cross-validation cross.validation is not yet implemented for the joint models; thus its logical value must be FALSE. The smoothing parameters κ for a joint model can be chosen by first fitting suitable shared frailty and Cox models with the cross-validation method.

The general joint frailty models (2) can be estimated if argument jointGeneral is TRUE. In this case, the additional gamma frailty term is assumed and the α parameter is not considered. These models can be applied only with the Gamma distribution for the random effects.

The type of approximation of the baseline hazard functions is defined by argument hazard and can be chosen among "Splines" for semiparametric functions using equidistant intervals, "Splines-per" using percentile intervals, "Piecewise-equi" and "Piecewise-per" for piecewise constant functions using equidistant and percentile intervals, respectively and "Weibull" for the parametric Weibull baseline hazard functions. If either "Splines" or "Splines-per" is chosen for the baseline hazard functions, arguments kappa for the positive smoothing parameters and n.knots should be given with the number of knots chosen between 4 and 20 which corresponds to n.knots+2 splines functions for approximation of the baseline hazard functions (the same number for hazard functions for both outcomes). If "Percentile-equi" or "Percentile-per" is chosen for the approximation, argument nb.int should be given with a 2-element vector of numbers of time intervals (1-20) for the two baseline hazard functions of the model.

Argument RandDist represents the type of the random effect distribution, either "Gamma" for the Gamma distribution or "LogN" for the normal distribution (log-normal joint model). If it is assumed that α in model ( 1) is equal to zero, argument Alpha should be set to "None".

In case of time dependent covariates, arguments betaknots and betaorder are used for the number of inner knots used for the estimation of B-splines (1, by default) and the order of B-splines (3 for quadratic B-splines, by default), respectively.

The rest of the arguments are allocated for the optimization algorithm. Argument maxit declares the maximum number of iterations for the Marquardt algorithm. For a joint nested frailty model, i.e., a model that allows joint analysis of recurrent and terminal events for hierarchically clustered data, argument initialize determines whether the parameters should be initialized with estimated values from the appropriate nested frailty models. Arguments init.B, init.Theta, init.Eta and init.Alpha are vectors of initial values for regression coefficients, variances of the random effects and for the α parameter (by default, 0.5 is set for all the parameters). Arguments init.Ksi and Ksi are defined for joint nested frailty models and correspond to initial values for the flexibility parameter and the logical value indicating whether to include this parameter in the model or not, respectively. The convergence thresholds of the Marquardt algorithm are for the difference between two consecutive log-likelihoods (LIMlogl), for the difference between the consecutive values of estimated coefficients (LIMparam) and for the small gradient of the log-likelihood (LIMderiv). All these threshold values are 10 -3 by default. Finally, argument print.times indicates whether to print the iteration process (the information note about the calculation process and time taken by the program after terminating), the default is TRUE.

The function frailtyPenal returns objects of class 'jointPenal' if joint models are estimated. It should be noted that using this function for univariate models, shared frailty models [START_REF] Rondeau | frailtypack: An R Package for the Analysis of Correlated Survival Data with Frailty Models Using Penalized Likelihood Estimation or Parametrical Estimation[END_REF] and Cox models can be applied as well and they result in objects of class 'frailtyPenal'. For both classes methods print(), summary() and plot() are available.

multivPenal function multivPenal(formula, formula.Event2, formula.terminalEvent, data, initialize = TRUE, recurrentAG = FALSE, n.knots, kappa, maxit = 350, hazard = "Splines", nb.int, print.times = TRUE)

This function allows to fit the multivariate frailty models (3). Argument formula must be a two-sided formula for a 'Surv' object, corresponding to the first type of the recurrent event (no interval-censoring is allowed). Arguments formula.Event2 refer to the second type of the recurrent event and formula.terminalEvent to the terminal event, and are equal to linear combinations related to the respective events. The rest of the arguments is analogical to frailtyPenal. Arguments n.knots (values between 4 and 20) and kappa must be vectors of length 3 for each type of event, first for the recurrent event of type 1, second for the recurrent event of type 2 and third for the terminal event. The function multivPenal return objects of class 'multivPenal' with print(), summary() and plot() methods available.

longiPenal function longiPenal(formula, formula.LongitudinalData, data, data.Longi, random, id, intercept = TRUE, link = "Random-effects", left.censoring = FALSE, n.knots, kappa, maxit = 350, hazard = "Splines", nb.int, init.B, init.Random, init.Eta, method.GH = "Standard", n.nodes, LIMparam = 1e-3, LIMlogl = 1e-3, LIMderiv = 1e-3, print.times = TRUE)

In this function for the joint analysis of a terminal event and a longitudinal outcome, argument formula refers to the terminal event and the left-hand side of the formula must be a 'Surv' object and the right-hand side indicate the explanatory variables. Argument formula.LongitudinalData is equal to the sum of fixed explanatory variables for the longitudinal outcome. For the model, two datasets are required: data containing information on the terminal event process and data.Longi with data related to longitudinal measurements. Random effects associated to the longitudinal outcome are defined with random using the appropriate names of the variables from data.Longi. If a random intercept is assumed, the option "1" should be used. For a random intercept and slope, arguments random should be equal to a vector with elements "1" and the name of a variable representing time points of the biomarker measurements. At the moment, more complicated structures of the random effects are not available in the package. The name of the variable representing the individuals in data.Longi is indicated by id. The logical argument intercept determines whether a fixed intercept should be included in the longitudinal part or not (default is TRUE). Two types of subject-specific link function can be selected and are defined with the argument link.

The default option "Random-effects" represents the link function of the random effects b i , otherwise the option is "Current-level" for the link function of the current level of the longitudinal outcome m i (t).

The initial values of the estimated parameters can be indicated by init.B for the fixed covariates (a vector of values starting with the parameters associated to the terminal event and then for the longitudinal outcome, interactions in the end of each component), init.Random for the vector of elements of the Cholesky decomposition of the covariance matrix of the random effects u i and init.Eta for the regression coefficients associated to the link function.

There are three methods of the Gaussian quadrature to approximate the integrals to choose from. The default "Standard" corresponds to the non-adaptive Gauss-Hermite quadrature for multidimensional integrals. The other possibility is the pseudo-adaptive Gaussian quadrature ("Pseudo-adaptive") [START_REF] Rizopoulos | Joint Models for Longitudinal and Time-to-Event Data: With Applications in R[END_REF]. Finally, the multivariate non-adaptive Gaussian quadrature using the algorithm implemented in a Fortran subroutine HRMSYM is indicated by "HRMSYM" [START_REF] Genz | Fully Symmetric Interpolatory Rules for Multiple Integrals over Infinite Regions with Gaussian Weight[END_REF]. The number of the quadrature nodes (n.nodes) can be chosen among 5, 7, 9, 12, 15, 20 and 32 using this argument (default is 9).

The function longiPenal returns objects of class 'longiPenal' with print(), summary() and plot() methods available.

trivPenal function 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 = 300, hazard = "Splines", nb.int, 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)

The function for the trivariate joint model comprises three formulas for each type of process. The first two arguments are analogous to frailtyPenal, argument formula, referring to recurrent events, is a two-sided formula for a 'Surv' object on the left-hand side and covariates on the right-hand side (with indication of the variable for the terminal event using method terminal) and argument formula.terminalEvent represents the terminal event and is equal to a linear combination of the explanatory variables. Finally, argument formula.LongitudinalData as in function longiPenal corresponds to the longitudinal outcome indicating the fixed effect covariates. The rest of the arguments are detailed in the descriptions of functions frailtyPenal and longiPenal. The function trivPenal returns objects of class 'trivPenal' with print(), summary() and plot() methods available.

Prediction

The current increase of interest in the joint modeling of correlated data is often related to the individual predictions that these models offer. Indeed, calculating the probabilities of a terminal event given a joint model results in precise predictions that consider the past of an individual. Moreover, there exist statistical tools that evaluate a model's capacity for these dynamic predictions. In the package frailtypack we provide the prediction function for dynamic predictions of a terminal event in a finite horizon, the epoce function for evaluating predictive accuracy of a joint model and the Diffepoce function for comparing the accuracy of two joint models.

Predicted probabilities with prediction function

In the package it is possible to calculate the prediction probabilities for the Cox proportional hazard, shared frailty (for clustered data, [START_REF] Rondeau | frailtypack: An R Package for the Analysis of Correlated Survival Data with Frailty Models Using Penalized Likelihood Estimation or Parametrical Estimation[END_REF]) and joint models. Among the joint models the predictions are provided for the standard joint frailty models (recurrent events and a terminal event), for the joint models for a longitudinal outcome and a terminal event and for the trivariate joint models (a longitudinal outcome, recurrent events and a terminal event). These probabilities can be calculated for a given prediction time and window or for a vector of times, with varying prediction time or varying window. For the shared frailty models for clustered events, marginal and conditional on a specific cluster predictions can be calculated and for the joint models only the marginal predictions are provided. Finally, for the joint frailty models the predictions are calculated in three settings: given the exact history of recurrences, given the partial history of recurrences (the first J recurrences) and ignoring the past recurrences. For the joint models with a longitudinal outcome (bivariate and trivariate) only the predictions considering the patient's complete history are provided. For all the aforementioned predictions the following function is proposed:

prediction(fit, data, data.Longi, t, window, group, MC.sample = 0)

Argument fit indicates the, a 'frailtyPenal', 'jointPenal', 'longiPenal' or 'trivPenal' object. The data with individual characteristics for predictions must be provided in the data frame data with information on the recurrent events and covariates related to recurrences and the terminal event and in case of 'longiPenal' and 'trivPenal' objects the data frame data.Longi needs to contain repeated measurements and covariates related to the longitudinal outcome. These two datasets must refer to the same individuals for which the predictions will be calculated. Moreover, the names and the types of variables should be equivalent to those in the dataset used for estimation. The details on how to prepare correctly the data are presented in appropriate examples (Section 5).

Argument t is a time or vector of times for predictions and window is a horizon or vector of horizons. The function calculates the probability of the terminal event between a time of prediction and a horizon (both arguments are scalars), between several times of prediction and a horizon (t is a vector and window a scalar) and finally, between a time of prediction and several horizons (t is a scalar and window a vector of positive values).

For all the predictions, confidence intervals can be calculated using the MC method with MC.sample number of samples (maximum 1000). If the confidence bands are not to be calculated argument MC.sample should be equal to 0 (the default).

Predictive accuracy measure with epoce and Diffepoce functions

Predictive ability of joint models can be evaluated with function epoce that computes the estimators of EPOCE, MPOL and CVPOL a . For a given estimation, the evaluation can be performed on the same data and then both MPOL and CVPOL a are calculated, as well as on a new dataset but then only MPOL is calculated as the correction for over-optimism is not necessary. The call of the function is:

epoce(fit, pred.times, newdata = NULL, newdata.Longi = NULL)

with fit an object of class 'jointPenal', 'longiPenal' or 'trivPenal', pred.times a vector of time for which the calculations are performed. In case of external validation, new datasets newdata and newdata.Longi should be provided (newdata.Longi only in case of 'longiPenal' and 'trivPenal' objects). However, the names and types of variables in newdata and newdata.Longi must be the same as in the datasets used for estimation.

To compare the predictive accuracy of two joint models fit to the same data but possibly with different covariates, the simple comparison of obtained values of EPOCE can be enhanced by calculating the 95% tracking interval of difference between the EPOCE estimators. For this purpose we propose the function:

Diffepoce(epoce1, epoce2)
where epoce1 and epoce2 are objects inheriting from the 'epoce' class.

Illustrating examples

The package frailtypack provides various functions for models for correlated outcomes and survival data. The Cox proportional hazard model, the shared frailty model for recurrent events (clustered data), the nested frailty model, the additive frailty model and the joint frailty model for recurrent events and a terminal event have already been illustrated elsewhere [START_REF] Rondeau | frailtypack: A Computer Program for the Analysis of Correlated Failure Time Data Using Penalized Likelihood Estimation[END_REF][START_REF] Rondeau | frailtypack: An R Package for the Analysis of Correlated Survival Data with Frailty Models Using Penalized Likelihood Estimation or Parametrical Estimation[END_REF].

In this section we focus on extended models for correlated data presented in Section 2 using three datasets: readmission, dataMultiv and colorectal. Although the joint frailty model has already been presented in the literature, we illustrate its usage as the form of the function has developed in the meantime.

Example on dataset readmission for joint frailty models

Dataset readmission comes from a rehospitalization study of patients after a surgery and diagnosed with colorectal cancer [START_REF] Gonzalez | Sex Differences in Hospital Readmission among Colorectal Cancer Patients[END_REF][START_REF] Rondeau | frailtypack: An R Package for the Analysis of Correlated Survival Data with Frailty Models Using Penalized Likelihood Estimation or Parametrical Estimation[END_REF]. It contains information on times (in days) of successive hospitalizations and death (or last registered time of follow-up for right-censored patients) counting from date of surgery, and patient characteristics: type of treatment, sex, Dukes' tumoral stage, comorbidity Charlson's index and survival status. The dataset includes 403 patients with 861 rehospitalization events in total. Among the patients 112 (28%) died during the study.

Standard joint frailty model

We adapt the joint model for recurrent events and a terminal event using the gap timescale from the example given in [START_REF] Rondeau | frailtypack: An R Package for the Analysis of Correlated Survival Data with Frailty Models Using Penalized Likelihood Estimation or Parametrical Estimation[END_REF]. The model modJoint.gap is defined:

R> library("frailtypack") R> data("readmission", package = "frailtypack") R> modJoint.gap <-frailtyPenal(Surv(time, event) ~cluster(id) + dukes + + charlson + sex + chemo + terminal(death), + formula.terminalEvent = ~dukes + charlson + sex + chemo, + data = readmission, n.knots = 8, kappa = c(2.11e+08, 9.53e+11))

This model includes Dukes' stage, Charlson's index, sex and treatment as covariates for both hospitalizations and death. The frailties are from the Gamma distribution (default option) and the baseline hazard functions are approximated by splines with 8 knots and the smoothing parameters for the penalized log-likelihood are 2.11e + 8 for the recurrent part and 9.53e + 11 for the terminal part. To find the optimal number of knots, we first fitted the model with a small number of knots (n.knots = 4) and increased the number of knots until the graph of the baseline hazard functions was not changing importantly anymore. The smoothing parameters are obtained from a shared frailty and Cox model with respectively recurrent and terminal event as the outcome using the cross-validation method.

With this model, it is found that chemotherapy is a prognostic factor only on death with a positive association (HR = 2.99, p < 0.001). Both Charlson's index (Index ≥ 3 vs. Index 0) and Dukes' stage (Stage C and D vs. Stages A, B) are positively related to the recurrent and terminal events. A detailed description of the output for the standard joint frailty models are presented in [START_REF] Rondeau | frailtypack: An R Package for the Analysis of Correlated Survival Data with Frailty Models Using Penalized Likelihood Estimation or Parametrical Estimation[END_REF].

To verify whether the model predicts correctly the number of observed events, we represent the martingale residuals for both events against the follow-up time. These residuals in a well adjusted model should have a mean equal to 0 and thus a smoothing curve added to a graph should be approximately overlapping with the horizontal line y = 0. The following code produces the plots given in Figure 3: For the rehospitalization process the mean of residuals is approximately 0 with the smooth curve close to the line y = 0, but in case of death this tendency is deviated by relatively higher values for short follow-up times. This may suggest, that the model may have underestimated the number of deaths in the early follow-up period. The identified individuals of which the residuals result in non-zero mean, have short intervals between their rehospitalization and death (1 day). Indeed, the removal of these individuals (50 patients) results in residuals with the mean close to 0 all along the follow-up period (plot not shown here).

The package frailtypack provides also the estimation of the random effects. The vector frailty.pred from a 'jointPenal' object contains the individual empirical Bayes estimates. They can be graphically represented for each individual with additional information on number of events (point size) to identify the outlying data.

R> plot(1:403, modJoint.gap$frailty.pred, xlab = "Id of patients", + ylab = "Frailty predictions for each patient", type = "p", axes = FALSE, + cex = as.vector(table(readmission$id)), pch = 1, ylim = c(-0.1, 7), + xlim = c(-2, 420)) R> axis(1, round(seq(0, 403, length = 10), digit = 0)) R> axis(2, round(seq(0, 7, length = 10), digit = 1))

Figure 4 shows the values of frailty prediction for each patient with an association to the number of events (the bigger the point, the greater the number of rehospitalizations). The frailties tend to have bigger values if the number of events of a given individual is high.

From the plot it can be noticed that there is an outlying frailty suggesting verification of the follow-up of the concerned individual.

Time-varying coefficients

In the framework of standard joint frailty models, it is possible to fit the models with timevarying effects of prognostic factors. Using function timedep in a formula of frailtyPenal, the time-dependent coefficients can be estimated using B-splines of order q (option betaorder) with m interior knots (option betaknots). In the example of the readmission dataset we are interested in verifying whether the variable sex has a time-varying effect on both recurrent and terminal events. Thus, we fit a model equivalent to modJoint.gap but with time dependent effects assuming quadratic B-splines and 3 interior knots: In the result, using the method print the estimated values of parameters with time-constant effects and graphics of log-hazard ratios for time-dependent variables for each event are obtained (Appendix B, Figure 10). For rehospitalizations, we found firstly a protective effect q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Id of patients for females (β(t) < 0) and later an increased risk (β(t) > 0). For death, at the beginning, the effect of sex was weakening the risk but shortly became non-influential (β(t) around 0).

The PH assumption for the variable sex can be checked using the LR test. We compare two models: modJoint.gap (related to the null hypothesis that the effect is constant in time:

H 0 : β R (t) = β R , β T (t) = β T )
and modJoint.gap.timedep (related to the alternative hypothesis of time-varying effects:

H 1 : β R (t) = β R , β T (t) = β T ):
R> LR.statistic <--2 * modJoint.gap$logLik + + 2 * modJoint.gap.timedep$logLik R> p.value <-signif(1 -pchisq(LR.statistic, df = 10), 5) Given the obtained p value = 0.049, the PH assumption for the variable sex is not satisfied (at the level 0.05). Next, we check whether sex with time-dependent effects is an influential prognostic factor. Thus, again, we use the LR test to compare two models: modelJoint.gap.nosex without the covariate sex (model related to the null hypothesis: H 0 : β R (t) = 0, β T (t) = 0) and modelJoint.gap.timedep (related to the alternative hypothesis: The test shows that the variable sex has a significant time-varying effect (p value < 0.001). In this example, the variable sex is found significant for both, the PH model and the non-PH model, but we showed that this variable does not satisfy the PH assumption.

H 1 : β R (t) = 0, β T (t) = 0): R> modJoint.

Dynamic predictions

Using the joint model modJoint.gap we can calculate predicted probabilities of death using prediction. These predictions may serve as a tool to compare two exemplars of patients with different histories of recurrences but the same values of the prognostic factors to study the effect of the events on the survival [START_REF] Mauguen | Dynamic Prediction of Risk of Death Using History of Cancer Recurrences in Joint Frailty Models[END_REF]. On the other hand, the influence of some explanatory covariates can be examined for patients that are considered to have the same history of recurrences. Here, we aim at evaluating the predictive effect of the Dukes' stage on survival considering the history of hospitalizations. We compare the predicted risk of death for two patients having the same characteristics (men with Charlson's index 0 and the chemotherapy treatment) and having two hospitalizations: 1 and 1.5 year after their surgeries. We set the time of prediction to 2 years and calculate the probability in a time window of 3 years (we apply a moving window with a step of 0.5 years). Patient 1 has Dukes' stage A and patient 2, Dukes' stage D. We focus on the predicted probabilities regarding the complete history of recurrences (Equation 8) and compare the results with those obtained using the incomplete history and ignoring the history.

To prepare data for the predictions, we start with an empty data frame with the variables of interest and the covariates: R> datapred <-data.frame(time = 0, event = 0, id = 0, dukes = 0, + charlson = 0, sex = 0, chemo = 0) R> datapred[, 4:7] <-lapply(datapred[, 4:7], as.factor) R> levels(datapred$dukes) <-c(1, 2, 3) R> levels(datapred$charlson) <-c(1, 2, 3) R> levels(datapred$sex) <-c(1, 2) R> levels(datapred$chemo) <-c(1, 2) Patient 1 with Dukes' stage A has two observed hospitalizations at the 365th and 548th day after the surgery:

R> datapred[1, ] <-c(365, 1, 1, 1, 1, 1, 2) R> datapred[2, ] <-c(548, 1, 1, 1, 1, 1, 2)
Patient 2 has the hospitalizations observed in the same times as Patient 1 but is assumed to have Dukes' stage B:

R> datapred[3, ] <-c(365, 1, 2, 3, 1, 1, 2) R> datapred[4, ] <-c(548, 1, 2, 3, 1, 1, 2)
We calculate the predictions for both patients: In the result, three types of predictions are calculated for 6 time horizons. As it has been already observed from the estimates of the model, an increased Dukes' stage (D) is positively associated with death. Figure 5 compares the predicted probabilities of death in the three settings. The predictions using the exact number of recurrences (2 hospitalizations) and at least 2 recurrences are very close to each other and are higher than the risk obtained without the information on recurrent events. Indeed, the significant estimate of variance of the frailty θ and significant positive estimate of α imply the positive association between the process of recurrences and death. Finally, all the predicted probabilities are higher for the patient with Dukes' stage D compared to the patient with Dukes' stage A.

Joint frailty model for clustered data

The joint models for clustered survival data can be estimated using again the frailtyPenal function. A dataset should include information on two survival outcomes for individuals from several groups. This model is presented using the readmission dataset with artificially created clusters on individuals. The first survival event will be the first observed rehospitalization and the second event, death. The framework of semi-competing risks is used here, thus individuals' follow-up stops at time of the rehospitalization, death or in case when none of these events are observed, the censoring time. n.knots = 8, kappa = c(1.e+10, 1.e+10), recurrentAG = TRUE, + Alpha = "None")

In the result, the estimates of prognostic factors for both types of events are obtained. The estimate of the variance of the frailty term indicates whether, at the cluster level, the processes are associated with each other and measures the heterogeneity between individuals (intracluster correlation). In the given example the estimate of the variance θ is significantly different from 0 (p value = 0.037), thus there is a positive association between the risk of hospitalizations and death via the non-observed frailty.

In case of the joint frailty models for clustered data it should be noted that sufficient amount of information must be provided, i.e., number of observations per cluster. Otherwise, given the complexity of the model, convergence might not be attained. The parameter α is assumed to be equal to 1 as these models are defined in the framework of semi-competing risks and not of recurrent events.

General joint frailty model

To estimate the general frailty model, argument jointGeneral must be equal to TRUE in function frailtyPenal. We applied this model to the original readmission dataset assuming two independent frailty terms using the following code:

R> modJoint.general <-frailtyPenal(Surv(time, event) ~cluster(id) + + dukes + charlson + sex + chemo + terminal(death), + formula.terminalEvent = ~dukes + charlson + sex + chemo, + data = readmission, jointGeneral = TRUE, n.knots = 8, + kappa = c(2.11e+08, 9.53e+11))

In the output of the function, estimations of the variances of both random effects are given. For the analyzed example, the estimated variance θ of the frailty u i associating recurrent events and death indicates strong relationship between the processes ( θ = 0.68, p value < 0.001). Moreover, the estimate of η implies small but significant dependence between the recurrent event gap times explained by the frailty v i (η = 0.01, p value < 0.001). This information complements the inference from the standard joint frailty model because it separates the correlation linked to the recurrent events with correlation between the recurrent events and the terminal event.

Example on dataset dataMultiv for multivariate joint frailty model

For the following example we use a generated dataset for 800 individuals from model (3) with 2 types of recurrent events and a terminal event. The random effects are assumed to follow the normal distribution u i ∼ N (0, 0.5), v i ∼ N (0, 0.5) and the correlation coefficient ρ = 0.5. The coefficients for the random effects are α 1 = α 2 = 1. The baseline hazard functions r

(1) 0 (t), r

(2) 0 (t) and λ 0 (t) follow the Weibull distribution and the time for right censoring is fixed at 5. The generated data includes 1652 observations. For a detailed description of the generation scenario see [START_REF] Mazroui | General Joint Frailty Model for Recurrent Event Data with a Dependent Terminal Event: Application to Follicular Lymphoma Data[END_REF].

The dataset includes the individuals' times of events with variables indicating the type of event: INDICREC for the recurrent event of type 1 (local recurrences), INDICMETA for the recurrent event of type 2 (metastases) and INDICDEATH for the censoring status (death). Additionally there are 3 binary covariates v1, v2 and v3.

Multivariate frailty model

We consider the multivariate frailty model for the exemplary dataset dataMultiv to study jointly local recurrences, metastases and death for patients diagnosed with cancer. To define the model, three formulas must be defined in the function with additional indication on variables including status of the second recurrent event (event2) and of the terminal event (terminal), both included in the first formula. All the baseline hazard functions must be of the same type (Weibull, splines or piecewise constant). We fit the model as follows (computational time 54 minutes on a personal computer with an Intel Core i7 3.40 GHz processor and 8 GB RAM running Windows 7): R> data("dataMultiv", package = "frailtypack") R> modMultiv. ---------- The output presents the results for prognostic factor estimates for each type of event. The estimates of parameters associated with the random effects are given by the variance of the frailty related to the first type of the recurrent events and the association with the terminal event (theta1), the variance of the frailty related to the second type of the recurrent events and the association with the terminal event (theta2) and the correlation coefficient between the frailties (rho). The sign and strength of the dependency between the recurrent event of type 1 (2) and the terminal event is represented by alpha1 (alpha2). In the analyzed example, both theta1 and theta2 are not significantly different from 0, thus there are no dependencies between the processes explained by the non-observed factors.

Example on dataset colorectal for models with longitudinal data

Datasets colorectal and colorectal.Longi represent a random selection of 150 patients from the multi-center randomized phase III clinical trial FFCD 2000-05 of patients diagnosed with metastatic colorectal cancer not amenable to curative intent surgery. The trial was conducted between 2002 and 2007 in France by Fédèration Francophone de Cancérologie Digestive (FFCD; [START_REF] Ducreux | Sequential Versus Combination Chemotherapy for the Treatment of Advanced Colorectal Cancer (FFCD 2000-05): An Open-Label, Randomised, Phase 3 Trial[END_REF]. The data contains a follow-up of tumor size measurements (sum of the longest diameters of target lesions) and times of apparition of new lesions as recurrent events. Moreover, some baseline characteristics (age, WHO performance status and previous resection), treatment arm (combination vs. sequential) and time of death (or last observed time for a right-censored individual) are included in the data. Dataset colorectal provides information on the recurrent event and death and dataset colorectal.Longi on the measurements of tumor size. A total of 906 tumor size measurements and 289 of recurrences were recorded for the patients included. Among them, 121 died during the study.

The variable tumor.size in colorectalLongi is the transformed sum of the longest diameters (SLD * ) of an individual's target lesions measured during a visit (SLD * = (SLD 0.3 -1)/0.3). The status of new lesions occurrence is registered in new.lesions in dataset colorectal. In this dataset start of time interval time0 (0 or time of previous recurrence) and time of event time1 (recurrence or censoring by terminal event) represent information for times of apparition of new lesions and for death (or right censoring).

We provide extracts of colorectal and colorectalLongi datasets to guide the users how to prepare suitable datasets for joint models with longitudinal data using the package. These functions require datasets in long format (one row per observation and usually several rows per individuals) and only in case of data in longiPenal the long format is one row per individual (as it contains information on the terminal event only). In the example of the colorectal dataset this is represented as follows:

R> data("colorectalLongi", package = "frailtypack") R> head(colorectalLongi, 10) On average, the tumor size significantly decreases in time in interaction with the treatment, this effect is more important in the C arm (-0.63, p value = 0.001). However, there is no effect of treatment arm on the risk of death (HR = 0.91, p value = 0.65). The age of patients at baseline does not have any effect neither on the tumor size nor on the risk of death. The performance status WHO 2 evaluated before the treatment is a prognostic factor both for the tumor size (0.74, p value < 0.001) and overall survival (HR = 2.23, p value = 0.001).

id
It should be noted that the model is fitted on a subset of the original trial in which, using the data of all patients, some of the prognostic effects were found to be different [START_REF] Król | 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[END_REF].

The processes are linked together via the random intercept and slope of the longitudinal trajectory. This association is significant for the random intercept implying that with the increase of individual deviation from the population average tumor size, the risk of death increases as well (η t1 = 0.32, p value < 0.001).

To verify the goodness-of-fit of the model, the estimated baseline hazard function, martingale residuals for the terminal event and residuals for the longitudinal outcome can be plotted using the following code (results not shown here): Next, we explore predictiveness of this model comparing it with a model in which the link function is represented by the current level of the biomarker m i (t). Thus, we fit a bivariate model with the same covariates and characteristics but with the option link = "Current-level" (model modLongi2). Then, for times of prediction between 1.2 and 2.5 years we calculate the EPOCE estimator. We compare the models by plotting the values of CVPOL a and the tracking intervals for the differences between the models.

R> modLongi2 <-longiPenal(Surv(time0, 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, + n.knots = 8, kappa = 0.93, method.GH = "Pseudo-adaptive", n.nodes = 7) R> time <-seq(1.2, 2.4, 0.1) R> epoce <-epoce(modLongi, time) R> epoce2 <-epoce(modLongi2, time) R> diff <-Diffepoce(epoce, epoce2) R> plot(temps, epoce$cvpol, ylab = "CVPOL", xlab = "time", pch = 6, + col = "darkcyan", type= "b", ylim = c(0.4, 1.0)) R> points(temps, epoce2$cvpol, pch = 15, col = "brown3", type = "b", + lty = 2) R> legend(1.9, 1, legend = c("modLongi", "modLongi2"), pch = c(6, 15), + col = c("darkcyan", "brown3"), bty = "n", cex = 0.7, lty = c(1, 2)) R> plot(diff)

The results of CVPOL a are presented in Figure 6. Model modLongi with the random effects as the link function has better predictive abilities than model modLongi2 with the current level of the biomarker as the link function until 1.9 year. After this time point, the tendency inverses. However, these difference are significant only at 2.4 years of treatment. Using method AtRisk for 'epoce' objects, we may verify for how many subjects at risk, CVPOL a is calculated at this prediction time (epoce$AtRisk). The significant difference between the models is found for the time point for which the number of subjects considered is relatively small (25). Thus, we can conclude that the models are close to each other and the choice of the link function does not strongly influence the predictiveness.

The comparison of the predictive abilities of bivariate models can be useful for the choice of the appropriate trivariate model. If we are interested in the prediction of death it is important to choose a model that provides the best predictive abilities for survival. For this reason, we choose a model with such a structure of the link function that provides the best predictive abilities, for both bivariate and trivariate models. First, we fit appropriate bivariate models (without the recurrent part) with different link structures and compare their predictive abilities. Then, we fit the trivariate model with the same link structure as in the bivariate model with the best predictive abilities. For the colorectal dataset the differences in CVPOL a are not of great importance, thus we apply a trivariate model with the random effects of the biomarker as the link function, as this model is less computationally intensive.

Joint model for longitudinal data, recurrent events and a terminal event

The package frailtypack allows the estimation of the trivariate models for a longitudinal biomarker, recurrent events and a terminal event. In the dataset colorectal, the recurrent events are represented by the appearance of new lesions during the treatment. Usually in clinical trials the size of new lesions is not registered and thus, their burden cannot be simply added to the measure of the tumor size of the target lesions. However, it is of interest to add the information on new lesions to a model as it influences overall survival. A trivariate model can be a solution for such a goal and can be implemented using function trivPenal.

We fit a model with the calendar timescale for recurrent events and the baseline hazard functions approximated by splines with 8 knots. The smoothing parameter values are found from the separate models (shared frailty model and Cox model). A random intercept and a random slope and left-censoring (threshold s = -3.33) are assumed for the biomarker. The pseudo-adaptive Gaussian quadrature with 9 nodes is chosen for calculation of integrals.

Firstly, we find initial values for the covariates. For the longitudinal outcome and the terminal event we use the estimates from the bivariate model, modLongi. The covariates related to the recurrent events are initialized using the results of a shared frailty model for the appearance of new lesions:

R> modShared <-frailtyPenal(Surv(time0, time1, new.lesions) ~cluster(id) + + age + treatment + who.PS, data = colorectal, recurrentAG = TRUE, + n.knots = 8, cross.validation = TRUE, kappa = 1000, RandDist = "LogN")

Then, using argument init.B we fit the trivariate model with appropriate initial values.

In init.B, the vector of initial values must follow the order: covariates for the recurrent events, terminal event and then biomarker (interactions at the end of each component). The trivariate model due to its complexity is computationally intensive but using the pseudoadaptive quadrature method the time of estimation is reduced (40 minutes on a personal computer with an Intel Core i7 3.40 GHz processor and 8 GB RAM running Windows 7).

The model is defined as follows:

R> modTrivariate <-trivPenal (Surv(time0,time1,new.lesions) --------------- three types of processes: longitudinal outcome, recurrences and terminal event. The treatment effect is not significant neither for the risk of appearance of new lesions nor for the risk of death. As in the bivariate model, the tumor size decreases on average more in the treatment arm C (-0.66, p value < 0.001). We find the baseline WHO performance status 2 having a significant effect on the average population biomarker, the risk of new lesions appearance and death. The coefficients for the association between the biomarker and the recurrences are found not significant and for the link between the biomarker and death only the coefficient related to the random intercept is significantly different from 0 (0.85, p value = 0.004). The variance of the frailty term is found to be not significant.

As for the bivariate model modLongi, goodness-of-fit of the trivariate model can be evaluated using the martingale residuals for the recurrent and terminal events and the residuals related to the biomarker. For this purpose, the following code can be used (Figure 11 The martingale residuals for both recurrences and death processes show slight signs of skewness, their means seem to decrease in time. Using the marginal Cholesky residuals and the conditional residuals we find that the model does not fit the longitudinal data very well.

We compare the trivariate model modTrivariate with the bivariate model modLongi in terms of predictive accuracy for OS using the Brier score. Firstly, we define the prediction time and horizons and calculate the predictions for the bivariate and trivariate models:

R> predtime <-1.0 R> window <-seq(0. the bivariate model has better predictive accuracy than the trivariate model. However, the results are presented for the apparent error calculated using the data from the estimation. For verification of the models for predictions for new patients (not used in the data for the estimation) a cross-validation procedure would be required.

Dynamic predictions

Both, for the bivariate models for longitudinal data and a terminal event and for the trivariate models, the package provides individual predictions of the terminal event. To create patients' profiles, two datasets must be provided, one including the history of the biomarker and the values of covariates at measurement times (dataLongi) and a second one with the information on covariates related to the terminal event (data). In case of the trivariate model, this dataset should also include the history of recurrences.

For the example of the colorectal dataset, we create two profiles of patients that differ from each other by the trajectory of the tumor size. The first patient has a progressive disease with a tumor size increasing in time and the second patient a response to the treatment with a diminishing tumor. Both patients have the same baseline characteristics. In the data creation step, the types of variables must be appropriately defined and coherent with the dataset used for estimation. Firstly, we prepare the data for the bivariate model setting:

R> datapredj_longi <-data.frame(id = 0, year = 0, tumor.size = 0, + treatment = 0, age = 0, who.PS = 0, prev.resection = 0) R> datapredj_longi$treatment <-factor(datapredj_longi$treatment, + levels = 1:2) R> datapredj_longi$age <-factor(datapredj_longi$age, levels = 1:3) R> datapredj_longi$who.PS <-factor(datapredj_longi$who.PS, levels = 1:3) R> datapredj_longi$prev.resection <-factor(datapredj_longi$prev.resection, + levels = 1:2)

For Patient 1 we assume 5 measurements indicating increasing tumor burden:

R> datapredj_longi[1, ] <-c(1, 0, 1.2, 2, 1, 1, 1) R> datapredj_longi[2, ] <-c(1, 0.3, 1.1, 2, 1, 1, 1) R> datapredj_longi[3, ] <-c(1, 0.6, 1.4, 2, 1, 1, 1) R> datapredj_longi[4, ] <-c(1, 0.9, 2.2, 2, 1, 1, 1) R> datapredj_longi[5, ] <-c(1, 1.5, 3.0, 2, 1, 1, 1)
On the contrary, Patient 2 is assumed to have a decreasing size of the tumor:

R> datapredj_longi[6, ] <-c(2, 0, 1.2, 2, 1, 1, 1) R> datapredj_longi[7, ] <-c(2, 0.3, 1.1, 2, 1, 1, 1) R> datapredj_longi[8, ] <-c(2, 0.5, 0.7, 2, 1, 1, 1) R> datapredj_longi[9, ] <-c(2, 0.7, 0.3, 2, 1, 1, 1) R> datapredj_longi[10, ] <-c(2, 0.9, 0.1, 2, 1, 1, 1)
Next, for both patients we prepare the data with information on covariates included in the survival part in the bivariate model:

R> datapredj <-data.frame(id = 0, treatment = 0, age = 0, who.PS = 0, + prev.resection = 0) R> datapredj$treatment <-factor(datapredj$treatment, levels = 1:2) R> datapredj$age <-factor(datapredj$age, levels = 1:3) R> datapredj$who.PS <-factor(datapredj$who.PS, levels = 1:3) R> datapredj$prev.resection <-factor(datapredj$prev.resection, levels = 1:2

) R> datapredj[1, ] <-c(1, 2, 1, 1, 1) R> datapredj[2, ] <-c(2, 2, 1, 1, 1)
We calculate the estimated probabilities of the terminal event given that the patients are alive at time of prediction 1 year and a horizon varying from 0.5 to 2.5 years. We compare the predicted risk of death for the patients by plotting and smoothing the estimations. Additionally, the 95% MC confidence intervals are calculated in order to facilitate the interpretation.

R> pred.joint <-prediction(modLongi, datapredj, datapredj_longi, 1.0, + seq(0.5, 2.5, 0.2), MC.sample = 500) R> plot(pred.joint, conf.bands = TRUE)

The left graph in Figure 8 presents the dynamic predictions for the patients. The patient with a decreasing tumor size (profile 2) has lower probability of death than the patient with the tumor size that increases during the treatment (profile 1). However, considering the MC confidence intervals, this difference is not significant. Thus, in the analyzed example, the biomarker itself does not influence the risk of death significantly. It is of interest if addition of the history of the recurrent event would increase the difference between the profiles. For this purpose, we modify the data datapredj by adding the history of recurrences and calculate the analogous dynamic predictions using the trivariate model modTrivariate. We assume that patient 1 experiences the occurrence of new lesions twice and patient 2 only once.

The right graph in Figure 8 shows the dynamic predictions of death using the trivariate model. As in the case of the bivariate model, patient 1 has an increased probability of death compared to patient 1. Moreover, the difference between the patients is more accentuated than in the bivariate model which is not able to include information on the history of recurrences. By considering the information on the appearance of new lesions, the predicted probabilities are influenced, although the differences between the patients are not significant according to the MC confidence intervals.

Conclusions

Joint models are now a well recognized statistical tool for a complex analysis of correlated data. They provide unbiased estimates compared to the univariate models. The package frailtypack provides several functions for applications of joint models for a longitudinal outcome, recurrent events and a terminal event. Methods of estimation applied in the package produce reliable results proven by simulations for all the implemented models [START_REF] Rondeau | Joint Frailty Models for Recurring Events and Death Using Maximum Penalized Likelihood Estimation: Application on Cancer Events[END_REF], Mazroui et al. 2012[START_REF] Król | 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[END_REF]. Moreover, each function is furnished with several options for better adjustment of a model and accompanied by tools helping with the diagnostic of a model (residuals, individual predictions of random effects) and comparison with other models (LCV, AIC).

The increasing interest in individual predictions of death in the clinical perspective motivated the implementation of model-based dynamic predictions. Function prediction allows the users to calculate the estimated probability of death given the history of a patient. The history, depending on a model, can be a part of the information on observed recurrences, complete information on observed recurrences, complete past measurements of a biomarker or the entire available history of the recurrences and biomarker. Graphical representation of the predictions is a useful tool for comparing probability of death of patients, e.g., with the same characteristics but different past of recurrences and/or biomarker.

Further developments of the frailtypack package will concern the extensions of the existent functions. These extensions will be related to the random effects (distributions, number of random effects for the longitudinal part), association functions for models with a biomarker (provide more forms to choose from), stratification (increase the number of strata) and dynamic predictions for the joint model with two types of recurrent events and a terminal event. Moreover, in order to increase the variety of possible application, it will be of interest to implement left truncation and interval-censoring in the proposed functions.

In the near future, we will develop dynamic predictions of a recurrent event. The conditional probability of a new recurrence occurring in a finite time horizon given the history of the observed events and/or the biomarker and given that an individual is alive at the prediction time, will constitute a useful tool for clinicians, e.g., to explore the chances of developing a new relapse given a patient's characteristics.

Survival functions

Time

Baseline survival function recurrent event terminal event q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 1 2 3 4 -4 -2 0 2 4 Martingale Residuals -Death time q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 1 2 3 4 -4 -2 0 2 4 Martingale Residuals -Recurrences time q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q qq q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -3 -

2 -1 0 1 2 3 -4 -3 -2 -1 0 1 2

Marginal Cholesky residuals

Sample Quantiles q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -10 -5 0 5

  2 and the parameters to estimate are ξ = {r

Figure 1 :

 1 Figure1: Three settings to take into account patient history of recurrent events in prediction. Illustration with two recurrent events.

  R> plot(aggregate(readmission$t.stop, by = list(readmission$id), + FUN = max)[2][, 1], modJoint.gap$martingale.res, ylab = "", + xlab = "time", main = "Rehospitalizations", ylim = c(-4, 4)) R> lines(lowess(aggregate(readmission$t.stop, by = list(readmission$id), + FUN = max)[2][, 1], modJoint.gap$martingale.res, f = 1), lwd = 3, + col = "grey") R> plot(aggregate(readmission$t.stop, by = list(readmission$id), + FUN = max)[2][, 1], modJoint.gap$martingaledeath.res, ylab = "", + xlab = "time", main = "Death", ylim = c(-2, 2)) R> lines(lowess(aggregate(readmission$t.stop, by = list(readmission$id), + FUN = max)[2][, 1], modJoint.gap$martingaledeath.res, f = 1), lwd = 3, + col = "grey")

Figure 3 :

 3 Figure 3: Martingale residuals for rehospitalizations and death against the follow-up time (in days). The grey line corresponds to a smooth curve obtained with lowess.

  R> modJoint.gap.timedep <-frailtyPenal(Surv(time, event) ~cluster(id) + + dukes + charlson + timedep(sex) + chemo + terminal(death), + formula.terminalEvent = ~dukes + charlson + timedep(sex) + chemo, + data = readmission, n.knots = 8, kappa = c(2.11e+08, 9.53e+11), + betaorder = 3, betaknots = 3)

FrailtyFigure 4 :

 4 Figure 4: Individual prediction of the frailties. The size of points corresponds to number of an individual's recurrent events.

Figure 5 :

 5 Figure 5: Predicted probabilities of death for two patients sharing the same history of recurrences and characteristics except of the Dukes' stage. Patient 1 has Dukes' stage A and patient 2, Dukes' stage D. Dashed lines represent the MC confidence intervals.

  spli <-multivPenal(Surv(TIMEGAP, INDICREC) ~cluster(PATIENT) + + v1 + v2 + event2(INDICMETA) + terminal(INDICDEATH), formula.Event2 = + v1 + v2 + v3, formula.terminalEvent = ~v1, data = dataMultiv, + n.knots = c(8, 8, 8), kappa = c(1, 1, 1), initialize = FALSE) Option initialize indicates whether to initialize the parameters (including parameters for the baseline hazard functions) using the estimates of separate models: shared frailty models (for the two types of recurrent events) and a Cox proportional hazard model (for the terminal event). The output of function multivPenal is presented below: R> modMultiv.spli Call: multivPenal(formula = Surv(TIMEGAP, INDICREC) ~cluster(PATIENT) + v1 + v2 + event2(INDICMETA) + terminal(INDICDEATH), formula.Event2 = ṽ1 + v2 + v3, formula.terminalEvent = ~v1, data = dataMultiv, initialize = FALSE, n.knots = c(8, 8, 8), kappa = c(1, 1, 1)) Multivariate joint gaussian frailty model for two survival outcomes and a terminal event using a Penalized Likelihood on the hazard function Recurrences 1: -

  R> plot(modLongi, main = "Hazard function") R> plot(aggregate(colorectalSurv$time1, by = list(colorectalSurv$id), + FUN = max)[2][, 1], modLongi$martingaledeath.res, ylab = "", + xlab = "time", main = "Martingale Residuals -Death", + ylim = c(-4.2, 4.2)) R> lines(lowess(aggregate(colorectalSurv$time1, by = list(colorectalSurv$id), + FUN = max)[2][, 1], modLongi$martingaledeath.res, f = 1), lwd = 3, + col = "grey") R> qqnorm(modLongi$marginal_chol.res, main = "Marginal Cholesky residuals", + xlab = "") R> qqline(modLongi$marginal_chol.res) R> plot(modLongi$pred.y.cond, modLongi$conditional.res, xlab = "Fitted", + ylab = "Conditional residuals", + main = "Conditional Residuals vs. Fitted Values")

Figure 6 :

 6 Figure6: Estimated cross-validated EPOCE computed from the joint models until 2.5 years (left plot) and differences in the estimates with 95% tracking intervals (right plot).

  data.Longi = colorectalLongi, random = c("1", "year"), id = "id", link = "Random-effects", left.censoring = -3.33, recurrentAG = TRUE, n.knots = 7, kappa = c(0.01, 0.7), init.B = c(-0.18, -0.22, -0.24, -0.22, 0.35, -0.23, -0.1, -0.09, -0.12, 0.8, -0.23, 3.02, -0.3, 0.05, -0.02, -0.29, 0.11, 0.74, -0.63), method.GH = "Pseudo-adaptive", n.nodes = 7)Calendar timescaleTrivariate Joint Model for Longitudinal Data, Recurrent Events and a Terminal Event Parameter estimates using a Penalized Likelihood on the hazard functions and assuming left-censored longitudinal outcome Association function: random effects Longitudinal outcome:

Figure 7 :

 7 Figure 7: Error of prediction at t = 1 year and varying window w from 0.1 to 1.5. modLongi -tumor size and death, modTrivariate -tumor size, occurrence of new lesions and death.

Figure 10 :

 10 Figure 10: Log-hazard ratios for sex for modJoint.gap.timedep model.

Figure 11 :

 11 Figure 11: Goodness-of-fit of trivariate model modTrivariate. From the top-left plot: baseline hazard and survival functions, martingale residuals (the grey line corresponds to a smooth curve obtained with lowess), Q-Q normal plot of the marginal Cholesky residuals of the biomarker and the conditional residuals against the fitted values of the biomarker.

  The possible prediction settings including the longitudinal data and considering the whole information available. The top setting corresponds to the model (

	Setting of complete history
	5 longitudinal measures (o) before t	Prediction
		of death
	t	t+w
	Setting of complete history 2 recurrent events (x) and 5	Prediction
	longitudinal measures (o) before t	of death
	t	t+w
	Figure 2:

  New dataset readm_cluster includes clusters with 97 to 107 individuals per group. For definition of a clustered joint model two inner functions are required in frailtyPenal, num.id for the individual level and cluster for the groups:

	R> readmission <-transform(readmission, group = id %% 6 + 1 )
	R> readm_cluster <-subset(readmission,
	+	(t.start == 0 & event == 1) | event == 0)
	R> joi.clus <-frailtyPenal(Surv(t.start, t.stop, event) ~cluster(group) +
	+	num.id(id) + dukes + sex + chemo + terminal(death),
	+	formula.terminalEvent = ~dukes + sex + chemo, data = readm_cluster,
	+
		We consider 6 clusters defined by a new
	variable group:

Acknowledgments

The authors thank the Fédération Francophone de Cancérologie Digestive and Gustave Roussy for sharing the data of the FFCD 2000-05 trial supported by an unrestricted Grant from Sanofi. The authors thank as well the IReSP for financing the project AAP 2012 Plan Cancer

Joint model for longitudinal data and a terminal event

Firstly, we estimate the bivariate joint model for longitudinal data and a terminal event (4). We consider a left-censored biomarker, i.e., the transformed tumor size measurements are not observed below a threshold -3.33 (which corresponds to 'zero' measures in the nontransformed data). The value of the smoothing parameter κ is chosen using a corresponding reduced model, i.e., a Cox model for the terminal event. For a model with a random intercept and slope for the biomarker and the link function being the random effects of the biomarker we use the following form of the longiPenal function using the colorectalLongi dataset and a subset of colorectal containing only the information on the terminal event (colorectalSurv):

R> colorectalSurv <-subset(colorectal, new.lesions == 0) R> modLongi <-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 = 8, kappa = 0.93, method.GH = "Pseudo-adaptive", n.nodes = 7) R> modLongi Call: longiPenal(formula = Surv(time1, state) ~age + treatment + who.PS + prev.resection, formula.LongitudinalData = tumor.size ~year * treatment + age + who.PS, data = colorectalSurv, data.Longi = colorectalLongi, random = c("1", "year"), id = "id", link = "Random-effects", left.censoring = -3.33, n.knots = 8, kappa = 0.93, method.GH = "Pseudo-adaptive", n.nodes = 7) Gaussian quadrature method: Pseudo-adaptive with 7 nodes

In the output of the fitted model, we obtain the estimates for the covariates related to the + data.Longi = colorectalLongi, predtime, window) R> pred_trivariate <-prediction(modTrivariate, data = colorectal, + data.Longi = colorectalLongi, predtime, window)

We prepare data for the Brier score computations: data_surv including data on survival (variables of time and state), survj_biv and survj_tri for predicted probabilities of survival of subjects that are alive at time of prediction using the bivariate and trivariate model, respectively:

R> data_surv <-colorectalSurv[, c(1, 3, 9)] R> predictions_bivariate <-as.matrix(pred_bivariate$pred) R> survj_biv <-+ as.matrix(predictions_bivariate[data_surv$time1 >= predtime, ]) R> survj_biv <-1 -cbind(0, survj_biv) R> predictions_trivariate <-as.matrix(pred_trivariate$pred) R> survj_tri <-+ as.matrix(predictions_trivariate[data_surv$time1 >= predtime, ]) R> survj_tri <-1 -cbind(0, survj_tri)

Prediction error curves can be obtained by calculating the Brier score for several time points. Using the joint model estimated with frailtypack, the Brier score can be computed using the package pec (version 2.4.9, available at https://CRAN.R-project.org/src/contrib/ Archive/pec/) and its modified function pec included in the supplementary material. The following code is used to calculate and plot the Brier score for the bivariate and trivariate models: R> library("pec") R> library("prodlim") R> source("ipcw_2-9_modif.R") R> source("pecMethods_2-9_modif.R") R> BrierScore <-pec(list("Bivariate" = survj_biv, "Trivariate" = survj_tri), + formula = Surv(time1, state) ~1, + data = data_surv[data_surv$time1 >= predtime, ], + cens.model = "marginal", data.cens = data_surv, exact = FALSE, + times = fwindow, ptime = predtime, reference = FALSE) R> plot(fwindow, BrierScore$AppErr$Bivariate[-1], pch = 6, main = "", + ylab = "Prediction error", xlab = "Years", ylim = c(0, 0.3), + xlim = c(0.9, 2.7), col = "black", type = "l", lwd = 2, axes = FALSE, + lty = 1) R> points(fwindow, BrierScore$AppErr$Trivairate[-1], pch = 15, col = "blue", + type = "l", lwd = 2, lty = 2) R> legend(2.1, 0.3, legend = c("modLongi", "modTrivariate"), bty = "n", + cex = 1.1, lty = c(1, 2), lwd = 2, col = c("black", "blue")) R> axis(1, at = seq(0.9, 2.7, by = 0.3)) R> axis(2, at = seq(0, 0.3, by = 0.1)) R> abline(v = 1, lty = 2, lwd = 2, col = "gray35")

Dynamic predictions for survival using modTrivariate

Time t Predicted cumulative probability of event between 1 and time t profile 1 profile 2 Figure 8: Predicted probabilities of death for two patients sharing the same characteristics and different history using information from the biomarker only (on the left) and the biomarker and recurrent events (on the right). The prediction time is set to 1 year and the horizon changed from 0.5 to 2.5 years. Dashed lines represent the MC confidence intervals.

R> datapredj <-data.frame(time0 = 0, time1 = 0, new.lesions = 0, id = 0, + treatment = 0, age = 0, who.PS = 0, prev.resection = 0) R> datapredj$treatment <-factor(datapredj$treatment, levels = 1:2) R> datapredj$age <-factor(datapredj$age, levels = 1:3) R> datapredj$who.PS <-factor(datapredj$who.PS, levels = 1:3) R> datapredj$prev.resection <-factor(datapredj$prev.resection,

Then, we calculate the predictions and plot the results with the 95% MC confidence intervals: R> pred.joint2 <-prediction(modTrivariate, datapredj, datapredj_longi, 1.0, + seq(0.5, 2.5, 0.2), MC.sample = 500) R> plot(pred.joint2, conf.bands = TRUE)

A. Summary of the package frailtypack