Olivier Bouaziz

Assessing model prediction performance for the expected cumulative number of recurrent events

Keywords: Recurrent events, Prediction assessment, Right-censoring, Terminal event, Brier Score

In a recurrent event setting, we introduce a new score designed to evaluate the prediction ability, for a given model, of the expected cumulative number of recurrent events. This score allows to take into account the individual history of a patient through its external covariates and can be seen as an extension of the Brier Score for single time to event data but works for recurrent events with or without a terminal event. Theoretical results are provided that show that under standard assumptions in a recurrent event context, our score can be asymptotically decomposed as the sum of the theoretical mean squared error between the model and the true expected cumulative number of recurrent events and an inseparability term that does not depend on the model. This decomposition is further illustrated on simulations studies. It is also shown that this score should be used in comparison with a reference model, such as a nonparametric estimator that does not include the covariates. Finally, the score is applied for the prediction of hospitalisations on a dataset of patients suffering from atrial fibrillation and a comparison of the prediction performances of different models, such as the Cox model or the Aalen Model, is investigated.

Introduction

Recurrent event data are often encountered in follow-up studies. They can be seen as a generalisation of the standard time to event data, where individuals may experience the same event repeatedly over time. Typical examples may include HIV studies where patients can experience repeated opportunistic infections, remission data from Leukemia patients who can experience multiple relapses, repeated seizures for epileptic patients, or hospitalisation data where the events of interest are the hospitalisations. In those studies, the focus might be on assessing the effect of covariates on the risk of recurrences or on predicting the future recurrences. The first model to deal with recurrent event data was the Andersen-Gill model [START_REF] Andersen | Cox's regression model for counting processes: a large sample study[END_REF] which was further extended by [START_REF] Lin | Semiparametric regression for the mean and rate functions of recurrent events[END_REF] to account for possibly dependent jumps of the recurrent event process. Further models were developed such as in [START_REF] Prentice | On the regression analysis of multivariate failure time data[END_REF], [START_REF] Cook | Marginal analysis of recurrent events and a terminating event[END_REF], [START_REF] Ghosh | Marginal regression models for recurrent and terminal events[END_REF], [START_REF] Ghosh | Semiparametric analysis of recurrent events data in the presence of dependent censoring[END_REF] or [START_REF] Andersen | Modeling marginal features in studies of recurrent events in the presence of a terminal event[END_REF] where the last four papers incorporate the presence of a terminal event in the estimation procedure, or using random effects such as in [START_REF] Hougaard | Analysis of multivariate survival data[END_REF], [START_REF] Liu | Shared frailty models for recurrent events and a terminal event[END_REF], or [START_REF] Rondeau | Joint frailty models for recurring events and death using maximum penalized likelihood estimation: application on cancer events[END_REF]. In particular, in [START_REF] Cook | Marginal analysis of recurrent events and a terminating event[END_REF], [START_REF] Ghosh | Marginal regression models for recurrent and terminal events[END_REF], [START_REF] Ghosh | Semiparametric analysis of recurrent events data in the presence of dependent censoring[END_REF], [START_REF] Andersen | Modeling marginal features in studies of recurrent events in the presence of a terminal event[END_REF] the authors focused on the estimation of the expected cumulative number of recurrent events. This is a marginal quantity that computes the expectation of the number of experienced events of an individual before any time point. This quantity is particularly interesting as it summarises the evolution of the recurrent event process with time. In the presence of a terminal event, it also includes the fact that when the terminal event occurs the patient can no longer experience any further recurrent events.

In some studies the focus is more on the predictiveness ability of a model rather than on the interpretation of the covariates effects. This is the case when clinicians aim at predicting the future repeated events in order to offer the best medical care. Being able to predict the future recurrences of any patients on a short time period also allows to predict the future burden of the disease over the patient's life. Moreover, a predictive model can be an important tool for making medical decisions but also for communicating with the patients about the future course of his/her disease. For instance, in [START_REF] Schroder | Recurrent event survival analysis predicts future risk of hospitalization in patients with paroxysmal and persistent atrial fibrillation[END_REF] the authors studied patients with atrial fibrillation, a well known cardiac disease, in an attempt to predict the future hospitalisations of patients due to their disease. Since the patients suffering from this disease are usually old (the median age in the study was 63 years) and since atrial fibrillation can be a severe disease in some cases, those patients were also at risk of death. Several covariates were collected and a prediction of the expected cumulative number of recurrent events over time was performed using a Cox model with dependence on prior counts.

While such models are certainly of interest for clinicians, it is important to propose relevant diagnosis tools that can evaluate the prediction performance of the proposed model. There already exists several indicators for prediction performances in the standard context of time to event data with only one event per individual. The Brier score was developed in [START_REF] Graf | Assessment and comparison of prognostic classification schemes for survival data[END_REF] and in [START_REF] Gerds | Consistent estimation of the expected brier score in general survival models with right-censored event times[END_REF], which basically is a score for computing the mean squared error of the time to event in the presence of censoring. This score was further developed to deal with random effect models in [START_REF] Van Oirbeek | Exploring the clustering effect of the frailty survival model by means of the brier score[END_REF], or to evaluate the performance of dynamic prediction models in [START_REF] Schoop | Measures of prediction error for survival data with longitudinal covariates[END_REF] where the information available from a longitudinal covariate is updated at each time point. Note also that other types of predictive accuracy measures exist, called discrimination measures, such as the C-index (see [START_REF] Harrell | Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors[END_REF], [START_REF] Gerds | Estimating a time-dependent concordance index for survival prediction models with covariate dependent censoring[END_REF]) or the time dependent ROC curve and area under the curve (see for instance [START_REF] Heagerty | Survival model predictive accuracy and roc curves[END_REF]).

In this paper, the aim is to derive a predictive accuracy measure for recurrent events models where the focus is on predictiveness rather than discrimination. The quantity of interest is solely the expected cumulative mean number of recurrent events. Since no mean squared error measure, such as the Brier score, exists in the context of recurrent events, the goal of this work is to fill in this gap by deriving a new score of this type for recurrent events, which also accommodates for the presence of a terminal event. In this work, we show that this score reduces to the Brier score when only one event per individuals can occur and hence can be seen as a direct generalisation of the standard Brier score. Also, since our prediction criterion focuses on the marginal quantity of the expected cumulative number of recurrent events, it provides a summary score that takes into account the prediction of all recurrent events. In the context of a terminal event, it also incorporates the quality of prediction of the terminal event.

In Section 2.1, we introduce the general prediction criterion for recurrent events, denoted MSE. In Sections 2.2 and 2.3, we separate the modelling assumptions in cases where no terminal events are observed and in the presence of a terminal event. In Section 2.4, we present some existing estimators for the expected cumulative number of recurrent events. In Section 3, we derive the main theoretical results of this paper. We first introduce a theoretical criterion and show that it can be decomposed into an inseparability term and an imprecision term, similarly to the results in [START_REF] Gerds | Consistent estimation of the expected brier score in general survival models with right-censored event times[END_REF]. The former does not depend on the model and cannot be removed while the latter is exactly the mean squared error between the recurrent event process and the prediction model of the expected cumulative mean number of recurrent events. We then show that our prediction criterion asymptotically converges towards the theoretical criterion. In Section 4, we demonstrate that when individuals can only experience one event, our prediction criterion is equivalent to the standard Brier score. In Section 5.1, a simulation study is conducted. First, the decomposition between inseparability and imprecision terms is illustrated. As the inseparability is, by far, the dominant term, we then recommend to consider as a prediction score, the difference of MSEs between the considered model and a reference model. Second, we illustrate how this prediction score can be used in order to compare prediction models. In Section 6, the atrial fibrillation dataset is studied. We show that the model with dependence on prior counts (which is a multi-state model), stratified with respect to atrial fibrillation type, provides the best prediction performance among all other models considered.

2 Prediction criterion for the expected cumulative number of recurrent events

The prediction criterion in the general framework

In this section we present a prediction criterion for two different recurrent event settings under right-censoring. A scenario with right-censoring only (Section 2.2) and a scenario with the inclusion of a terminal event (Section 2.3) are investigated. In each case, a counting process of interest N * (t) is defined which counts the number of recurrent events that have occurred before time t. We assume that a multivariate external time dependent covariate vector X(t) (see [START_REF] Kalbfleisch | The statistical analysis of failure time data[END_REF] for the definition of external covariates) is observed and we define M a class of bounded functions depending on t and X(t). For each t ≥ 0, we define X t , the support of the process {X(u) : 0 ≤ u} at the time point u = t. In both scenarios, censoring can occur such that the observed recurrent event process is N (t) := N * (t ∧ C), where C is a censoring variable and a ∧ b represents the minimum between a and b. We also define τ the endpoint of the study. On the basis of i.i.d. replications (N i (t), X i (t) : 0 ≤ t ≤ τ ), we will assume that an estimator µ ∈ M of the expected cumulative number of recurrent events µ * (t | X(t))

:= E[N * (t) | X(u) : 0 ≤ u ≤ t]
is available. We will say that this estimator is consistent if there exists µ ∈ M such that for all t ≤ τ ,

sup x∈Xt | µ(t | x) -µ(t | x)| → 0, in probability as n → ∞.
The main goal of this paper is to develop a new mean squared error criterion designed to evaluate the performance of this estimator in different scenarios. In Sections 2.2 and 2.3 we propose to evaluate the prediction ability of a given estimator µ ∈ M, through the following criterion:

MSE(t, µ) = 1 n n i=1 t 0 dN i (u) 1 -Ĝ(u-| X i (u)) -µ(t | X i (t)) 2 , ( 1
)
where Ĝ is an estimator of G, the conditional cumulative distribution function of the censoring variable C given X(•). The notation Ĝ(u-| X i (u)) indicates the left limit of the function Ĝ at u. We will assume uniform consistency of this censoring estimator in the following way.

Assumption 1 Let G be a model for the conditional censoring distribution. We say that Ĝ is a uniformly consistent estimator for G ∈ G if for all t ≤ τ ,

sup x∈Xt | Ĝ(t | x) -G(t | x)| → 0, in probability as n → ∞.
Presentations of different estimators for G, depending on the considered scenario, are discussed in Sections 2.2 and 2.3. We now introduce a theoretical criterion that would be available if the censoring distribution was known. For some function µ ∈ M, let:

MSE(t, µ) = E t 0 dN (u) 1 -G(u-| X(u)) -µ(t | X(t)) 2 . (2)
The crucial idea behind our criterion (1) comes from the fact that

E t 0 dN (u) (1 -G(u-| X(u))) = E[µ * (t | X(t))], (3)
a relationship that is proved for each scenario in Sections 2.2 and 2.3. In Section 3, we provide theoretical results for both scenarios that justify the appropriateness of the proposed criterion.

In particular, Proposition 1 of Section 3 shows that the theoretical criterion can be decomposed in the following way:

MSE(t, µ) = E µ * (t | X(t)) -µ(t | X(t)) 2 + A(t),
with A(t) not depending on µ. The first term is an imprecision term and the second term is an inseparability (or residual) term that does not depend on the chosen model. It should be noted that this kind of result is similar to the imprecision/inseparability decomposition of the Brier score (see [START_REF] Gerds | Consistent estimation of the expected brier score in general survival models with right-censored event times[END_REF]). On the other hand, Proposition 2 of Section 3 states that if µ is a consistent estimator for some µ ∈ M then as n tends to infinity, our empirical criterion ( 1) is asymptotically equivalent to the theoretical criterion (2) evaluated at µ. Then, this latter criterion will be minimal if µ = µ * from Proposition 1. In case µ is not equal to µ * , a bias will occur which will be equal to the squared expectation between µ and µ * . The precise assumptions under which Proposition 1 and Proposition 2 hold are specified in each situation.

In what follows, the two scenarios, with right censoring only and with this inclusion of a terminal event, are presented. It should be noted that the estimated quantity µ * , is based on the rate function, as defined in [START_REF] Lin | Semiparametric regression for the mean and rate functions of recurrent events[END_REF] or [START_REF] Cook | The statistical analysis of recurrent events[END_REF]. In this modelling approach, the conditioning is performed on the covariates but, contrary to the intensity function, we do not condition on the entire history of the process. As a result, our criterion (1) works for general recurrent event processes and, in particular, they are not restricted to Poisson processes.

Situations with right-censoring and no terminal event

As before, we introduce the counting process of interest N * (t) which counts the number of recurrent events that have occurred before time t. Let X(t) be some multivariate external covariate (see [START_REF] Kalbfleisch | The statistical analysis of failure time data[END_REF]) that is allowed to depend on t. For each t ≥ 0, we define X t , the support of the process {X(u) : 0 ≤ u} at the time point u = t. The regression modelling approach considered here is defined by:

E[dN * (t) | X(t)] = λ * (t | X(t))dt, (4)
where λ * (t | X(t)) is the true rate function. Note that this definition is very general and does not make any assumption on the recurrent event process such as independent increments or a Poisson assumption. See [START_REF] Lin | Semiparametric regression for the mean and rate functions of recurrent events[END_REF] for more details on this model. Let t 0 λ * (u | X(u))du be the true cumulative rate function. In the absence of a terminal event, this cumulative function has a direct interpretation as the expected cumulative number of recurrent events given the covariate process:

t 0 λ * (u | X(u))du = µ * (t | X(t))
where as defined at the beginning of Section 2, µ

* (t | X(t)) = E[N * (t) | X(u) : 0 ≤ u ≤ t].
In the presence of censoring, a variable C is observed such that the observed recurrent event process is now N (t) = N * (t ∧ C). We assume independent censoring (see [START_REF] Lin | Semiparametric regression for the mean and rate functions of recurrent events[END_REF]) in the following way:

E[dN * (t) | X(t)] = E[dN * (t) | C, X(t)].
This assumption implies that C does not convey any additional information on the probability of a jump of the recurrent event process. Under this assumption, we have

E[dN (t) | I(C ≥ t), X(t)] = I(C ≥ t)λ * (t | X(t))dt, (5)
where I(•) denotes the indicator function. This last equation justifies the use of our criterion [START_REF] Andersen | Modeling marginal features in studies of recurrent events in the presence of a terminal event[END_REF] since

E[dN (t) | X(t)] = (1-G(t-| X(t)))λ * (t | X(t)
)dt and therefore Equation (3) holds. Next, we assume Assumption 1 and we make the following additional assumption.

1. ∀t ∈ [0, τ ], P[C ≥ t | X(t)] ≥ c almost surely, 2. N (τ )
is almost surely bounded by a constant.

Note that condition 1. was also assumed in [START_REF] Gerds | Consistent estimation of the expected brier score in general survival models with right-censored event times[END_REF]. It is stronger than simply assuming P[C ≥ τ ] ≥ c. Condition 2. is standard for recurrent event data, see for instance [START_REF] Lin | Semiparametric regression for the mean and rate functions of recurrent events[END_REF]. Finally, note that through Equation ( 5), conditions 1. and 2. imply that

E[µ * (τ | X(τ ))] < ∞.
On the basis of i.i.d. replications (N i (t), X i (t) : 0 ≤ t ≤ τ ), let µ ∈ M be an estimator of µ * where M is a class of models that are assumed to be bounded. We propose to evaluate the prediction ability of this estimator through criterion MSE(t, µ) defined in Equation [START_REF] Andersen | Modeling marginal features in studies of recurrent events in the presence of a terminal event[END_REF]. This criterion involves an estimator of G, the conditional cumulative distribution function of the censoring variable. If C and X(•) are independent, one can estimate G using the empirical cumulative distribution function of the censored variables since all these variables are observed. If C depends on X the conditional distribution of C must also be modelled. This can be done using kernel based estimators such as the Nadaraya-Watson estimator for the binary response variable I(C ≤ t) or extensions of this model. For instance, in [START_REF] Hall | Methods for estimating a conditional distribution function[END_REF], a local logistic method and an adjusted Nadaraya-Watson estimator are proposed. If the dimension of the covariates that are assumed to depend on the censoring distribution is too large, then a dimension reduction technique can first be employed, for example through a Single-Index-Modelling approach (see [START_REF] Delecroix | Efficient estimation in conditional single-index regression[END_REF]).

Theoretical results on the validity of this criterion are derived in Section 3.

Situations with right-censoring and a terminal event

In this section, we introduce a terminal event T * such that the counting process of interest N * verifies (see e.g. [START_REF] Cook | The statistical analysis of recurrent events[END_REF]):

E[dN * (t) | I(T * ≥ t), X(t)] = I(T * ≥ t)λ * (t | X(t))dt, (6)
where λ * (t | X(t)) is the true rate function. The difference in Model [START_REF] Cook | The statistical analysis of recurrent events[END_REF] with respect to Model (4) is that the recurrent event process N * is stopped by the terminal event. This often occurs in real-data analysis where death is typically the terminal event. Under this model, we observe that

µ * (t | X(t)) = E[N * (t) | X(u) : 0 ≤ u ≤ t] = t 0 S(u | X(u))λ * (u | X(u))du,
where S(t | X(t)) := P[T * ≥ t | X(t)] is the conditional survival function of the terminal event. This implies that in order to define an estimator of µ

* (t | X(t)) = E[N * (t) | X(u) : 0 ≤ u ≤ t]
one usually needs to also model the hazard rate for the terminal event and to derive an estimator of the conditional survival function. As a result our prediction criterion will both take into account the predictive performance of the survival function and of the rate function of N * since, if one of those two estimators behaves poorly, the resulting estimator for µ * is likely to perform badly as well.

The independent censoring assumption is now expressed in the following way:

E[dN * (t) | T * , X(t)] = E[dN * (t) | T * , C, X(t)].
We denote T = T * ∧ C the minimum between terminal event and censoring, Y (t) = I(T ≥ t) the observed at-risk process and N (t) = N * (T ∧ t) the observed counting process. Under the independent censoring assumption, it can be shown that

E[dN (t) | Y (t), X(t)] = Y (t)λ * (t | X(t))dt. (7)
We assume Assumption 1 and we make the following additional assumption.

Assumption 3 We assume that there exists a constant τ > 0 and a constant c > 0 such that

1. ∀t ∈ [0, τ ], P[T ≥ t | X(t)] ≥ c almost surely, 2. N (τ )
is almost surely bounded by a constant.

We also assume that T * is independent of C conditionally on X(•).

Those conditions are standard in the context of regression for recurrent events with a terminal event, see [START_REF] Ghosh | Marginal regression models for recurrent and terminal events[END_REF] for example. Using Equality [START_REF] Cook | Marginal analysis of recurrent events and a terminating event[END_REF] one can easily observe that

E[dN (t) | X(t)] = S(t | X(t))(1 -G(t-| X(t))λ * (t | X(t)
)dt under the independent censoring hypothesis in Assumption 3. We then directly see that Equation (3) holds.

On the basis of i.i.d. replications (N i (t), X i (t) : 0 ≤ t ≤ τ ), let µ ∈ M be an estimator of µ * where M is a class of models that are assumed to be bounded. We propose to evaluate the prediction ability of this estimator through criterion MSE(t, µ) defined in Equation ( 1). This criterion involves an estimator of G, the conditional cumulative distribution function of the censoring variable. If C and X(•) are independent, one can estimate G using the Kaplan-Meier estimator by considering C to be the variable of interest that is incompletely observed due to the terminal event T * . If C depends on X the conditional distribution of C must also be modelled. Several possible models are presented in [START_REF] Gerds | Consistent estimation of the expected brier score in general survival models with right-censored event times[END_REF] such as the Cox model, the Aalen additive model, or the kernel type model of [START_REF] Dabrowska | Uniform consistency of the kernel conditional kaplan-meier estimate[END_REF]. Alternatively, a single-index approach for right-censored data, as in [START_REF] Bouaziz | Conditional density estimation in a censored singleindex regression model[END_REF], or the random survival forest method developed in [START_REF] Ishwaran | Random survival forests[END_REF] can be used.

Theoretical results on the validity of this criterion are derived in Section 3.

Examples of estimators for the expected cumulative number of recurrent events

In this section we present some estimators for the expected cumulative mean number. In order to achieve this goal, one possibility is to first model the rate function λ * and then to use a plug-in estimator to derive the final estimator of µ * . For example, the Cox (see [START_REF] Cox | Regression models and life-tables[END_REF]), Aalen (see [START_REF] Scheike | The additive nonparametric and semiparametric aalen model as the rate function for a counting process[END_REF]) or Accelerated Failure Time (see [START_REF] Lin | Accelerated failure time models for counting processes[END_REF]) models can be used for the estimation of λ * . We note λ and Λ such estimators of λ and Λ. In the context of right-censored data and no terminal event, the estimator of µ * can then be expressed in the following way (see [START_REF] Cook | The statistical analysis of recurrent events[END_REF]):

µ(t | X(t)) = t 0 d Λ(u | X(u)). (8)
In the presence of a terminal event, a common approach is to first model the hazard rate of the terminal event (using again a Cox model for instance) and to derive an estimator of the survival function

S(t | X(t)) = exp(- t 0 λ T * (u | X(u))du)
where λ T * is the estimator of the hazard rate of the terminal event. Then, the final estimator of µ * is (see [START_REF] Cook | The statistical analysis of recurrent events[END_REF] and [START_REF] Andersen | Modeling marginal features in studies of recurrent events in the presence of a terminal event[END_REF]):

µ(t | X(t)) = t 0 S(u | X(u))d Λ(u | X(u)). (9)
Alternative approaches that directly model µ * also exist. In [START_REF] Ghosh | Marginal regression models for recurrent and terminal events[END_REF] the authors consider the following Cox type model: µ * (t | X(t)) = µ 0 (t) exp(X(t) β). A more general approach consists in using a Single-Index-Model for estimating µ * : in [START_REF] Bouaziz | Semiparametric inference for the recurrent events process by means of a single-index model[END_REF] the authors assume the existence of a nonparametric function g and a parameter β such that µ * (t | X(t)) = g(X(t) β). Those two approaches provide a direct estimator of the quantity of interest µ * . However, in the presence of a terminal event, it is no longer possible to disentangle the effects on the recurrent event process or the terminal event. As a result, the regression parameters should be interpreted with caution. See also [START_REF] Cook | The statistical analysis of recurrent events[END_REF] for a discussion about this issue.

Theoretical results

In this section we provide theoretical results on the proposed criterion (1) in the two different contexts that have been previously considered: the scenario with right-censoring only and the scenario with right-censoring and a terminal event. Two results are obtained. The first one is concerned with the theoretical criterion [START_REF] Andersen | Statistical models based on counting processes[END_REF]. It shows that this criterion applied to a function µ ∈ M reduces to the mean squared error between µ and µ * and a term that does not depend on µ. The second result shows the asymptotic consistency between MSE(t, μ) and MSE(t, µ) when the estimator µ is an asymptotically consistent estimator of µ.

Proposition 1 For the scenarios with right-censoring only (Section 2.2) and with right-censoring and terminal event (Section 2.3) we respectively assume Assumption 2, Assumption 3. In both scenarios we also assume independent censoring. We then have for µ ∈ M,

MSE(t, µ) = E µ * (t | X(t)) -µ(t | X(t)) 2 + A(t), (10)
where A(t) ≥ 0 for all t ≥ 0 and A(t) does not depend on µ.

Proposition 2 For the scenarios with right-censoring only (Section 2.2) and with right-censoring and terminal event (Section 2.3) we respectively assume Assumption 2 or Assumption 3. In both scenarios we also assume Assumption 1 and independent censoring. Then, if the estimator µ ∈ M is consistent for µ ∈ M, we have

sup t≤τ MSE(t, μ) -MSE(t, µ) → 0, in probability as n → ∞.
The proofs are provided in the Appendix, in Sections 8.1 and 8.2, with an explicit expression of A(t) in both scenarios.

Link with the Brier score

The Brier score (see [START_REF] Gerds | Consistent estimation of the expected brier score in general survival models with right-censored event times[END_REF]) is a popular criterion to evaluate the prediction performance of a regression model for the conditional survival function in the context of right-censoring when only a single event can be observed per individual. We show in this section that if we use our criterion when individuals can only experience one event at most, then our theoretical criterion denoted MSE and defined as follows

MSE (t, π) = E 1 - t 0 dN (u) 1 -G(u-| X(u)) -π(t | X(t)) 2 , (11)
where π = 1 -µ, reduces to the theoretical Brier score up to a term that does not depend on the model M. Note that when individuals can only experience one event, the recurrent event process reduces to

N * (t) = I(T * ≤ t) and µ * (t | X) = E[N * (t) | X]
is the conditional cumulative distribution function of T * . Since the Brier score has been designed for the prediction of the conditional survival function, we have simply rewritten our criterion in Equation ( 2) such that π represents the model for the conditional survival function.

We first recall that the theoretical Brier score is defined as (see Equation (1) from [START_REF] Gerds | Consistent estimation of the expected brier score in general survival models with right-censored event times[END_REF]):

MSE Brier (t, π) = E [(I(T * > s) -π(t | X)) 2 ],
where the expectation is taken with respect to the joint distribution of T * and X. For simplicity, the covariate X is not time dependent in the formula, as presented in the paper of [START_REF] Gerds | Consistent estimation of the expected brier score in general survival models with right-censored event times[END_REF], but the results presented in this section are still valid for time dependent covariates. In their work, the authors show similar results as Propositions 1 and 2 of the present paper when the aim is to provide a prediction of the survival function

S(t | X) = P[T * > t | X]
. Note that we have suppressed the dependency with respect to S in the definition of the Brier score to stay consistent with the notations used throughout this paper. Also, in the definition of the Brier score, π plays the role of 1 -µ in the present paper, that is, it is the limiting function of a proposed conditional survival estimator Ŝ(t | X). We have the following result.

Proposition 3

We assume that only one event per individual can be experienced, that is, N (t) = I(T ≤ t, ∆ = 1), with T = T * ∧ C is the observed time, T * is the true event time, C is the censoring variable and ∆ = I(T * ≤ C) is the censoring indicator. Then, under independent censoring, we have:

MSE (t, π) = MSE Brier (t, π) + B(t),
where B(t) ≥ 0 for all t ≥ 0 and B(t) does not depend on π.

The proof is provided in the Appendix, in Section 8.3, with an explicit expression of B(t). Since B(t) does not depend on the model π, those two criterions are completely equivalent. In particular, if one considers the difference of a regression model to a reference model, such as a model that does not include covariates, then the two criterions will provide exactly the same values since the B(t) term will cancel out in the difference. As we will see in the next section, comparing a model to a reference is typically what we recommend in practice. Since recurrent events are a generalisation of the single-event per individual situation, our criterion can be seen as an extension of the Brier score for recurrent events.

Simulations

A scenario with right-censoring and no terminal events

For i = 1, . . . , n, we first simulate a two-dimensional covariate vector X i = (X i,1 , X i,2 ) with X i,1 a Bernoulli variable with parameter 0.5 and X i,2 a Gaussian variable with expectation 2 and standard deviation 0.5. Conditional on X i , the recurrent events are generated from a non-homogeneous Poisson process with rate λ * (• | X i ) that follows a Cox model with Weibull baseline and a two-dimension time independent covariate. More specifically,

λ * (t | X i ) = λ 0 (t) exp(θ 0 X i ), λ 0 (t) = α β t β α-1
, with α = 2 the shape parameter, β = 0.39 the scale parameter and θ 0 = (log(2), log(0.5)) .

Recurrent events generated under this simulation setting follow Equation ( 4) with the true expected number of recurrent events equal to

µ * (t | X i ) = t 0 λ * (t | X i )du = t β α exp(θ 0 X i ). ( 12
)
We further simulate a censoring variable C i that follows a uniform distribution on [0, 3]. Using those parameters, we observe 0 or 1 recurrent event for 30% of the individuals, less or equal than 5 events for 54% of the individuals, and less or equal than 12 events for 77% of the individuals.

On average, we observe approximately 8 recurrent events per individual.

Based on a single simulated sample, we first illustrate Propositions 1 and 2 when the class of models M that contains µ assumes no effect of the covariates on the occurrence of the recurrent events. For this purpose, we independently simulate a training and a test samples. The training sample is used to compute μtrain based on Equation [START_REF] Cox | Regression models and life-tables[END_REF] where Λ does not depend on X and is simply the Nelson-Aalen estimator:

Λ(t) = 1 n n i=1 t 0 dN i (u) 1 -Ĝ(u-) ,
with Ĝ the Kaplan-Meier estimator of C. We then compute MSE(t, µ train ) from Equation (1) based on the test sample of size n test , that is, the computation is performed on a sample (N 1 (•), X 1 ), . . . , (N ntest (•), X ntest ) independent of the training sample. This quantity should provide an accurate estimation of MSE(t, µ) from Proposition 2. We then compute the imprecision

term E µ * (t | X) -µ(t | X) 2
in Equation ( 10) using the true value of µ * (t | X(t)), replacing µ by μtrain and replacing the expectation by its empirical sum. In other words, we compute

1 n test ntest i=1 µ * (t | X i ) -μtrain (t | X i ) 2 ,
which should give a very accurate estimation of the imprecision term. The A(t) term is exactly computed based on its explicit expression (see Section 8.4 in the Appendix). The decomposition of the MSE between its imprecision and inseparability terms is displayed in Figure 1 using n train = 200 and n test = 1, 000. In Equation ( 1) the Kaplan-Meier estimator of G was computed from the combination of the training and test samples. The solid line represents the estimated MSE while the dotted and dashed lines represent the inseparability and imprecision terms, respectively. The inseparability term is seen to be very close to the MSE. In contrast, the imprecision term, which clearly is not null here since the estimated model uses no covariates, is relatively small as compared to the other two terms. This plot suggests that it might be difficult to compare different models as the inseparability term is dominant in the decomposition of the MSE, which implies that two MSEs computed from two different models will tend to look very similar (for instance, for t = 2.5, the inseparability term represents approximately 84% of the value of the MSE). As a result, we advocate the use of a reference or null model and to compute the score of a given model as the difference between the MSE of the reference and the MSE of this model. Therefore, this score will represent the prediction gain of the model as compared to the null model. A typical choice of the null model is the one that uses no covariates. The score is defined as

Score(t, µ, µ 0 ) = MSE(t, µ 0 ) -MSE(t, µ), ( 13
)
where µ is the evaluated model and µ 0 is the reference model. Those models will usually be implemented based on a training sample. The idea behind this score is that the inseparability term will cancel out in the difference, and the score is therefore equal to the difference between the imprecision terms of the two models. An illustration of this score is presented in Figure 2. Using the same simulation setting as before, we compare the performance of four different models based on the Cox and Aalen models, implemented using either only the first covariate X i,1 or the two covariates X i,1 and X i,2 . Figure 2 displays the prediction scores for 100 training samples of size 50 and a unique test sample of size 1, 000. The reference model is the one that uses no covariates and is estimated from the Nelson-Aalen estimator. Roughly, we see that all models have a better prediction performance than the Nelson-Aalen estimator as time increases especially from time equal to 1.5 and time equal to 2, for the models with one covariate and the models with two covariates, respectively. The models with two covariates clearly outperform the ones with one covariate with a slightly better performance for the Cox model as compared to the Aalen model. This is further illustrated in Table 1 where we compare the mean score of those four different models based on 500 training samples of size 20 and 50 and one single test sample of size 1, 000. We clearly observe that the correctly specified Cox model with two covariates outperforms all other models on average, for all time points and sample sizes. However, it tends to have a slightly bigger standard deviation, especially for n train = 20 and t = 2 or t = 2.9.

n train = 20 n train = 50 t = 1 t = 2 t = 2.9 t = 1 t = 2 t =

A scenario with right-censoring and a terminal event

We now consider a simulation scenario which also includes a terminal event. The recurrent event process and its covariates are simulated in the same manner as in the previous section, with the same parameter values. The censoring variable is simulated following a uniform variable on [0, 8]. The terminal event is simulated according to a Cox model with baseline following a Weibull distribution with shape parameter equal to 5 and scale parameter equal to 1.8. This Cox model also includes the same two covariates as for the recurrent event process with the same effects on the hazard function (i.e. the effects are equal to log(2), log(0.5) for the Bernoulli and Gaussian covariates respectively). This setting leads to 8.5 events per individual on average, with 26%, 50%, 77% of individuals that experience less than or equal to 3, 7 and 12 events, respectively. On average 28% of individuals are censored.

We estimate the expected cumulative number of recurrent events based on Equation ( 9) where in the formula, the Breslow estimator is used to estimate the conditional survival function of the terminal event, if the estimation model for the terminal event includes covariates. In other words:

Aalen
S(t | X) = exp - t 0 exp(X βT * )d Λ T * 0 ,
with βT * is the estimated regression parameter from the Cox model for the terminal event and ΛT * 0 its corresponding baseline estimator known as the Breslow estimator. If the terminal event model does not contain any covariates, then the Kaplan-Meier estimator is used instead. As previously, we use the score defined in Equation [START_REF] Ghosh | Semiparametric analysis of recurrent events data in the presence of dependent censoring[END_REF] to evaluate the quality of prediction of a model where the reference model µ 0 is defined as

µ 0 (t) = t 0 Ŝ(u)d Λ(u), (14)
with Ŝ the Kaplan-Meier estimator of the terminal event and Λ the Nelson-Aalen estimator of the recurrent event process. The same score is also used for the prediction of the survival function with the Kaplan-Meier estimator as the reference model. We consider four different regression models: a correctly specified model that includes the two covariates for both the recurrent event process and the terminal event in two Cox models, a model where the Gaussian covariate is missing for the Cox model of the terminal event (but the Cox model of the recurrent event is correctly specified) and a model where the Gaussian covariate is missing for both Cox models.

In Figure 3, we simulated 100 training samples each of size 800 and we evaluated the prediction score on a unique test sample of size 1, 000. In the bottom panel, we see that including the two covariates in the survival model increases the prediction performance as compared to the model with only one covariate. Also, for both models, the gain in terms of prediction is more important for small time points and is reduced after time 2 approximately, as compared to the Kaplan-Meier estimator. This loss in terms of prediction efficiency of the survival function for large time points impacts the prediction of the expected cumulative number of recurrent events. In the top panel, we see that adding the correct covariates in the Cox models of the survival function and of the recurrent event models increases the prediction performances. After time 2, the gain in terms of the performance prediction of the expected cumulative number of recurrent events slightly decreases due to the loss of efficiency in the prediction of the survival function. Table 2 provides the mean score of those three different models based on 500 training samples of size 100, 200, 400 and 800 and one single test sample of size 1, 000. We see the same trend as in Figure 3 for all sample sizes. Clearly, increasing the sample size does not provide much gain in terms of average especially for small time points but it does reduce the variability of the predictors. In Figure 4, we also looked at models where the terminal event is treated as a censored variable. One model does not include any covariate and is estimated using the Nelson-Aalen estimator. The other includes the two correct covariates and is based on the Cox model. The reference model is the same as previously (see Equation ( 14)), that is with no covariates but with the terminal event properly taken into account. We clearly see that both models provide very poor predictions. This is due to the fact that when the terminal event is treated as a censoring variable, the corresponding Nelson-Aalen estimator predicts too many recurrences over time. 2: Means and standard deviations (in bracket) over 500 simulations for the prediction score of the expected number of recurrent events in the presence of a terminal event. The results are presented for the same three models as in Figure 3. The reference model uses no covariates and was estimated from the non-parametric estimator in Equation [START_REF] Ghosh | Marginal regression models for recurrent and terminal events[END_REF]. Large values indicate better predictive performances.

Real data analysis: the Atrial Fibrillation dataset

In this section, we analyse a dataset on patients with atrial fibrillation (AF). The aim is to compare different regression models for the prediction of the expected cumulative number of atrial fibrillation hospitalisations, using the prediction score developed in this work. Patients were enrolled from January 1st 2008 to December 1st 2012 in the "Atrial Fibrillation Sur- ignored and treated as a censoring variable. The data were generated as in Figure 3. The expected cumulative number recurrent of events was predicted using the Cox model with one or two covariates. The reference model uses no covariates and was estimated from the Nelson-Aalen estimator from the non-parametric estimator (see Equation ( 14)).

vey-Copenhagen (ATLAS-CPH)" from both the in-and outpatient clinics at the Department of Cardiology at University Hospital Copenhagen, Hvidovre, Denmark. All patients were previously diagnosed with AF and were categorised at baseline, into either suffering from paroxysmal atrial fibrillation (PAF) or persistent atrial fibrillation (PeAF). PAF was defined as at least one recorded AF episode with spontaneous conversion to sinus rhythm, no valvular AF, and excluding other temporal forms of AF. PeAF was defined as at least one recorded episode of AF lasting more than 7 days, or where either medical or electrical cardioversion was needed to restore sinus rhythm (in accordance with the Danish Cardiology Society AF guidelines at this time). Other inclusion criteria were age > 18 years, recent (< 1 month) AF documented via either standard 12-lead electrocardiogram (ECG) or home monitoring and ability to give oral and written consent. In total, 174 patients were enrolled with 50 PAF patients and 124 PeAF patients. Time is measured in days, with a mean follow-up duration of 1 279 days. In terms of observed events, the patients experienced a total of 325 AF hospitalisations, with 305 AF hospitalisations in the PeAF group and 20 in the PAF group. A terminal event was defined as either progression to permanent AF or as the occurence of death. In the dataset, 45 patients experienced a terminal event and the remaining 129 patients were censored. Finally, in top of the AF type, the dataset also includes 11 additional variables: gender, age, alcohol consumption (with two levels 0 -5 and > 5), tobacco consumption (with three levels "never smoked", "ex-smoker", "current smoker"), presence of hypertension, heart failure, valvular heart disease, ischemic heart disease, diabetes, COPD, antiarrythmic medication. The data are presented in great details in [START_REF] Schroder | Recurrent event survival analysis predicts future risk of hospitalization in patients with paroxysmal and persistent atrial fibrillation[END_REF]. Note also that the data are fully available from the Plos One website.

In [START_REF] Schroder | Recurrent event survival analysis predicts future risk of hospitalization in patients with paroxysmal and persistent atrial fibrillation[END_REF], the authors analysed the data using a multi-state approach with four possible states: no experience of recurrent events yet, 1 recurrent event, 2 or more recurrent events and the absorbing state for the terminal event. The transition intensities were assumed to be proportional with each other using a Cox model, where the number of previous recurrent events was included in the model. Those types of multi-state models with terminal event are described for instance in [START_REF] Cook | The statistical analysis of recurrent events[END_REF] (see Section 6.6.4 of their book). Those analyses showed a high significant effect of the AF type, the number of previous recurrent events (p-values < 10 -4 ) and of age (p-value= 0.0253) for the risk of future AF hospitalisations. The effect of diabetes had a p-value equal to 0.0955. All other variables were assessed as non significants (p-values> 0.2). A Cox model was also implemented for the terminal event using the multi-state approach (that is including the effect of previous AF hospitalisations through a proportional effect) with all variables. Only the age variable was significant (in the multivariate Cox model, the hazard rate was equal to 1.05 and the p-value was equal to 0.0016). In this previous work, the authors then decided to only include the covariates AF type and age, with a proportional effect of the number of previous AF hospitalisations for the modelisation of the recurrent event process. For the terminal event model, they only included the age variable. Based on those models, it is then possible to produce predictions for the expected cumulative number of future AF hospitalisations, for a given patient based on his/her characteristics. Using the prediction score developed in this paper, we will compare the performance of the model used in [START_REF] Schroder | Recurrent event survival analysis predicts future risk of hospitalization in patients with paroxysmal and persistent atrial fibrillation[END_REF] with several other possible models. Since diabetes is a known risk factor for AF, we will also consider models with this variable, along with AF type and age. As in the simulation section, the prediction score will be computed from a training and a test samples, but this time using 10-fold cross validation, that is one tenth of the observations are used for the test sample and the remaining observations are used for the model estimations and the procedure is repeated and averaged ten times.

We first evaluate the prediction performance for the terminal event with the Cox models with age only and with age, AF type and diabetes, the Aalen model with age and the random survival forests with age. The reference model is taken as the Kaplan-Meier estimator and the score is computed using formula [START_REF] Ghosh | Semiparametric analysis of recurrent events data in the presence of dependent censoring[END_REF], where the prediction criterion MSE is computed from the Kaplan-Meier estimator of the censoring variable. The random survival forests were implemented from the rfsrc package (see [START_REF] Ishwaran | Random survival forests[END_REF]). The results are presented in Table 4 and Figure 5. We clearly see that the survival random forests perform poorly, especially before time 1 000 where the Kaplan-Meier shows a better performance. The Cox model with age, AF type and diabetes shows a better performance for all time points and the Aalen and Cox models with age show very similar performances and outperform all four models. Other models were also investigated with the different combinations of all three variables with each algorithm and the results were similar and are therefore omitted. In the following, we now decide to use the Cox model with age for the modelling of the terminal event. The prediction performance for the recurrent event process is now investigated. We consider the following models: four multivariate Cox models based on the age, diabetes and AF type variables, the stratified Cox model with respect to AF type with the variables age and diabetes, the Aalen model with all three variables, the Cox multi-state model with the three variables that was used in [START_REF] Schroder | Recurrent event survival analysis predicts future risk of hospitalization in patients with paroxysmal and persistent atrial fibrillation[END_REF] and the same Cox multi-state model but stratified with respect to AF type. The reference model is taken as the non-parametric estimator (see Equation ( 14)) and the score is again computed using formula [START_REF] Ghosh | Semiparametric analysis of recurrent events data in the presence of dependent censoring[END_REF]. The results are displayed in Table 5 and Figure 6 (in the figure only five different models are represented). We observe that the Cox model with the age variable has a poor predictive performance. From this model, adding the diabetes or AF type improves the model, with a much bigger gain with the AF type variable. Further, combining all three variables in the same model provides a substantial gain with respect to all previous models. On the other hand, the two multi-state models provide only a minor improvement of the predictions with a slight advantage for the stratified model. Finally, the predictions for some of these models on the expected cumulative number of AF hospitalisations are displayed in Figure 7. In this figure, the predictions are made for two 60 year old patients with diabetes, one with persistent AF and the other with paroxysmal AF. While the different models does not vary much in their predictions for the paroxysmal AF patient, they offer different results for the persistent AF patient. According to the results from Table 5 and Figure 6, the Cox multi-state model stratified with respect to AF type has the greatest prediction performance and therefore should be chosen. After 1 500 days after AF diagnosis, this model predicts an expected number of AF hospitalisations equal to 1.03 approximately. On the other hand, if one uses the multi-state Cox model, the prediction is equal to 0.98, if one uses the multivariate Cox model with all three variables, the prediction is equal to 0.76. With the Kaplan-Meier estimator as the reference, four different models are compared. For ease of visualisation, we describe those models in increasing order of their scores at time 14)) as the reference, six different models are compared. All the models use the Cox model with age as the unique covariate for the estimation of the survival function. For ease of visualisation, we describe the six models (for the recurrent events) in increasing order of their scores at time t = 2000: the Cox Figure 7: Expected cumulative number of recurrent events predicted from four different models for two 60 year old patients with diabetes, one with persistent AF and the other with paroxysmal AF. Those models had similar prediction performances (see Table 5 and Figure 6).

Conclusion

In this work a new prediction criterion was proposed in the context of recurrent event data. The criterion evaluates the prediction performance of the expected cumulative number of recurrent events, while taking into account censoring and a possible terminal event. We showed that it can be decomposed into an inseparability and imprecision terms in the same manner as in [START_REF] Graf | Assessment and comparison of prognostic classification schemes for survival data[END_REF]. Moreover, the simulations revealed that the inseparability term was largely dominant in the decomposition. As a result, we recommend to use the prediction score defined in Equation [START_REF] Ghosh | Semiparametric analysis of recurrent events data in the presence of dependent censoring[END_REF], as the difference between the prediction criterion of a given model and of a reference model, typically a model that does make use of the covariates, such that the score provides the absolute gain from the covariates in the proposed model. An alternative score could be derived by computing the relative gain as proposed in [START_REF] Steyerberg | Assessing the performance of prediction models: a framework for some traditional and novel measures[END_REF]. This produces a score that ranges from 0% to 100% and shares similarities with the Pearson's R 2 statistic. However, care should be taken with such a score, due to the fact that we normalise with respect to the prediction criterion of the reference model, which itself can be decomposed into imprecision and inseparability. This criterion could therefore be misleading due to the magnitude of the inseparability term which is unknown in practice.

The proposed prediction criterion is simple to compute and has the advantage to include all the recurrent events. As a result, it can be seen as an overall performance measure that provides information about the global predictive ability of the proposed model in a recurrent event context. Nevertheless, it would be possible to modify the criterion if one is interested into evaluating the performance of a model to only predict further recurrent events after a fixed number of events have already been experienced by a patient. This would amount to condition on a given number of experienced recurrent events in a multi-state framework. This type of criterion would be similar to the one developed in [START_REF] Schoop | Measures of prediction error for survival data with longitudinal covariates[END_REF] which conditions on being alive up to a time t * and evaluate the prediction of the model for a time s > t * . Another improvement would be to allow for frailty models in the manner of [START_REF] Van Oirbeek | Exploring the clustering effect of the frailty survival model by means of the brier score[END_REF]. A marginal score that integrates the frailty variable could be derived. Such a score would provide an overall evaluation of the frailty model and would be a natural extension of the score proposed in this paper. Alternatively, a conditional score could be proposed for the conditional (with respect to the frailty) expected cumulative number of recurrent events. More work is needed to develop these two scores.

8 Appendix: proofs of the convergence of the prediction criterion for the expected cumulative number of recurrent events under the two scenarios

In the proof of Proposition 1 we will use a different decomposition of the observed recurrent event process N , depending on the scenario under study. With a slight abuse of notations we use the same name, M , for the residual process in both scenarios.

• In the case of right-censoring only (Section 2.2) we have:

N (t) - t 0 I(C ≥ u)λ * (u | X(u))du = M (t), with E[dM (t) | {I(C ≥ u), X(u) : 0 ≤ u}] = 0.
• In the case of right-censoring and terminal event (Section 2.3) we have:

N (t) - t 0 I(T ≥ u)λ * (u | X(u))du = M (t), with E[dM (t) | {I(T ≥ u), X(u) : 0 ≤ u}] = 0.

Proof of Proposition 1

In both scenarios, we directly have:

MSE(t, µ) = E µ(t | X(t)) -µ * (t | X(t)) 2 + E   t 0 dN (u) 1 -G(u-| X(u)) -µ * (t | X(t)) 2   + 2E t 0 dN (u) 1 -G(u-| X(u)) -µ * (t | X(t)) µ * (t | X(t)) -µ(t | X(t)) .
Using the fact that E[

t 0 dN (u)/(1 -G(u-| X(u))) | X(t)] = µ * (t | X(t)), we conclude that MSE(t, µ) = E (µ(t | X(t)) -µ * (t | X(t))) 2 + A(t),
where

A(t) = E t 0 dN (u) 1 -G(u-| X(u)) 2 -E (µ * (t | X(t))) 2 . ( 15
)
We now show, in the case of right-censoring only, that A(t) ≥ 0. Using the decomposition given at the beginning of Section 8, we have,

E t 0 dN (u) 1 -G(u-| X(u)) 2 = E t 0 dM (u) 1 -G(u-| X(u)) 2 + E t 0 λ * (u | X(u))I(C ≥ u)du 1 -G(u-| X(u)) 2 , since E t 0 dM (u) 1 -G(u-| X(u)) t 0 λ * (v | X(v))I(C ≥ v)dv 1 -G(v-| X(v)) = E t 0 E[dM (u) | {I(C ≥ v), X(v) : 0 ≤ v} 1 -G(u-| X(u)) t 0 λ * (v | X(v))I(C ≥ v)dv 1 -G(v-| X(v)) = 0.
Then,

E t 0 λ * (u | X(u))I(C ≥ u)du 1 -G(u-| X(u)) 2 = 2E 0≤u<v≤t λ * (u | X(u))λ * (v | X(v)) 1 -G(u-| X(u)) dudv E (µ * (t | X(t))) 2 = 2E 0≤u<v≤t λ * (u | X(u))λ * (v | X(v))dudv .
As a result,

A(t) = E t 0 dM (u) 1 -G(u-| X(u)) 2 + 2E 0≤u<v≤t λ * (u | X(u))λ * (v | X(v))G(u-| X(u)) 1 -G(u-| X(u)) dudv ≥ 0.
The proof in the presence of a terminal event is similar. We find

A(t) = E t 0 dM (u) 1 -G(u-| X(u)) 2 + 2E 0≤u<v≤t 1 -S(u | X(u))(1 -G(u-| X(u))) 1 -G(u-| X(u)) × λ * (u | X(u))λ * (v | X(v))S(v | X(v))dudv ≥ 0.
and we can deal with all four terms in the same fashion. For instance, for the first term,

0≤u,v≤t Xu,v Ĝ(u-| x) -G(u-| x) E[dN (u)dN (v) | X(u) = x, X(v) = y] (1 -Ĝ(u-| x))(1 -Ĝ(v-| y))(1 -G(u-| x))(1 -G(v-| y)) dF X(u),X(v) (x, y) ≤ t 0 Xu Ĝ(u-| x) -G(u-| x) E[dN (u) | X(u) = x] (1 -Ĝ(u-| x))(1 -G(u-| x)) dF X(u) (x), using the fact that t 0 dN (v)/((1 -Ĝ(v-| y))(1 -G(v-| y))) is bounded. Then, since t 0 E[dN (u)/(1 -G(u-| X(u))) | X(u) = x] = µ * (t | X(t)) and (1 -Ĝ(u-| x)
) -1 is asymptotically uniformly bounded, we conclude that |C 1 (t)| tends toward 0 in probability using the uniform consistency of the censoring estimator.

For C 2 (t) we use the consistency of µ and the fact that E[µ * (t | X(t))] is finite to prove that |C 2 (t)| tends towards 0 in probability.

For C 3 (t), we directly write (

µ(t | x)) 2 -(µ(t | x)) 2 = ( µ(t | x) -µ(t | x))( µ(t | x) + µ(t | x)
) and we use the fact that µ(t | x) is bounded and the consistency of µ to prove that |C 3 (t)| tends towards 0 in probability.

Similarly to C(t) we obtain the following bound:

D(t) ≤ |D 1 (t)| + |D 2 (t)| + |D 3 (t)| with D 1 (t) = 1 n n i=1 0≤u,v≤t dN i (u)dN i (v) (1 -Ĝ(u-| X i (u)))(1 -Ĝ(v-| X i (v))) -ξ 1 (t), D 2 (t) = - 2 n n i=1 t 0 dN i (u) 1 -Ĝ(u-| X i (u)) µ(t | X i (t)) -ξ 2 (t), D 3 (t) = 1 n n i=1 µ(t | X i (t)) 2 -ξ 3 (t).
We now use the bound

|D 1 (t)| ≤ |D 1,1 (t)| + |D 1,2 (t)| + |D 1,3 (t)| with D 1,1 (t) = 1 n n i=1 0≤u,v≤t dN i (u)dN i (v) (1 -G(u-| X i (u)))(1 -G(v-| X i (v))) - 0≤u,v≤t Xu,v E[dN (u)dN (v) | X(u) = x, X(v) = y] (1 -G(u-| x))(1 -G(v-| y)) dF X(u),X(v) (x, y), D 1,2 (t) = 1 n n i=1 0≤u,v≤t χ(u, v, X i (u), X i (v))dN i (u)dN i (v) D 1,3 (t) = - 0≤u,v≤t Xu,v χ(u, v, x, y)E[dN (u)dN (v) | X(u) = x, X(v) = y]dF X(u),X(v) (x, y) and χ(u, v, x, y) = Ĝ(u-| x) -G(u-| x) + Ĝ(v-| y) -G(v-| y) + Ĝ(v-| y) -G(v-| y) + G(u-| x) G(v-| y) -Ĝ(v-| y) + Ĝ(v-| x) G(u-| y) -Ĝ(u-| y) × 1 (1 -Ĝ(u-| x))(1 -Ĝ(v-| y))(1 -G(u-| x))(1 -G(v-| y)) •
The term |D

D 2,1 (t) = - 2 n n i=1 t 0 dN i (u) 1 -G(u-| X i (u)) µ(t | X i (t)) + 2 Xt µ(t | x)µ * (t | x)dF X(t) (x),
| X i (u)) -Ĝ(u-| X i (u)))dN i (u) (1 -G(u-| X i (u)))(1 -Ĝ(u-| X i (u))) µ(t | X i (t)).
The D 2,1 (t) term converges towards 0 in probability from the law of large numbers. For D 2,2 (t), D 2,3 (t) and D 2,4 (t) we use the consistency of µ, the convergence in probability of i t 0 dN i (u))(1 -G(u-| X i (u)))/n, the boundedness of E[µ * (t | X(t))], the uniform consistency of Ĝ and the asymptotic boundedness of µ and (1 -Ĝ(u-| x)) -1 to prove that all three terms converge towards 0 in probability.

Finally, for D 3 (t), we write Each of the three terms converges towards 0 in probability using the law of large numbers for the first term and the uniform consistency of µ for the other two.

D 3 (t) = 1 n n i=1 µ(t | X i (t)) 2 - Xt µ(t | x) 2 dF X(t) (x) + 1 n n i=1 µ(t | X i (t)) 2 -µ(t | X i (t))

Proof of Proposition 3

First, note that the Brier score can be written in the following way:

MSE Brier (t, π) = E[S(t | X)] -2E[S(t | X)π(t | X)] + E[(π(t | X)) 2 ].
We now study, our prediction score MSE (t, π). Using standard martingale properties (see for instance Now, using Equation ( 16), we can rewrite B(t) in the following way:

B(t) = -E t 0 S(u | X)λ * (u | X)du + E t 0 S(u | X) 1 -G(u-) λ * (u | X)du = E t 0 G(u-) 1 -G(u-) S(u | X)λ * (u | X)du .
This shows that B(t) ≥ 0 and that this quantity does not depend on π.

Computation of A(t) in the simulation section

While the inseparability term cannot be computed on real data, it is possible to obtain its expression when the distribution of all variables are known. In this section, we provide the explicit expression of the inseparability term A(t), in the simulation setting of Section 5.1.

To this end, we use Equation [START_REF] Graf | Assessment and comparison of prognostic classification schemes for survival data[END_REF]. We first notice that, for u < v, dN

A 2 (t) = t 0 E[dN (u)] (1 -G(u-)) 2 • ( (17)
) 18
We now compute

v 0 λ 0 (u) 1 -G(u-) du = 2 × γ β 2 v 0 u γ -u du = 2 × γ β 2 × γ log γ γ -v -v ,
where we replaced λ 0 by the hazard of a Weibull distribution with shape parameter α = 2, scale parameter β and the last equality was obtained from the change of variables w = γ -u. We then need to compute the following integral in A 1 (t):

t 0 γ log γ γ -v -v λ 0 (v)dv = 2 β 2 t 0 γ log γ γ -v -v vdv = 2 β 2 γt 2 2 log(γ) -γ t 0 v log(γ -v)dv - t 3 3 .
The last integral is computed using integration by parts and then by using the change of variables w = γ -v:

t 0 v log(γ -v)dv = t 0 v 2 2 dv γ -v + t
A 2 (t) = t 0 λ 0 (t) 1 -G(t-) dtE[exp(θ 0 X)] = 2 × γ β 2 × γ log γ γ -t -t E[exp(θ 0 X)].
Finally, according to Equation [START_REF] Graf | Assessment and comparison of prognostic classification schemes for survival data[END_REF], we need to compute A 3 (t) = E (µ * (t | X(t))) 2 . From Equation ( 12), we directly have

A 3 (t) = t β 2α E[exp(2θ 0 X i )].
To conclude, A(t) = A 1 (t) + A 2 (t) -A 3 (t) and the terms involved in this equation including E[exp(θ 0 X i )] or E[exp(2θ 0 X i )] can easily be computed using Monte-Carlo simulations.

Figure 1 :

 1 Figure1: Decomposition of the MSE (solid line) in Proposition 1 as the sum between the the inseparability term A(t) (dotted line) and the imprecision term (dashed line). The data were simulated from a Cox model with two covariates and the expected cumulative number of recurrent events was predicted using the Nelson-Aalen estimator. The train sample (ntrain = 200) is used for the computation of the Nelson-Aalen estimator, the test sample (ntest = 1, 000) is used for the computation of the MSE.

Figure 2 :

 2 Figure 2: Prediction scores (see Equation (13)) using four different models. The data were generated from a Cox model with two covariates and the expected cumulative number recurrent of events was predicted using the Cox model with one covariate, the Cox model with two covariates, the Aalen model with one covariate, the Aalen model with two covariates, respectively. The reference model uses no covariates and was estimated from the Nelson-Aalen estimator. The prediction scores are computed for 100 training samples of size ntrain = 50 and a unique test sample of size ntest = 1, 000.

  One cov. for recurrent events One cov. for terminal event Two cov. for recurrent events One cov. for terminal event Two cov. for recurrent events Two cov. for terminal

Figure 3 :

 3 Figure3: Prediction scores for the recurrent events (top panels) and the survival function (bottom panels) using different models. The data were generated from a Cox model with two covariates (n = 200) for the recurrent event process and with the same two covariates for the terminal event. The expected cumulative number recurrent of events and the survival function of the terminal event were predicted using the Cox model with one or two covariates. The reference model uses no covariates and was estimated from the non-parametric estimator in Equation (14) and from the Kaplan-Meier estimator in the top and bottom panels, respectively. The prediction scores are computed for 100 training samples of size ntrain = 800 and a unique test sample of size ntest = 1, 000.

×10

  Nelson-Aalen estimator (terminal event treated as censoring) Two cov. for recurrent and terminal events (terminal event treated as censoring)

Figure 4 :

 4 Figure 4: Prediction scores for the expected cumulative number of recurrent events when the terminal event is

Figure 5 :

 5 Figure 5: Prediction scores for the survival function of the terminal event in the atrial fibrillation dataset, computed using 10-fold cross validation. With the Kaplan-Meier estimator as the reference, four different models are compared. For ease of visualisation, we describe those models in increasing order of their scores at time t = 2000: the Cox model with covariates age, AF type and diabetes (score = 0.067), the Random Survival Forest (RSF) model with covariate age (score = 0.069), the Aalen model with covariate age (score = 0.094) and the Cox model with covariate age (score = 0.105).

Figure 6 :

 6 Figure 6: Prediction scores for the expected cumulative number of recurrent events in the atrial fibrillation dataset. With the non-parametric estimator (see Equation (14)) as the reference, six different models are compared. All the models use the Cox model with age as the unique covariate for the estimation of the survival function. For ease of visualisation, we describe the six models (for the recurrent events) in increasing order of their scores at time t = 2000: the Cox model with covariate age (score = 0.911), the Cox model with covariates age and diabetes (score = 1.572), the Cox model with covariate age and AF type (score = 4.340), the multi-state (MSM) Cox model with covariates age, AF type and diabetes (score = 5.241) and the Cox model with covariates age, AF type and diabetes (score = 5.431). The MSM Cox model assumes that the transition intensities from 0 event to 1 and to "one event or more" to a new event are proportionals.

D 2 , 2 ( 1 -

 221 G(u-| X i (u)) (µ(t | X i (t)) -µ(t | X i (t))) D 2,3 (t) = 2 Xt ( µ(t | x) -µ(t | x))µ * (t | x)dF X(t) (x),

2 +µ(t | x) 2 -

 22 Xt µ(t | x) 2 dF X(t) (x).

1 -= 1 - 2 -

 112 [START_REF] Andersen | Statistical models based on counting processes[END_REF]), we directly have thatE[dN (t) | X] = H(t | X)λ * (t | X)dt, where H(t | X) = P[T > t | X] = S(t | X)(1 -G(t | X)) under independent censoring and λ * is the hazard rate of T * . As a consequence, | X)λ * (u | X)du = 1 -S(t | X),(16)since S(u | X)λ * (u | X) is equal to the conditional density function of T * . Also, it is important to notice thatE t 0 dN (u) 1 -G(u-) 2 = E t 0 dN (u) (1 -G(u-)) 2 = E t 0 S(u | X) 1 -G(u-) λ * (u | X)du ,where the first equality is due to the fact that N can only jump once and thus (t 0 dN (u)/(1 -G(u-))) 2 is simply equal to ∆I(T ≤ t)/(1 -G(T -)) 2 . Now, MSE (t, π) = E 1 -G(u-) π(t | X) + E[(π(t | X)) 2 ] 2E[(1 -S(t | X))] + E 2E[S(t | X)π(t | X)] + E[(π(t | X)) 2 ] = MSE Brier (t, π) + B(t),withB(t) = -E[1 -S(t | X)] + E t 0 S(u | X) 1 -G(u-) λ * (u | X)du .

1 -

 1 (u)dN (v) = I(C ≥ v)dN * (u)dN * (v) and E[dN (u)dN (v)] = E[I(C ≥ v)E[dN * (u)dN * (v) | X, C] = E[I(C ≥ v)λ 0 (u)λ 0 (v) exp(2θ 0 X)]dudv = (1 -G(v-))λ 0 (u)λ 0 (v)E[exp(2θ 0 X)]dudv,where we used the fact thatE[dN * (u)dN * (v) | X, C] = E[dN * (u) | X]E[dN * (v) | X] = λ 0 (u)λ 0 (v) exp(2θ0X)dudv, since under our simulation scheme, N * is independent of C and dN * (u) is independent of dN * (v) conditionally on X, for u = v. Let γ = 3, such that C follows a uniform distribution on [0, γ]. For t < γ, we have: G(u-| X(u))

2 = A 1

 21 (t) + A 2 (t),whereA 1 (t) = 2 0<u<v<t λ 0 (u)λ 0 (v) 1 -G(u-) dudvE[exp(2θ 0 X)],

Table 1 :

 1 Means and standard deviations (in bracket) over 500 simulations for the prediction score of the expected number of recurrent events. Large values indicate better predictive performances.

	2.9

Table

		n train = 100			n train = 200
	t = 1	t = 2	t = 2.9	t = 1	t = 2	t = 2.9
	One cov.-one cov. 0.68 (0.05) 6.47 (0.43)	6.14 (0.83) 0.69 (0.03) 6.54 (0.24)	6.33 (0.45)
	Two cov.-one cov. 1.67 (0.07)	12 (0.92)	6.43 (2.71) 1.67 (0.04) 12.11 (0.6)	6.84 (1.77)
	Two cov.-two cov. 1.67 (0.07) 12.61 (0.59) 10.07 (1.18) 1.68 (0.04) 12.75 (0.3) 10.51 (0.59)
		n train = 400			n train = 800
	t = 1	t = 2	t = 2.9	t = 1	t = 2	t = 2.9
	One cov.-one cov. 0.69 (0.01) 6.57 (0.14)	6.42 (0.28) 0.69 (0.01) 6.59 (0.08)	6.47 (0.17)
	Two cov.-one cov. 1.68 (0.02) 12.11 (0.39) 6.92 (1.22) 1.68 (0.01) 12.15 (0.27) 7.07 (0.85)
	Two cov.-two cov. 1.68 (0.02) 12.82 (0.17) 10.69 (0.35) 1.68 (0.01) 12.84 (0.09) 10.78 (0.23)

Table 3 :

 3 Means and standard deviations (in bracket) over 500 simulations for the prediction score of the survival function. The results were multiplied by 10 4 and are presented for the model with one or two covariates. The reference model uses no covariates and is estimated from the Kaplan-Meier estimator. Large values indicate better predictive performance.

Table 4 :

 4 t = 2000: the Cox model with covariates age, AF type and diabetes (score = 0.067), the Random Survival Forest (RSF) model with covariate age (score = 0.069), the Aalen model with covariate age (score = 0.094) and the Cox model with covariate age (score = 0.105). Means and 80% intervals (in curly bracket) over 10-folds cross validation for the prediction score of the survival function of the terminal event in the atrial fibrillation dataset. With the Kaplan-Meier estimator as the reference, four different models are compared at three different time points: the Aalen and Cox models with covariate age, the Cox model with covariates age, AF type and diabetes and the Random Survival Forest model with covariate age.

			t = 1000		t = 1500	t = 2000
	Aalen with age	0.011 [0.002, 0.017]	0.015 [0.001, 0.039]	0.038 [-0.001, 0.035]
	Cox with age	0.010 [-0.006, 0.023] 0.016 [-0.001, 0.042] 0.037 [-0.005, 0.053]
	Cox with age, AF type, diabetes 0.006 [-0.009, 0.019] 0.011 [-0.007, 0.04]	0.027 [-0.002, 0.04]
	RSF with age	0.005 [-0.021, 0.032] 0.010 [-0.014, 0.039] 0.010 [-0.031, 0.045]
			t = 1000		t = 1500	t = 2000
	Cox with age		0.036 [-0.248, 0.289] 0.089 [-0.690, 0.629] 0.283 [-0.719, 0.985]
	Cox with age and diabetes		0.166 [-0.342, 0.538] 0.232 [-0.372, 1.051] 0.532 [-0.640, 1.495]
	Cox with age and AF type		0.611 [0.104, 1.138]	1.170 [0.278, 2.081]	1.702 [-0.027, 3.989]
	Cox with age, AF type and diabetes		0.840 [-0.123, 1.539] 1.420 [-0.033, 2.983] 2.121 [-0.485, 5.683]
	Aalen with age, AF type and diabetes		0.818 [0.143, 1.312]	1.482 [0.290, 2.992]	2.140 [-0.190, 5.506]
	Cox with age, diabetes and strata(AF type)	0.836 [-0.133, 1.521] 1.416 [-0.053, 3.000] 2.118 [-0.468, 5.667]
	Cox (MSM) with age, AF type and diabetes	0.847 [-0.041, 1.459]	1.439 [0.084, 2.893]	2.136 [-0.310, 5.526]
	Cox (MSM/strata) with age, AF type and diabetes 0.852 [-0.014, 1.475]	1.454 [0.108, 3.073]	2.147 [-0.298, 5.756]

Table 5 :

 5 Means and 80% intervals (in curly bracket) over 10-folds cross validation for the prediction score of the expected number of recurrent events in the atrial fibrillation dataset. With the non-parametric estimator (see Equation (14)) as the reference, eight different models (for the recurrent events) are compared at three different time points: five Cox and Aalen models with covariates age, diabetes and AF type, one Cox model stratified with respect to AF type and two multi-state Cox models, Cox (MSM) and Cox (MSM/strata). The difference between the last two models is that the first one assumes the baseline transition intensities from 0 event to 1 and to "one event or more" to a new event to be proportionals while the last one uses two different baselines functions. All the models use the Cox model with age as the unique covariate for the estimation of the survival function.

		2.0	Models
			Cox (MSM) with age, AF type and diabetes
			Cox with age
	Prediction score for recurrent events	0.5 1.0 1.5	Cox with age and AF type Cox with age and diabetes Cox with age, diabetes and AF type
		0.0
			0	500	1000	1500	2000
					Time

  model with covariate age (score = 0.911), the Cox model with covariates age and diabetes (score = 1.572), the Cox model with covariate age and AF type (score = 4.340), the multi-state (MSM) Cox model with covariates age, AF type and diabetes (score = 5.241) and the Cox model with covariates age, AF type and diabetes (score = 5.431). The MSM Cox model assumes that the transition intensities from 0 event to 1 and to "one event or more" to a new event are proportionals.

				AF type: persistent					AF type: paroxysmal
		1.2
								Models
	Expected cumulmative number of recurrent events	0.3 0.6 0.9						Cox (MSM) with age, AF type and diabetes Cox (MSM/strata) with age, AF type and diabetes Cox with age, diabetes and AF type Cox with age, diabetes and strata(AF type)
		0.0
		0	500	1000	1500	2000	0	500	1000	1500	2000
						Time

  u,v,x,y |χ(u, v, x, y)| converges towards 0 from the uniform consistency of Ĝ while the other term converges towards a bounded quantity from the law of large numbers. The same argument applies to |D 1,3 (t)| which also converges towards 0 in probability.For D 2 (t) we write |D 2 (t)| ≤ |D 2,1 (t)| + |D 2,2 (t)| + |D 2,3 (t)| + |D 2,4 (t)| with

1,1 (t)| converges towards 0 in probability from the strong law of large numbers. The term |D 1,2 (t)| is bounded by sup u,v,x,y |χ(u, v, x, y)| 1 n n i=1 0≤u<v≤t dN i (u)dN i (v), sup

  On the other hand, computation of A 2 (t) is straightforward, using the relation (see Section 2.2)

											2
											2	log(γ -t)
						=	γ γ-t	(γ -w) 2 2w	dw +	t 2 2	log(γ -t)
						=	γ 2 2	log	γ γ -t	-γt +	γ 2 4	-	(γ -t) 2 4	+	t 2 2	log(γ -t)
						=	t 2 -γ 2 2	log(γ -t) -	γt 2	-	t 2 4	+	γ 2 2	log(γ).
	Gathering all the parts in A 1 (t), we have:
	A 1 (t)
	=	8γ β 4	γt 2 2	log(γ) -	γ 2	(t 2 -γ 2 ) log(γ -t) +	γ 2 t 2	+	γt 2 4	-	γ 3 2	log(γ) -	t 3 3	E[exp(2θ 0 X)].

E[dN (t) | X(t)] = (1 -G(t-))λ * (t | X(t))dt. For t < γ, we have

Proof of Proposition 2

We start by proving that E[µ * (τ | X(τ ))] < ∞ in the presence of a terminal event (the scenario without terminal event follows from the same arguments). We have for all t ∈ [0, τ ] :

From the same assumption, N (τ ) is almost surely bounded by a constant. As a consequence,

is almost surely bounded, where the equality has been proved in Section 2.3. The rest of the proof of Proposition 2 is identical in both scenarios.

We first note

we let X u,v denote the support of the joint distribution (X(u), X(v)) and we note

We then introduce the quantity

Write:

By decomposing the square term into three other terms, we bound C(t) in the following way:

For C 1 (t) we have