<figure><img src="image_1.png" /></figure>

HAL
open science

# Random survival forests for the analysis of recurrent events for right-censored data, with or without a terminal event

Juliette Murris, Olivier Bouaziz, Michal Jakubczak, Sandrine Katsahian, Audrey Lavenu

<figure><img src="image_2.png" /></figure>

## ► To cite this version:

Juliette Murris, Olivier Bouaziz, Michal Jakubczak, Sandrine Katsahian, Audrey Lavenu. Random survival forests for the analysis of recurrent events for right-censored data, with or without a terminal event. 2024. hal-04612431

HAL Id: hal-04612431

https://hal.science/hal-04612431

Preprint submitted on 14 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Random survival forests for the analysis of recurrent events for right-censored data, with or without a terminal event

Juliette Murris, Olivier Bouaziz, Michal Jakubczak, Sandrine Katsahian,
Audrey Lavenu

► To cite this version:

Juliette Murris, Olivier Bouaziz, Michal Jakubczak, Sandrine Katsahian, Audrey Lavenu. Random survival forests for the analysis of recurrent events for right-censored data, with or without a terminal event. 2025. hal-04612431v2

HAL Id: hal-04612431

https://hal.science/hal-04612431v2

Preprint submitted on 9 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Random survival forests for the analysis
of recurrent events for right-censored
data, with or without a terminal event

Juliette Murris¹⑤, Olivier Bouaziz²⑤, Michal Jakubczak³,
Sandrine Katsahian⁴⑤, and Audrey Lavenu⁵①*

¹HeKA, Inria, Inserm, Centre de recherche des Cordeliers, Université Paris Cité, Paris, France, & R&D, Pierre Fabre, Boulogne-Billancourt, France, e-mail: juliette.murris-ext@aphp.fr

²Université Paris Cité, CNRS, MAP5, F-75006 Paris, France, e-mail: olivier.bouaziz@parisdescartes.fr

³R&D, Pierre Fabre, Boulogne-Billancourt, France, & Ardigen S.A., Kraków, Poland, e-mail: michal.jakubczak-ext@pierre-fabre.com

⁴HeKA, Inria, Inserm, Centre de recherche des Cordeliers, Université Paris Cité, Paris, France, Unité de Recherche Clinique, Hôpital Européen Georges-Pompidou, Assistance Publique – Hôpitaux de Paris (AP-HP), Centre, Paris, France, & Centre d'Investigation Clinique 1418 Épidémiologie Clinique, Paris, France, e-mail: sandrine.katsahian@aphp.fr

⁵Faculté de Médecine, Université de Rennes, Rennes, France, Institut de Recherche Mathématique de Rennes (IRMAR), Rennes, France, & Centre de Investigation Clinique 1414, Inserm, Université de Rennes, Rennes, France, e-mail: audrey.lavenu@univ-rennes.fr

**Abstract:** Random survival forests (RSF) have emerged as valuable tools in medical research. They have shown their utility in modelling complex relationships between predictors and survival outcomes, overcoming linearity or low dimensionality assumptions. Nevertheless, RSF have not been adapted to right-censored data with recurrent events (RE). This work introduces RecForest, an extension of RSF and tailored for RE data, leveraging principles from survival analysis and ensemble learning. RecForest adapts the splitting rule to account for RE, with or without a terminal event, by employing the pseudo-score test or the Wald test derived from the marginal Ghosh-Lin model. The ensemble estimate is constructed by aggregating the expected number of events from each tree. Performance metrics involve a concordance index (C-index) tailored for RE analysis, along with an extension of the mean squared error (MSE). A comprehensive evaluation was conducted on both simulated and open-source data. We compared RecForest against the non-parametric mean cumulative function and the Ghosh-Lin model.

<footnote>*SK and AL share last authorship.</footnote>
Random survival forests for the analysis of recurrent events for right-censored data, with or without a terminal event

Juliette Murris<sup>1,2*</sup>, Olivier Bouaziz<sup>3</sup>, Michal Jakubczak<sup>2,4</sup>, Sandrine Katsahian<sup>1,5,6†</sup> and Audrey Lavenu<sup>7,8,9‡</sup>

<sup>1</sup> HeKA, Inria, Inserm, Centre de recherche des Cordeliers, Université Paris Cité, Paris, France

<sup>2</sup> R&D, Pierre Fabre, Boulogne-Billancourt, France

<sup>3</sup> Université Paris Cité, CNRS, MAP5, F-75006 Paris, France

<sup>4</sup>Ardigen S.A., Kraków, Poland

<sup>5</sup> Unité de Recherche Clinique, Hôpital Européen Georges-Pompidou, Assistance Publique – Hôpitaux de Paris (AP-HP), Centre, Paris, France

<sup>6</sup> Centre d'Investigation Clinique 1418 Épidémiologie Clinique, Paris, France

<sup>7</sup> Faculté de Médecine, Université de Rennes, Rennes, France

<sup>8</sup> Institut de Recherche Mathématique de Rennes (IRMAR), Rennes, France

<sup>9</sup> Centre de Investigation Clinique 1414, Inserm, Université de Rennes, Rennes, France

*JM is the corresponding author (email: murris.juliette@gmail.com); † SK and AL share last authorship.
Across the simulations and application, RecForest consistently outper-
forms, exhibiting C-index values ranging from 0.64 to 0.80 and lowest
MSE metrics. As analysing time-to-recurrence data is critical in medi-
cal research, the proposed method represents a valuable addition to the
analytical toolbox in this domain.

**Keywords and phrases:** Random forests, Recurrent events, Survival analyses, Terminal events, High-dimensional data.

# 1. Introduction

Recurrent events refer to instances where individuals may experience mul-
tiple occurrences of the same event over time. In medical research, patients
may face recurrent disease relapses, frequent hospitalizations, or repeated
surgeries. While traditional survival analyses focus solely on the first occur-
rence of an event, specific statistical models have been developed to capture
the complexity of recurrence in a survival framework. Intensity models rely
on instantaneous hazards at each time point and account for dependence
amongst event occurrences captured by time-varying covariates (Andersen
and Gill (1982); Prentice, Williams and Peterson (1981)). Besides, marginal
models centre on the overall distribution of event times and the cumulative
event counts (Wei, Lin and Weissfeld (1989); Cook, Lawless and Lee (2010)).
For a more in-depth exploration of these models concerning recurrent events,
comprehensive discussions can be found in works by Amorim and Cai (2015)
and Ozga, Kieser and Rauch (2018).

Time-to-event analyses are systematically challenging due to the presence
of censoring, i.e. when the precise timing of an event remains unknown or
unobserved. Above methodologies strictly assume the censoring process to
be uninformative, hence independent of the underlying event process. Never-
theless, a terminal event may occur in competition, preventing further events
of interest from happening. A terminal event is then a specific type of event
considered as a termination point for the study period, making the censor-
ing process no longer uninformative. Strategies for handling terminal events
include ignoring them, although this approach is acknowledged to be flawed,
or accounting for competing risks. Several pertinent statistical models en-
able to analyse both recurrent events and competing risks (Charles-Nelson,
Katsahian and Schramm (2019)).

Navigating medical data introduces numerous challenges, including high-
dimensionality, variable selection, and multicollinearity. To address these,
survival time-to-first-event approaches have integrated statistical and ma-
chine learning techniques. In practice, various algorithms now have their
**Abstract. Background:** Random survival forests (RSF) have emerged as valuable tools in medical research. They have shown their utility in modelling complex relationships between predictors and survival outcomes, overcoming linearity or low dimensionality assumptions. Nevertheless, RSF have not been adapted to right-censored data with recurrent events (RE). **Methods:** This work introduces RecForest, an extension of RSF and tailored for RE data, leveraging principles from survival analysis and ensemble learning. RecForest adapts the splitting rule to account for RE, with or without a terminal event, by employing the pseudo-score test or the Wald test derived from the marginal Ghosh-Lin model. The ensemble estimate is constructed by aggregating the expected number of events from each tree. Performance metrics involve a concordance index (C-index) tailored for RE analysis, along with an extension of the mean squared error (MSE). A comprehensive evaluation was conducted on both simulated and open-source data. We compared RecForest against the non-parametric mean cumulative function and the Ghosh-Lin model. **Results:** Across the simulations and application, RecForest consistently outperforms, exhibiting C-index values ranging from 0.60 to 0.82 and lowest MSE metrics. **Conclusions:** As analysing time-to-recurrence data is critical in medical research, the proposed method represents a valuable addition to the analytical toolbox in this domain. The RecForest implementation is publicly available as an R package on CRAN.

**Key words:** Random forests, Recurrent events, Survival analyses, Terminal events, High-dimensional data
survival counterparts that are effectively employed to answer medical questions in real-world applications (Huang et al. (2023)). For instance, penalized regression methods, such as LASSO (Least Absolute Shrinkage and Selection Operator), Ridge, and Elastic-Net, have been tailored for Cox models, facilitating variable selection and regularization (Cox (1972); Tibshirani (1997)). Support-vector machines introduced by Van Belle et al. (2011), renowned for their capacity to handle high-dimensional data and non-linearity, have also been extended to survival endpoints. Likewise, random survival forests (RSF) from Ishwaran et al. (2008) embody a powerful ensemble learning technique handling interactions. The RSF algorithm has been extended to model several phenomena, such as competing risks, or longitudinal data (Ishwaran et al. (2014); Devaux et al. (2023)). However, within the survival framework, no machine learning approach has hitherto been extended to recurrent events (Murris et al. (2023)). To address these unmet needs and confront the aforementioned challenges, we introduce the first RSF capable of handling recurrent events, with or without a terminal event. Illustrated in Figure 1, our method entails a 5-step approach that i) discerns the relevance of recurrent and terminal events, ii) grows trees to construct a coherent RSF, iii) thoroughly assesses performance, iv) provides relevant variable importance, and v) enables predictions on new data.

In this paper, we consider $n$ individuals. Let $N^*(t)$ be the number of recurrent events that occur in the time interval $[0, t]$, $D$ the survival time and $C$ the censoring time. The data is made of $(N(\cdot), \Upsilon, \delta)$ where $N(t) = N^*(t \wedge C)$, $\Upsilon = D \wedge C$, $\delta = I(D \le C)$, where $a \wedge b = \min(a, b)$ and $I(\cdot)$ is the indicator function. For $i = 1, \dots, n$, $(N_i(\cdot), \Upsilon_i, \delta_i)$ are assumed to be independent replicates of $(N(\cdot), \Upsilon, \delta)$. The marginal mean frequency function is $\mu(t) = \mathbb{E}[N(t)]$. An estimator of $\mu$ in the absence of a terminal event is the Nelson-Aalen estimator from Lawless and Nadeau (1995), that writes

$$ (1) \qquad \hat{\mu}(t) = \hat{R}(t) = \int_{0}^{t} \frac{dN(u)}{Y(u)} $$

with $N(t) = \sum_i N_i(t)$, and $Y(t) = \sum_i Y_i(t)$ the number of individuals at risk at time $t$. In presence of a terminal event, we have $\mu(t) = \int_0^t S(u)dR(u)$ where $S(t) = \mathbb{P}(D \ge t)$ and $dR(t) = \mathbb{E}[dN^*(t)|D \ge t]$ (Cook and Lawless (1997); Ghosh and Lin (2000)). The associated estimator writes

$$ (2) \qquad \hat{\mu}(t) = \int_{0}^{t} \hat{S}(u)d\hat{R}(u) = \int_{0}^{t} \hat{S}(u) \frac{\sum_{i} Y_{i}(u)dN_{i}(u)}{\sum_{i} Y_{i}(u)} $$
# 1. Introduction

Recurrent events refer to instances where individuals may experience multiple occurrences of the same event over time. In medical research, patients may face recurrent disease relapses, frequent hospitalizations, or repeated surgeries. While traditional survival analyses focus solely on the first occurrence of an event, specific statistical models have been developed to capture the complexity of recurrence in a survival framework. Intensity models rely on instantaneous hazards at each time point and account for dependence amongst event occurrences captured by time-varying covariates (1,2). Besides, marginal models centre on the overall distribution of event times and the cumulative event counts (3,4). For a more in-depth exploration of these models concerning recurrent events, comprehensive discussions can be found in works by Amorim (2015) and Ozga (2018) (5,6).

Time-to-event analyses are systematically challenging due to the presence of censoring, i.e. when the precise timing of an event remains unknown or unobserved. Above methodologies strictly assume the censoring process to be uninformative, hence independent of the underlying event process. Nevertheless, a terminal event may occur in competition, preventing further events of interest from happening. A terminal event is then a specific type of event considered as a termination point for the study period, making the censoring process no longer uninformative. Strategies for handling terminal events include ignoring them, although this approach is acknowledged to be flawed, or accounting for competing risks. Several pertinent statistical models enable to analyse both recurrent events and competing risks (7).

Navigating medical data introduces numerous challenges, including high-dimensionality, variable selection, and multicollinearity. To address these, survival time-to-first-event approaches have integrated statistical and machine learning techniques. In practice, various algorithms now have their survival counterparts that are effectively employed to answer medical questions in real-world applications (8). For instance, penalized regression methods, such as LASSO (Least Absolute Shrinkage and Selection Operator), Ridge, and Elastic-Net, have been tailored for Cox models, facilitating variable selection and regularization (9,10). Support-vector machines, renowned for their capacity to handle high-dimensional data and non-linearity, have also been extended to survival endpoints (11). Likewise, random survival forests (RSF) embody a powerful ensemble learning technique handling interactions (12). The RSF algorithm has been extended to model several phenomena, such as competing risks, or longitudinal data (13,14). However, within the survival framework, no machine learning approach has hitherto been extended
<figure><img src="image_3.png" /><figcaption>FIG 1. Scheme of the use of RecForest for survival data with recurrent events in presence or absence of a terminal event</figcaption></figure>
<figure><img src="image_2.png" /><figcaption>Figure 1. Scheme of the use of RecForest for survival data with recurrent events in presence or absence of a terminal event</figcaption></figure>

to recurrent events (15). To address these unmet needs and confront the aforementioned challenges, we introduce the first RSF capable of handling recurrent events, with or without a terminal event. Illustrated in Figure 1, our method entails a 5-step approach that i) discerns the relevance of recurrent and terminal events, ii) grows trees to construct a coherent RSF, iii) thoroughly assesses performance, iv) provides relevant variable importance, and v) enables predictions on new data.

In this paper, we consider n individuals. Let $T_{ij}$ be the time of the j-th recurrent event for subject i = 1,...,n, and $N_i^*(t)$ denote the true recurrent event counting process up to t, defined by

$$N_i^*(t) = \sum_{j=1}^{\infty} 1(T_{ij} \le t).$$

Let $D_i$ denote the time of a terminal event (e.g., death), with $N_i^*(t) = N_i^*(D_i)$ for all $t \ge D_i$, and let $C_i$ be the independent right-censoring time. Define $\gamma_i = \min(D_i, C_i)$, and let $\delta_i = I(D_i < C_i)$ be the terminal event indicator. The observed counting process is then
And $\hat{S}(t)$ is the Kaplan-Meier estimator of $S(t)$ based on $(\Upsilon_i, \delta_i)$ (Kaplan and Meier (1958)).

From the above considerations arises the evaluation of the provided estimations. Within survival framework, a widely common metric is an extension of the area under the ROC curve known as the concordance index (C-index). The principle of the C-index from Harrell, Lee and Mark (1996) and its derivatives as per Uno et al. (2011) is to measure the ability of a model to correctly order pairs of survival times. Recent developments from Kim, Schaubel and McCullough (2018) have expanded the application of the C-index to the recurrent event framework, incorporating the number of subsequent event occurrences. However, the number of events over time is only comparable if individuals have similar follow-up, which is hardly the case in real-world settings. Therefore, we proposed a generalized C-index by introducing event occurrence rate. Additionally, we employ the mean-square error, recently adapted to account for recurrent events by Bouaziz (2024).

Based on the non-parametric estimators and ensemble method principles, the objective of this work is to introduce a new ensemble approach, called RecForest, for the analysis of recurrent events in a survival framework, with or without a terminal event. The overall methodology based on survival decision trees and novel associated evaluation metrics is described in Section 2. Section 3 displays an extended simulation scheme for the comprehension of the proposed methodology. Illustrative examples based on open-source data are used for concrete application in Section 4.

## 2. Methodology

The proposed Algorithm 1 is an extension of the RSF introduced by Ishwaran et al. (2008). The first step is drawing bootstrap samples to prevent overfitting and capture inherent variability within the original dataset. Then, survival trees are constructed on each bootstrap sample. Unlike the original RSF, our approach accommodates for subsequent events by integrating statistical considerations tailored for recurrent events analysis. As a last step, the algorithm aggregates the results over the constructed recursive survival trees to obtain a comprehensive estimate.

Next subsections describe in further details how survival trees grow for constructing the random forest. Additionally, we provide adequate metrics for the evaluation. Finally, we expound on the computation of variable importance.
$$N_i(t) = N_i^*(\min(t, \gamma_i)),$$

and the observed data consists of $(N_i(t), \gamma_i, \delta_i, X_i)$ for $i = 1, \dots, n$, where $X_i$ are covariates.

The marginal mean frequency function is defined $\mu(t) = E[N^*(t)]$, representing the true cumulative number of recurrent events up to time $t$. To estimate $\mu(t)$, we adopt a unified estimator that accounts for both right-censoring and, when present, terminal events. Let $Y_i(t) = 1(\gamma_i \ge t)$ denote the at-risk indicator for subject $i$ at time $t$, and let $\hat{S}(t)$ denote the Kaplan-Meier estimator of the survival function of the terminal event $D$. The general form of the estimator is:

$$\hat{\mu}(t) = \int_{0}^{t} \left( \frac{1}{n} \sum_{i=1}^{n} \frac{Y_i(u)}{\hat{S}(u)} \right)^{-1} \left( \frac{1}{n} \sum_{i=1}^{n} \frac{Y_i(u)}{\hat{S}(u)} \right) dN_i(u),$$

or more compactly:

$$\hat{\mu}(t) = \int_{0}^{t} \frac{\sum_{i=1}^{n} \frac{Y_i(u)}{\hat{S}(u)} dN_i(u)}{\sum_{i=1}^{n} \frac{Y_i(u)}{\hat{S}(u)}}.$$

This inverse probability-weighted estimator accounts for the informative censoring induced by a terminal event (16). In the absence of a terminal event, the survival probability simplifies to $\hat{S}(t) = 1$ for all $t$, and the estimator reduces to the Nelson-Aalen estimator for the cumulative mean function (17):

$$\hat{\mu}(t) = \int_{0}^{t} \frac{\sum_{i=1}^{n} Y_i(u) dN_i(u)}{\sum_{i=1}^{n} Y_i(u)}.$$

From the above considerations arises the evaluation of the provided estimations. Within survival framework, a widely common metric is an extension of the area under the ROC curve known as the concordance index (C-index). The principle of the C-index and its derivatives is to measure the ability of a model to correctly order pairs of survival times (18,19). Recent developments have expanded the application of the C-index to the recurrent event framework, incorporating the number of subsequent event occurrences (20). However, the number of events over time is only comparable if individuals have similar follow-up, which is hardly the case in real-world settings. Therefore, we proposed a generalized C-index
<caption>**Algorithm 1 Overview of RecForest algorithm**</caption>

**Require:** Draw $B > 0$ bootstrap samples from the learning data
**for** Each node of survival tree **b** **do**
>mtry predictors are randomly selected with mtry ∈ N, mtry ≤ p;
A greedy algorithm for optimal threshold research is used to maximize the test statistic;
The tree grows until the stopping rule is met based on the minimal number of events
*minsplit* and the minimal number of individuals in terminal nodes *nodesize*;
Estimate $\hat{\mu}_b$ is computed;
**end for**
Estimate *M* is computed over the *B* trees.

## 2.1. Growing trees with recurrent events

### 2.1.1. Splitting rules

At each node $h \in \mathcal{H}$, the ongoing subsample is split into two daughter nodes denoted $h^{(+)}$ and $h^{(-)}$. The aim of the split is to make the daughter nodes as different as possible with regards to the outcome. The splitting rule requires that each of the *mtry* randomly drawn variable is dichotomized. For continuous variables, random split points, quartiles, and deciles are considered. Let $\mathbf{x}_h = \{A, B\}$ be the dichotomized vector of a variable inherited from $h$. With no terminal event, we compare the marginal mean frequency functions $\mu_A(t)$ and $\mu_B(t)$. The null hypothesis is their equality. In absence of a terminal event, we use the two-sample test akin to the log-rank test from Lawless and Nadeau (1995). The test statistic writes $U(t) = \int_0^t \frac{Y_A(u)Y_B(u)}{Y_A(u)+Y_B(u)}(d\hat{\mu}_A(u) - d\hat{\mu}_B(u))$. To incorporate the presence of a terminal event, we employ the marginal Gosh-Lin (GL) model from Ghosh and Lin (2002) within the single variable $\mathbf{x}_h$. Acknowledging there are no further recurrence after the terminal event, the marginal mean up to $t$ associated with $\mathbf{x}_h$ is defined as $\mu_{\mathbf{x}_h}(t) = E[N^*(t)|\mathbf{x}_h] = \mu_0(t) \times \exp(\beta \mathbf{x}_h)$ with $\mu_0$ left unspecified and $\beta$ the regression coefficient. To accommodate longitudinal variables, the GL model considers a rate function $d\mu_{\mathbf{x}_h}(t) = d\mu_0(t) \times \exp(\beta \mathbf{x}_h(t))$. The Wald test statistic is then extracted from $\mu_{\mathbf{x}_h}$ and $d\mu_{\mathbf{x}_h}$ to test the null hypothesis of $\beta = 0$.

The variable selected for node $h$ is the one that maximizes the adequate test statistic to generate $h^{(+)}$ and $h^{(-)}$, based on the presence of a terminal event and/or longitudinal variables.
by introducing event occurrence rate. Additionally, we employ the mean-square error, recently adapted to account for recurrent events (21).

Based on non-parametric estimators and ensemble method principles, our contribution is to introduce a new ensemble approach, called RecForest, for the analysis of recurrent events in a survival framework, with or without a terminal event. A review of existing methods and related work is provided in Section 2 to contextualize our work. The overall methodology, based on survival decision trees and novel associated evaluation metrics, is detailed in Section 3. Section 4 displays an extended simulation study to evaluate the proposed methodology, and Section 5 provides an illustrative medical application using open-source data.

# 2. Background and related work

## 2.1. Survival analysis for recurrent events

Survival analysis traditionally focuses on modeling the time until a single event occurs, such as death or disease recurrence. However, many real-world phenomena involve multiple occurrences of the same event over time, necessitating specialized methods for handling recurrent event data. Recurrent event processes arise in various fields, particularly in biomedical research, where individuals may experience repeated hospitalizations, disease relapses, or treatment cycles (22–24).

### 2.1.1. Counting processes and event history

Recurrent event data can be represented using counting processes, where $N_i(t)$ counts the number of events for subject $i$ up to time $t$. The event history of a process, also called filtration, captures all past events up to time $t$ (22). This historical perspective is fundamental in influencing model selection between conditional and marginal approaches. Conditional models account for past events through covariates, whereas marginal models assume independence between occurrences.

### 2.1.2. Recurrent event models

Recurrent event models can be broadly categorized based on their underlying assumptions about event dependence and the risk structure. The Andersen-Gill model extends the Cox proportional hazards model by considering the counting process representation of recurrent events (25). This approach assumes that each event is independent, conditional on observed covariates, making it suitable for scenarios where the dependence between recurrent events is explained entirely by time-dependent covariates. In contrast, the Prentice-Williams-Peterson model introduces a stratified approach where each event is treated as a separate
### 2.1.2. Terminal node estimator

Let b be a bootstrap sample from original data on which a tree is grown and x a p-dimensional vector of covariates dropped down the tree. The node-specific event count $N_b(t|x)$ is the number of recurrent events before censoring or a terminal event at time t. The associated number of individuals at risk $Y_b(t|x)$ is the number of individuals that were not censored, or that did not encounter a terminal event by time t. We then define a tree-specific estimate as follows

$$ (3) \qquad \hat{\mu}_b(t|\mathbf{x}) = \hat{R}_b(t|\mathbf{x}) = \int_0^t \frac{N_b(du|\mathbf{x})}{Y_b(du|\mathbf{x})} $$

In case of the presence of a terminal event,

$$ (4) \qquad \hat{\mu}_b(t|\mathbf{x}) = \int_0^t \hat{S}_b(u|\mathbf{x}) d\hat{R}_b(u|\mathbf{x}) $$

Individuals from the same terminal node share similar features inherited from their tree path, along with identical estimates. As per the splitting rule, the terminal node estimator depends on the presence of a terminal event in the original sample.

### 2.1.3. Pruning trees

A pruning strategy is essential to help find a trade-off to prevent overfitting and improve generalization performance of trees, within a reasonable computational time. Aligned with Devaux et al. (2023), we suggest two stopping rules for each terminal node: (i) a minimal number of events called *minsplit*, and (ii) a minimal number of individuals called *nodesize*. The validation of either stopping rule designates the current node h as terminal.

### 2.1.4. Handling missing data

To tackle eventual missing data, we include an adaptive-tree imputation which addresses missing data during the tree-growing stage by selectively drawing from available, non-missing, in-bag data (Ishwaran et al. (2008); Chen and Xu (2023)). At each node $h_b$ from tree b, the method entails imputing random non-missing information specifically from the selected variables. The imputed data is then utilized for making splits within the node $h_b$. Imputed values are reset to missing as the tree progresses to subsequent nodes.
stratum, allowing for varying baseline hazards across different recurrences (26). This stratification accounts for event order, making the model more suitable when the hazard function is expected to change with each recurrence.

An alternative perspective is provided by frailty models, which introduce an unobserved random effect to account for unmeasured heterogeneity across individuals (27). These models assume that individuals with a high frailty term are more likely to experience recurrent events, capturing dependencies that are not explained by observed covariates. The frailty component is typically modeled using a gamma distribution, leading to a shared frailty model where the recurrent events of an individual are correlated through the random effect. Such models are particularly useful in medical studies where genetic predisposition or underlying susceptibility factors contribute to repeated event occurrences (5).

Marginal models provide another avenue for analyzing recurrent events by focusing on the mean number of occurrences rather than the hazard function. The Wei, Lin, and Weissfeld model treats each recurrence as an independent process, estimating covariate effects without making assumptions about within-subject correlation (28). This model is particularly useful when the goal is to assess the overall effect of treatment or exposure on event frequency rather than the timing of each recurrence. The mean cumulative function approach further complements marginal models by estimating the expected cumulative number of events over time, providing an intuitive interpretation of recurrence patterns in a population (17).

The choice of model depends on the research question and data structure (5,7,29). Conditional models like Andersen-Gill are well-suited for time-to-event analyses where covariates explain recurrence dependencies, while marginal models such as Wei-Lin-Weissfeld are preferred when assessing cumulative event counts. Frailty models bridge these approaches by introducing individual-level heterogeneity, making them valuable in scenarios where unobserved factors influence event recurrence. The application of these models spans medical research, reliability analysis, and epidemiology, ensuring that recurrent event processes are appropriately accounted for in survival studies.

### 2.1.3. With terminal events

The presence of a terminal event, such as death, introduces a competing risk scenario where further recurrences become impossible (30,31). Traditional recurrent event models assume non-informative censoring, but when a terminal event is present, this assumption is violated, leading to biased estimates if
## 2.2. From trees to random forests

### 2.2.1. Ensemble estimates

Once all B trees are grown from the independent bootstrap samples, the ensemble estimate $\hat{M}$ is the aggregation of all B tree-specific estimates. We define $\hat{M}$ as

$$ (5) \qquad \hat{M}(t|\mathbf{x}) = \frac{1}{B} \sum_{b=1}^{B} \hat{\mu}_b(t|\mathbf{x}) $$

### 2.2.2. OOB ensemble estimates

By standard bootstrap theory, each bootstrap sample leaves out circa 37% of the data (Ishwaran et al. (2008)). This is the so-called out-of-bag (OOB) sample. OOB data is used to build OOB ensembles. Let $O_i \subseteq \{1, \dots, B\}$ be the index set of trees where for $b \in O_i$, $c_{i,b} = 0$, which means the individual i is in the OOB sample. The OOB ensemble estimate $\hat{M}^{OOB}$ of aggregated tree-specific estimates for i which is OOB writes

$$ (6) \qquad \hat{M}^{OOB}(t|\mathbf{x}_i) = \frac{1}{|O_i|} \sum_{b \in O_i} \hat{\mu}_b(t|\mathbf{x}_i) $$

OOB ensemble estimates are typically used for reporting errors.

## 2.3. Performance

Performance metrics below indicate the ability of the model to predict well from training data to unseen data. In our case, unseen data are either from the OOB sample or external validation data.

### 2.3.1. Assessing performance with relevant metrics

For the assessment of performance, we introduce an extended version of the C-index and employ the mean-square error (MSE), a derived score, and their integrated versions (Figure 2).

*Concordance index.* Kim, Schaubel and McCullough (2018) adapted the C-index to recurrent events and considered the number of events over time across individuals. This metric hence suffers from the potential bias in case of
not appropriately addressed. To correct for this, an adjusted mean cumulative function estimator has been proposed, which incorporates survival probabilities (16,17).

Other approaches for handling terminal events include multi-state models, which explicitly define transitions between recurrent and terminal states, and joint modeling techniques that simultaneously estimate survival and recurrent event processes while accounting for their dependence (32).

## 2.2. Random survival forest algorithm

Random Survival Forests (RSF) were introduced as an extension of classical random forests to survival data (12). Unlike traditional parametric models, RSF does not assume a specific form for the survival function, providing greater flexibility across different data distributions.

### 2.2.1. Principle of Random Survival Forests

RSF is based on ensemble learning, where multiple survival trees are built from bootstrap samples of the training data. Each tree is constructed via recursive partitioning using a splitting criterion adapted to censored data. The two main splitting rules commonly used are the log-rank splitting, and the likelihood ratio statistics.

Like classical random forests, RSF leverages bootstrap sampling to generate multiple trees. Observations not included in a given bootstrap sample (out-of-bag or OOB) are used to obtain an internal estimation of the model's error, eliminating the need for a separate validation set.

### 2.2.2. Terminal node estimator and prediction

In RSF, survival is estimated within each terminal node of the tree. Unlike standard decision trees that assign a class label to a terminal node, RSF constructs a non-parametric Nelson-Aalen estimator for the cumulative hazard function. The survival function for an individual in a terminal node is derived as:

$$ \hat{S}(t) = \exp\left(-\sum_{u \le t} \hat{H}(u)\right), $$

where $\hat{H}(u)$ is the cumulative hazard function estimated using the individuals within that node. The final survival prediction for a new observation is obtained by averaging the survival estimates across all trees in the forest. This ensemble approach helps smooth individual tree predictions and improves robustness.
<figure><img src="image_4.png" /><figcaption>FIG 2. Illustration of the performance metrics with true and predicted cumulative number of events over time</figcaption></figure>

substantial variability in the follow-up times. Individuals with longer follow-up times or a higher number of events might indeed disproportionately influence the C-index calculation. To address this issue, we suggest using occurrence rates by computing event rates per unit time.

The proposed C-index is defined as the proportion of all concordant pairs of individuals where predicted occurrence rates are correctly ordered with respect to observed occurrence rates (as shown in Figure 2a). As occurrence rates can be calculated for all individuals, including censored ones, the proposed C-index is not partial and considers all individuals in the computation. In this work, the C-index then writes

$$ (7) \qquad \hat{C}_{\text{rec}} = \frac{\sum_{i=1}^{n} \sum_{j=1}^{n} \mathbb{1}_{r_i > r_j} \times \mathbb{1}_{\hat{r}_i > \hat{r}_j}}{\sum_{i=1}^{n} \sum_{j=1}^{n} \mathbb{1}_{r_i > r_j}} $$

with $r_i = \frac{N_i(T_i)}{T_i}$ and $\hat{r}_i = \frac{\hat{\mu}(T_i|\mathbf{x}_i)}{T_i}$ the observed and predicted event occurrence rates, respectively. Like other C-indices, the value of the above C-index falls within the range of 0 to 1, where 1 indicates perfect concordance, and values close to 0.5 suggest randomness in the model.

*Mean-squared error and derived score.* No MSE measure has been adapted to recurrent events framework until very lately. Bouaziz (2024) filled this gap and suggested a generalization of the Brier score from Graf et al. (1999). For our problematic, for each tree b, we define

$$ (8) \qquad \widehat{MSE}_b(t, \hat{\mu}_b) = \frac{1}{n} \sum_{i=1}^{n} \left( \int_{0}^{t} \frac{dN_i(u)}{\hat{G}_c(u|\mathbf{x})} - \hat{\mu}_b(t|\mathbf{x}) \right)^2 $$

Where $\hat{G}_c(u|\mathbf{x}) = 1 - \hat{G}(u - |\mathbf{x})$ is an estimator of $G_c(u|\mathbf{x}) = 1 - G(u - |\mathbf{x})$ the conditional cumulative distribution function of the censoring variable C
### 2.2.3. Extensions and recent work

Since their introduction, RSF has been extended to address various specific needs. One of the key areas of development has been the adaptation of RSF to competing risks scenarios, where multiple event types exist, and the occurrence of one event precludes the occurrence of another. Ishwaran (2014) proposed an adaptation of RSF for competing risks using Gray's test, which allows for the estimation of cumulative incidence functions in a non-parametric manner (13). This extension has significantly improved the applicability of RSF in medical and epidemiological studies where multiple event types, such as different causes of death, must be considered.

Another significant extension of RSF involves dynamic prediction models. Pickett (2021) introduced landmark approaches that enable real-time updates of survival predictions as new covariate information becomes available (33). This extension is particularly useful in clinical settings where patient characteristics evolve over time. More recently, Devaux (2023) combined dynamic prediction techniques with longitudinal data (14). This framework enhances RSF's ability to incorporate time-dependent covariates, making it a powerful tool for personalized medicine and disease progression modeling.

Bayesian approaches have also been integrated into RSF to improve flexibility and incorporate prior information into survival analysis. Chipman (2010) introduced Bayesian Additive Regression Trees (BART), which were later adapted for survival data and incorporated into RSF (34). These Bayesian RSF models allow for probabilistic inference and uncertainty quantification, making them particularly useful in high-dimensional and sparse data settings. Additionally, Bayesian methods facilitate variable selection and regularization, helping to address overfitting issues in complex survival datasets.

Beyond these advances, RSF has also been extended for high-dimensional variable selection. Ishwaran (2010) introduced RSF-VH (Variable Hunting), which ranks variables based on their minimal depth within survival trees, providing a non-parametric approach to feature selection (35). Furthermore, Genuer (2010) proposed a structured variable selection process that balances interpretability and predictive accuracy (36). Another key contribution introduced regularization techniques to refine feature selection and improve the model's robustness against irrelevant or redundant variables (37).

Finally, ongoing research is focusing on the computational efficiency and scalability of RSF, particularly for large-scale biomedical datasets. Efforts are being made to parallelize RSF algorithms and optimize tree
given $\mathbf{x}$. We assume $C$ and $\mathbf{x}$ to be independent. With no terminal event, $\hat{G}$ is the empirical cumulative distribution function of the censored variable. In the presence of a terminal event, $\hat{G}$ is the Kaplan-Meier estimator of $C$. As suggested in Figure 2b), the general prediction criterion denoted $\widehat{MSE}$ over our random forest hence writes

$$ (9) \qquad \widehat{MSE}(t, \hat{M}) = \frac{1}{B} \sum_{b=1}^{B} \widehat{MSE}_{b}(t, \hat{\mu}_{b}) $$

As pointed out in Bouaziz (2024), two different models may lead to similar MSE values over time underlining the difficulty in assessing which model is better. A score is thus introduced to represent the prediction gain compared to a reference estimator and we define for each tree b:

$$ (10) \qquad Score_{b}(t, \hat{\mu}_{b}, \hat{\mu}_{b,0}) = \widehat{MSE}_{b}(t, \hat{\mu}_{b,0}) - \widehat{MSE}_{b}(t, \hat{\mu}_{b}) $$

Where $\hat{\mu}_b$ is the evaluated estimator and $\hat{\mu}_{b,0}$ the reference estimator over the $b$ samples. In our case, the reference estimator is the tree-specific non-parametric either the Nelson-Aalen or the Ghosh-Lin estimator described above. The ensemble score illustrated in Figure 2c) writes

$$ (11) \qquad Score(t, \hat{M}) = \frac{1}{B} \sum_{b=1}^{B} Score_{b}(t, \hat{\mu}_{b}, \hat{\mu}_{b,0}) $$

A higher score is associated with a better performance.

*Integrated counterparts.* Above MSE and derived score are time-dependent metrics. While they provide valuable insight of the performance for each time $t$, there is a need for the estimation of the expectation of single-time MSE and derived score over time (shaded areas in Figure 2). As demonstrated in Bouaziz (2024), above MSE reduces to the Brier score when individuals experience one event at most. In the spirit of the integrated version of the Brier score between two time points $\tau_1$ and $\tau_2$, we integrate the MSE and the score:

$$ (12) \qquad \begin{cases} \widehat{MSE}(\tau_1, \tau_2, \hat{M}) = \frac{1}{\tau_2 - \tau_1} \int_{\tau_1}^{\tau_2} \widehat{MSE}(t, \hat{M}) dt \\ IScore(\tau_1, \tau_2, \hat{M}) = \frac{1}{\tau_2 - \tau_1} \int_{\tau_1}^{\tau_2} Score(t, \hat{M}) dt \end{cases} $$

With $\tau_1 = 0$ and $\tau_2$ the maximum event time on the original sample.
construction to reduce computational overhead, making RSF more accessible for real-time clinical
decision-making and large-scale genomic studies.

3. Methods

The proposed Algorithm 1 of RecForest is an extension of the RSF introduced by Ishwaran (2008) (12).
The first step consists in drawing bootstrap samples to prevent overfitting and capture inherent variability
within the original dataset. Then, survival trees are constructed on each bootstrap sample. Unlike the
original RSF, our approach accommodates for subsequent events by integrating statistical considerations
tailored for recurrent events analysis. As a last step, the algorithm aggregates the results over the constructed
recursive survival trees to obtain a comprehensive estimate.

<caption>**Algorithm 1. Overview of RecForest algorithm**</caption>

(1) Draw *B* bootstrap samples from the learning data;
(2) Draw a survival tree *b* extended to recurrent events;
At each node,
* *mtry* predictors are randomly selected with *mtry* ∈ N, *mtry* ≤ *p*;
* A greedy algorithm for optimal threshold research is used to maximize the test statistic;
The tree grows until the stopping rule is met based on the minimal number of events *minsplit* and the
minimal number of individuals in terminal nodes *nodesize*;
Estimate *β*<sub>*b*</sub> is computed;
(3) Estimate *M̄* is computed over the *B* trees.

Next subsections describe in further details how survival trees grow for constructing the random forest.
Additionally, we provide adequate metrics for the evaluation. Finally, we expound on the computation of
variable importance.

**3.1. Growing trees with recurrent events**

3.1.1. Splitting rules

At each node $h \in \mathcal{H}$, the ongoing subsample is split into two daughter nodes denoted $h^{(+)}$ and $h^{(-)}$. The aim of the split is to make the daughter nodes as different as possible with regards to the outcome. The splitting rule requires that each of the *mtry* randomly drawn variable is dichotomized. For continuous variables, random split points, quartiles, and deciles are considered. Let $\mathbf{x}_h = A, B$ be the dichotomized vector of a variable inherited from $h$. For each split, we compare the marginal mean functions $\mu_A(t)$ and $\mu_B(t)$ between the two groups A and B. The null hypothesis is their equality.
### 2.3.2. OOB errors

OOB errors are used for tuning hyperparameters and evaluating predictive performances and are computed on OOB samples. They are also particularly useful in the absence of external validation data or when dealing with low-dimensional original samples, where allocating a portion for validation is hardly affordable. OOB predictions are calculated by average predictions from OOB trees, and the error rate is complementary to 1. In this work, we consider the IMSE to assess the OOB error:

$$ (13) \qquad \text{OOB error} = \widehat{\text{IMSE}}^{\text{OOB}}(t, \hat{M}^{\text{OOB}}) $$

In this way, models exhibiting lower OOB errors are consistently favored. Of note, computing OOB errors is not recommended when the number of trees is low as each one of them may underfit.

### 2.4. Variable importance

The importance of a variable ($VImp$) is evaluated by permutation, corresponding to the impact of random perturbations in the sample on the OOB error (Breiman (2001)). To quantify the $VImp$ of a covariate, a performance metric, as previously defined, is calculated following the permutation of values associated with this covariate. The $VImp$ is determined as the difference between the original and permuted performance metrics. For covariate $j$ and considering $K$ permutations, $VImp(j)$ writes

$$ (14) \qquad VImp(j) = \frac{1}{K} \sum_{k=1}^{K} (\hat{\theta} - \hat{\theta}_{k}^{j}) $$

With $\hat{\theta} = \{-\widehat{\text{IMSE}}, \hat{C}\}$ the original performance metric and $\hat{\theta}^j$ the permuted performance metric. High relative values of $VImp$ indicate a loss of performance and lower/null values are interpreted as no performance for such covariates.

## 3. Simulation study

We propose the following simulation settings to illustrate the use of RecForest, inspired by Ishwaran et al. (2014); Bouaziz (2024). Simulation scenarios will cover multiple cases with associated covariates, with or without a terminal event, low- and high-dimensional data, and with or without missing
In absence of a terminal event, we use the two-sample test akin to the log-rank test (38). The test statistic writes $U(t) = \int_0^t \frac{Y_A(u)Y_B(u)}{Y_A(u)+Y_B(u)} (d\hat{\mu}_A(u) - d\hat{\mu}_B(u))$, with $Y_A(t)$ and $Y_B(t)$ the number at risk in each group at time $t$.

In the presence of a terminal event, we follow the marginal model from Ghosh-Lin (GL) within the single variable $\mathbf{x}_h$ (39). Acknowledging there are no further recurrence after the terminal event, the marginal mean up to $t$ associated with $\mathbf{x}_h$ is defined as $\mu_{\mathbf{x}_h}(t) = \mathbb{E}[N^*(t)|\mathbf{x}_h] = \mu_0(t) \times \exp(\beta\mathbf{x}_h)$ with $\mu_0$ left unspecified and $\beta$ the regression coefficient. To accommodate longitudinal variables, the model follows the marginal approach of Ghosh-Lin and we adopt the marginal rate function $d\mu_{\mathbf{x}_h}(t) = d\mu_0(t) \times \exp(\beta\mathbf{x}_h(t))$, where $\mathbf{x}_h(t)$ may vary over time. The Wald test statistic is then extracted from $\mu_{\mathbf{x}_h}$ and $d\mu_{\mathbf{x}_h}$ to test the null hypothesis of $\beta = 0$. The variable selected for node $h$ is the one that maximizes the adequate test statistic to generate $h^{(+)}$ and $h^{(-)}$, based on the presence of a terminal event and/or longitudinal variables.

### 3.1.2. Terminal node estimator

Let $\mathcal{B}_b$ be a bootstrap sample drawn from original data on which a tree $b$ is grown, and let $\mathbf{x}$ be a $p$-dimensional covariate vector dropped down the tree. Individuals with similar covariate trajectories are grouped into the same terminal node, denoted $\mathcal{T}_b(\mathbf{x})$, which contains the subset of individuals assigned to the same terminal node as $\mathbf{x}$ in tree $b$. For each terminal node, we define a node-specific estimator $\hat{\mu}_{\mathcal{T}_b}(t|\mathbf{x})$ as the estimated number of recurrent events up to time $t$, based on the individuals in $\mathcal{T}_b(\mathbf{x})$. This estimator is given by:

$$ \hat{\mu}_{\mathcal{T}_b}(t|\mathbf{x}) = \int_0^t \frac{\sum_{i \in \mathcal{T}_b(\mathbf{x})} \frac{Y_i(u)}{\hat{S}_{\mathcal{T}_b(\mathbf{x})}(u)} dN_i(u)}{\sum_{i \in \mathcal{T}_b(\mathbf{x})} \frac{Y_i(u)}{\hat{S}_{\mathcal{T}_b(\mathbf{x})}(u)}} $$

Where $Y_i(t)$ is the at-risk indicator for subject $i \in \mathcal{T}_b(\mathbf{x})$, and $\hat{S}_{\mathcal{T}_b(\mathbf{x})}$ is the Kaplan-Meier estimator of the survival function, estimated using the subset $\mathcal{T}_b(\mathbf{x})$.

In the absence of a terminal event, we have $\hat{S}_{\mathcal{T}_b(\mathbf{x})}(t) = 1$ for all $t$, and the estimator simplifies to the Nelson-Aalen form:
<table><caption>TABLE 1
Summary of investigated scenarios</caption><thead><tr><th rowspan="2"></th><th colspan="2">Without a terminal event</th><th colspan="2">With a terminal event</th></tr><tr><th>Low dimensional<br>{n = 250, p = 20}</th><th>High dimensional<br>{n = 250, p = 300}</th><th>Low dimensional<br>{n = 250, p = 20}</th><th>High dimensional<br>{n = 250, p = 300}</th></tr></thead><tbody><tr><th>Complete</th><td>x</td><td>x</td><td>x</td><td>x</td></tr><tr><th>Missing</th><td>x</td><td>x</td><td>x</td><td>x</td></tr><tr><th>Random</th><td>x</td><td>x</td><td>x</td><td>x</td></tr></tbody></table>

data. 250 learning sets and one external validation set were generated for each scenario with $n = 250$ individuals and $p$ covariates. Table 1 below summarizes investigated scenarios. Next subsections further detail simulation parameters for each case. For each scenario, we grow 100 trees for RecForest. We set the following values: the minimal number of events is `minsplit` = 5, the minimal number of individuals is `nodesize` = 10, and the number of random predictors at each node is `mtry` = {1, $\sqrt{p}$, $\log(p)$}. We compared RecForest with non-parametric estimators, as well as a semi-parametric GL model where possible. Performances were measured on the external validation set using the C-index, MSE, the score and their integrated versions.

## 3.1. Simulation scheme

### 3.1.1. With and without a terminal event

For $i = 1, \dots, n$, $p_0$-dimensional covariate vector, $X_i = (X_{i,1}, \dots, X_{i,p_0})$, $X_{i,1:|\frac{p_0}{2}|} \mathcal{B}(0.5)$ Bernoulli variables and $X_{i,\lfloor \frac{p_0}{2} \rfloor+1:p_0} \mathcal{N}(2, 0.5)$ Gaussian variables are simulated. Recurrent events are generated from a non-homogenous Poisson process $\lambda(t|X_i) = \lambda_0(t) \exp(\beta^T X_i)$ with $\lambda_0(t) = \frac{\alpha}{\gamma}(\frac{t}{\gamma})^{\alpha-1}$ a Weibull baseline, $\alpha = 2$ the shape parameter and $\gamma = 0.39$ the scale parameter. The first ten associated coefficients are non-zero, with $\beta = (\beta_1, \dots, \beta_{p_0})^T$ and $\beta_1 = \log(5)$, $\beta_{2:4} = \log(1.3)$, $\beta_{5:p_0} = \log(0.7)$. The true expected number of events is $\mu^*(t|X_i) = \int_0^t \lambda(u|X_i)du = (\frac{t}{\gamma})^\alpha \exp(\beta^T X_i)$. The censoring process is then simulated based on a uniform distribution $U(0, 3)$. With a terminal event, the recurrent event process and covariates are simulated in the same way as above. The censoring process is simulated based on a uniform distribution $U(0, 8)$. The terminal event is simulated using a Cox model with shape parameter is 8 and scale parameter is 1.8. The same covariates as the recurrent event process are included with same coefficients $\beta$. We set $p_0 = 10$.
$$ \hat{\mu}_{\mathcal{T}_b}(t|\mathbf{x}) = \int_0^t \frac{\sum_{i \in \mathcal{T}_b(\mathbf{x})} Y_i(u) \, dN_i(u)}{\sum_{i \in \mathcal{T}_b(\mathbf{x})} Y_i(u)}. $$

This terminal-node specific estimator captures the recurrent event dynamics of individuals with similar covariate patterns, and forms the building block for constructing ensemble estimators.

### 3.1.3. Pruning trees

A pruning strategy is essential to help find a trade-off to prevent overfitting and improve generalization performance of trees, within a reasonable computational time. Aligned with Devaux (2023), we suggest two stopping rules for each terminal node: i) a minimal number of events called *minsplit*, and ii) a minimal number of individuals called *nodesize* (14). The validation of either stopping rule designates the current node *h* as terminal.

### 3.1.4. Handling missing data

To tackle eventual missing data, we include an adaptive-tree imputation which addresses missing data during the tree-growing stage by selectively drawing from available, non-missing, in-bag data (12,40). At each node $h_b$ from tree *b*, the method entails imputing random non-missing information specifically from the selected variables. The imputed data is then utilized for making splits within the node $h_b$. Imputed values are reset to missing as the tree progresses to subsequent nodes.

## 3.2. From trees to random forests

### 3.2.1. Ensemble estimates

Once all *B* trees are grown from the $\mathcal{B}_{b=1,...,B}$ independent bootstrap samples, we aggregate the predictions from each tree to form the final ensemble estimate. For a given covariate vector $\mathbf{x}$, let $\hat{\mu}_b$ denote the marginal mean function estimated from terminal nodes $\mathcal{T}_b(\mathbf{x})$. The ensemble estimate $\hat{\mathcal{M}}$ is the aggregation of all *B* tree-specific estimates:

$$ \hat{\mathcal{M}}(t|\mathbf{x}) = \frac{1}{B} \sum_{b=1}^{B} \hat{\mu}_b(t|\mathbf{x}). $$

This ensemble prediction leverages the variability captured by bootstrap sampling and the recursive partitioning structure of each tree to smooth out individual tree estimators and improve predictive stability.
### 3.1.2. Low- and high-dimensional scenarios

To define low- and high-dimensional scenarios, we introduce $q$ independent noise covariates randomly drawn from a standard normal distribution and add them to simulated datasets. We set $q = 10$ for low-dimensional scenarios, and $q = 290$ in high-dimensional scenarios. The total number of covariates for each scenario is $p = p_0 + q$.

*Complete.* When we analyze scenarios that involve all $p$ generated covariates, we refer to these as ‘complete’ datasets analyses.

*Missing.* To simulate real-world conditions where datasets may have missing values, we intentionally introduced missing data. Specifically, we randomly set 5% of the covariate $X_1$ to NA across all individuals. This was done in a completely random manner, ensuring that the missing data does not follow any pattern and is not dependent on any other variables or the values of $X_1$ itself. The missing data mechanism is completely at random. We refer to such scenarios as ‘Missing’.

*Random.* We created scenarios where the covariates are generated independently of the recurrent events to simulate a situation where no underlying factors influence the counting process. In such cases, $q$ independent noise covariates are randomly drawn from a standard normal distribution. We created scenarios where the covariates are generated independently of the recurrent events to simulate a situation where no underlying factors influence the counting process. We set $q = 20$ for low-dimensional scenarios, and $q = 300$ in high-dimensional scenarios. The total number of covariates for each scenario is $p = q$, ensuring that all covariates are unassociated with the event data and are purely random. We refer to these scenarios as ‘Random’.

## 3.2. Results

Performances were assessed in a framework of 250 training sets and one external validation set. The non-parametric estimator uses no covariates, regardless of the dimensionality by construction. For the GL model, no variable selection was performed, meaning all $p$ covariates were included in the model. This limits the analysis for the GL model to low-dimensional scenarios only.
Of note, in case of no terminal event, the ensemble estimator provides a nonparametric estimate of the cumulative expected number of recurrent events.

### 3.2.2. OOB ensemble estimates

By bootstrap sampling theory, each tree in the forest is trained on a subset of the data, leaving out roughly 37% of the original observations, called the out-of-bag (OOB) sample. These OOB samples are used to construct OOB ensemble estimates, which serve as internal validation predictions without the need for a separate test set.

Let $\mathcal{B}_b$ denote the bootstrap sample used to build tree $b$, and define $\mathcal{O}_i = \{b: i \notin \mathcal{B}_b\}$ as the set of trees for which individual $i$ is out-of-bag. The OOB ensemble estimator for the set of individuals with covariate vector $\mathbf{x}$ is given by:

$$ \hat{M}^{OOB}(t|\mathbf{x}) = \frac{1}{|\mathcal{O}_i|} \sum_{b \in \mathcal{O}_i} \hat{\mu}_b(t|\mathbf{x}), $$

where $\hat{\mu}_b(t|\mathbf{x})$ is the marginal mean function estimated from terminal nodes $\mathcal{T}_b(\mathbf{x})$, defined as in the prior section. The OOB ensemble prediction thus averages only the tree-specific predictions from trees where the subject was not included in training.

These OOB predictions are particularly useful for assessing model performance, tuning hyperparameters, or computing variable importance, especially in settings where external validation sets are unavailable or expensive to reserve.

## 3.3. Performance

Performance metrics below indicate the ability of the model to predict well from training data to unseen data. In our case, unseen data are either from the OOB sample or external validation data.

### 3.3.1. Assessing performance with relevant metrics

For the assessment of performance, we introduce an extended version of the C-index and employ the mean-square error (MSE), a derived score, and their integrated versions (Figure 2).
<table><caption>TABLE 2
Means and standard deviations of the C-index without a terminal event</caption><thead><tr><th>Scenario\Model</th><th>Np</th><th>GL</th><th>RecForest<br>mtry = 1</th><th>RecForest<br>mtry = &radic;p</th><th>RecForest<br>mtry = log(p)</th></tr></thead><tbody><tr><td></td><td></td><td colspan="4">Low dimensional {n = 250, p = 20}</td></tr><tr><td>Complete</td><td></td><td>0.55 (0.12)</td><td>0.68 (0.08)</td><td>0.71 (0.04)</td><td>0.70 (0.05)</td></tr><tr><td>Missing</td><td>0.55 (0.04)</td><td>0.52 (0.10)</td><td>0.65 (0.14)</td><td>0.69 (0.15)</td><td>0.67 (0.15)</td></tr><tr><td>Random</td><td></td><td>0.49 (0.05)</td><td>0.56 (0.15)</td><td>0.54 (0.15)</td><td>0.58 (0.18)</td></tr><tr><td></td><td></td><td colspan="4">High dimensional {n = 250, p = 300}</td></tr><tr><td>Complete</td><td></td><td>/</td><td>0.67 (0.21)</td><td>0.70 (0.11)</td><td>0.70 (0.17)</td></tr><tr><td>Missing</td><td>0.55 (0.04)</td><td>/</td><td>0.60 (0.29)</td><td>0.64 (0.18)</td><td>0.63 (0.24)</td></tr><tr><td>Random</td><td></td><td>/</td><td>0.51 (0.31)</td><td>0.55 (0.25)</td><td>0.56 (0.29)</td></tr></tbody></table>

Np = non-parametric estimator; GL = Gosh-Lin model with no variable selection.

RecForest was trained with fixed values for *minsplit* = 5 and *nodesize* = 10. 250 learning sets and one external validation set were generated for each scenario. Values closer to 1 indicate higher performance.

### 3.2.1. Without a terminal event

On average, 62% of the individuals experienced at least one recurrent event, 46% had at least two recurrent events, 26% had at least five recurrent events, and circa four recurrent events per individual. Table 2 below reports performances in terms of C-index values. Overall, performances based on C-index values are greater in scenarios with neither missing data, nor high-dimensionality, both in average and in variability. As expected, scenarios with random inputs lead to randomness with C-index values neighboring 0.50 for each model. The non-parametric estimator provides an average C-index of 0.55. The GL model seems to suffer from not being well-specified, with average C-index values ranging from 0.49 to 0.55 where assessable. RecForest consistently outperforms, irrespective of mtry with values ranging from 0.64 up to 0.71 (random scenarios are not deemed for comparing performance). Besides, it is not impacted by the introduction of massive noisy data, as C-index values remain similar across low- and -high-dimensional scenarios. Table 3 outlines performances in terms of integrated scores. As checked with C-indices, there is no expectations in the interpretation of the random scenarios, hence there are not displayed. The non-parametric estimator is the reference model in the computation of the score. Similar conclusions may be drawn from the different scenarios with the outperformance of RecForest. Yet, higher variability is observed, especially when introducing missing data.
<figure><img src="image_3.png" /><figcaption>Figure 2. Illustration of the performance metrics with true and predicted cumulative number of events over time</figcaption></figure>

### Concordance index

Kim (2018) adapted the C-index to recurrent events and considered the number of events over time across individuals (20). This metric hence suffers from the potential bias in case of substantial variability in the follow-up times (15). Individuals with longer follow-up times or a higher number of events might indeed disproportionately influence the C-index calculation. To address this issue, we suggest using occurrence rates by computing event rates per unit time.

The proposed C-index is defined as the proportion of all concordant pairs of individuals where predicted occurrence rates are correctly ordered with respect to observed occurrence rates (as shown in Figure 2(a)). As occurrence rates can be calculated for all individuals, including censored ones, the proposed C-index is not partial and considers all individuals in the computation. For each individual $i = 1, ..., n$, we define their observed rate as the total number of observed events divided by their observed time, $r_i = N_i(\gamma_i)/\gamma_i$ and $\hat{r}_i = \hat{M}(\gamma_i|\mathbf{x}_i)/\gamma_i$, respectively, where $\gamma_i = \min(D_i, C_i)$ is the observed time (i.e., minimum of terminal event or censoring), and $\hat{M}(\gamma_i|\mathbf{x}_i)$ is the predicted cumulative number of events. In this work, the C-index is calculated as the proportion of all comparable pairs $(i, j)$ where the ordering of $\hat{r}_j$ and $\hat{r}_i$ matches that of $r_i$ and $r_j$. The C-index then writes

$$ \hat{C}(\hat{M}) = \frac{\sum_{i=1}^{n} \sum_{j=1}^{n} I(r_i > r_j) \times I(\hat{r}_i > \hat{r}_j)}{\sum_{i=1}^{n} \sum_{j=1}^{n} I(r_i > r_j)}. $$

Like other C-indices, the value of the above C-index falls within the range of 0 to 1, where 1 indicates perfect concordance, and values close to 0.5 suggest randomness in the model.
<table><caption>TABLE 3
Means and standard deviations of the integrated score without a terminal event</caption><thead><tr><th>Scenario\Model</th><th>GL</th><th>RecForest<br>mtry = 1</th><th>RecForest<br>mtry = &radic;p</th><th>RecForest<br>mtry = log(p)</th></tr></thead><tbody><tr><td></td><td></td><td colspan="3">Low dimensional {n = 250, p = 20}</td></tr><tr><td>Complete</td><td>50.45 (41.00)</td><td>208.10 (102.42)</td><td>539.49 (451.96)</td><td>161.12 (87.65)</td></tr><tr><td>Missing</td><td>35.78 (17.91)</td><td>325.30 (189.63)</td><td>498.75 (415.28)</td><td>258.90 (112.14)</td></tr><tr><td></td><td></td><td colspan="3">High dimensional {n = 250, p = 300}</td></tr><tr><td>Complete</td><td>/</td><td>309.47 (134.57)</td><td>388.20 (226.95)</td><td>355.65 (117.78)</td></tr><tr><td>Missing</td><td>/</td><td>398.70 (229.57)</td><td>574.50 (318.75)</td><td>475.20 (213.80)</td></tr></tbody></table>

GL = Gosh-Lin model with no variable selection. RecForest was trained with fixed values for *minsplit* = 5 and *nodesize* = 10. 250 learning sets and one external validation set were generated for each scenario. Higher values indicate higher performance.

### 3.2.2. With a terminal event

On average, 44% of individuals experienced a terminal event during the observation period. Overall, performance discrepancies in terms of C-index values (Table 4) were observed when dealing with missing data or randomness, with performance being notably lower compared to complete data, as expected. The non-parametric estimator exhibited poor performance (C-index = 0.52 (0.03)). In both low and high-dimensional datasets, RecForest tends to perform better with higher mtry values, always reporting higher C-index values compared to non-parametric estimator and GL model. However, it is notable that in the high-dimensional scenarios, the performance drop is more significant. GL model performs better than non-parametric estimator only with complete data scenarios (C-index = 0.57 (0.07)). Integrated scores for evaluating approaches with a terminal event are displayed in Table 5. We observe similar results than without a terminal event. RecForest consistently yields integrated score values exceeding 300. The decrease in performance compared to the GL model is evident, with IScore values of 110.86 (75.14) and 112.81 (77.64) observed in complete and missing data scenarios, respectively.

In summary of the simulation study, our findings illustrate RecForest superior performance across all examined scenarios. Unlike the comparator GL model, RecForest effectively addresses both missing data and high-dimensionality. Furthermore, in random scenarios, RecForest outputs randomness, implying its reliability when the input lacks discernible patterns.
Mean-squared error and derived score

No MSE measure has been adapted to recurrent events framework until very lately. Bouaziz (2023) filled this gap and suggested a generalization of the Brier score from Graf (2000) (21,41). For our problematic, for each tree b and $\hat{\mu}_b$ the marginal mean function estimated, we define

$$ \overline{MSE}_b(t, \hat{\mu}_b) = \frac{1}{n} \sum_{i=1}^{n} \left( \int_{0}^{t} \frac{dN_i(u)}{\hat{G}_c(u|x)} - \hat{\mu}_b(t|x) \right)^2, $$

where $\hat{G}_c(u|x) = 1 - \hat{G}(u - |x)$ is an estimator of $G_c(u|x) = 1 - G(u - |x)$ the conditional cumulative distribution function of the censoring variable C given x. If there is no terminal event, the censoring distribution G is estimated using the empirical cumulative distribution function of the censoring times, since all censoring times are fully observed. In the presence of a terminal event, however, C becomes incompletely observed due to the competing risk of the terminal event D. In this case, G is estimated using the Kaplan-Meier estimator, treating terminal events as censoring events. Moreover, since assume C and x to be independent, the conditional survival function $G(t|x)$ simplifies to the marginal G(t), and we denote its estimate by $\hat{G}(t)$. As suggested in Figure 2(b), the general prediction criterion denoted $\overline{MSE}$ over our random forest hence writes

$$ \overline{MSE}(t, \hat{M}) = \frac{1}{B} \sum_{b=1}^{B} \overline{MSE}_b(t, \hat{\mu}_b). $$

However, as pointed out in Bouaziz (2023), two different models may lead to similar MSE values over time due to the inseparability term (noise inherent in the process), thus underlining the difficulty in assessing which model is better (21). Therefore, a score is introduced to represent the prediction gain compared to a reference estimator and we define for each tree b:

$$ Score_b(t, \hat{\mu}_b, \hat{\mu}_{b,0}) = \overline{MSE}_b(t, \hat{\mu}_{b,0}) - \overline{MSE}_b(t, \hat{\mu}_b), $$

where $\hat{\mu}_b$ is the evaluated estimator and $\hat{\mu}_{b,0}$ the reference estimator over the b samples. In our case, the reference estimator is the tree-specific non-parametric either the Nelson-Aalen or the Ghosh-Lin estimator described above. The ensemble score illustrated in Figure 2(c) writes
<table><caption>TABLE 4
Means and standard deviations of the C-index with a terminal event</caption><thead><tr><th>Scenario\Model</th><th>Np</th><th>GL</th><th>RecForest<br>mtry = 1</th><th>RecForest<br>mtry = &radic;p</th><th>RecForest<br>mtry = log(p)</th></tr></thead><tbody><tr><td></td><td></td><td colspan="4">Low dimensional {n = 250, p = 20}</td></tr><tr><td>Complete</td><td></td><td>0.57 (0.07)</td><td>0.79 (0.05)</td><td>0.80 (0.04)</td><td>0.82 (0.04)</td></tr><tr><td>Missing</td><td>0.52 (0.03)</td><td>0.51 (0.11)</td><td>0.73 (0.19)</td><td>0.71 (0.13)</td><td>0.75 (0.11)</td></tr><tr><td>Random</td><td></td><td>0.45 (0.10)</td><td>0.53 (0.16)</td><td>0.51 (0.11)</td><td>0.50 (0.19)</td></tr><tr><td></td><td></td><td colspan="4">High dimensional {n = 250, p = 300}</td></tr><tr><td>Complete</td><td></td><td>/</td><td>0.71 (0.19)</td><td>0.69 (0.13)</td><td>0.74 (0.10)</td></tr><tr><td>Missing</td><td>0.52 (0.03)</td><td>/</td><td>0.64 (0.20)</td><td>0.68 (0.13)</td><td>0.71 (0.11)</td></tr><tr><td>Random</td><td></td><td>/</td><td>0.49 (23.10)</td><td>0.48 (12.30)</td><td>0.50 (20.09)</td></tr></tbody></table>

Np = non-parametric estimator; GL = Gosh-Lin model with no variable selection. RecForest was trained with fixed values for *minsplit* = 5 and *nodesize* = 10. 250 learning sets and one external validation set were generated for each scenario. Values closer to 1 indicate higher performance.

<table><caption>TABLE 5
Means and standard deviations of the integrated score with a terminal event</caption><thead><tr><th>Scenario\Model</th><th>GL</th><th>RecForest<br>mtry = 1</th><th>RecForest<br>mtry = &radic;p</th><th>RecForest<br>mtry = log(p)</th></tr></thead><tbody><tr><td></td><td></td><td colspan="3">Low dimensional {n = 250, p = 20}</td></tr><tr><td>Complete</td><td>110.86 (75.14)</td><td>315.76 (119.41)</td><td>446.14 (410.88)</td><td>410.90 (115.54)</td></tr><tr><td>Missing</td><td>112.81 (77.64)</td><td>368.62 (211.11)</td><td>406.20 (275.57)</td><td>392.85 (138.23)</td></tr><tr><td></td><td></td><td colspan="3">High dimensional {n = 250, p = 300}</td></tr><tr><td>Complete</td><td>/</td><td>547.89 (229.37)</td><td>589.14 (472.33)</td><td>628.67 (122.41)</td></tr><tr><td>Missing</td><td>/</td><td>392.34 (39.16)</td><td>578.52 (336.71)</td><td>512.85 (441.29)</td></tr></tbody></table>

GL = Gosh-Lin model with no variable selection. RecForest was trained with fixed values for *minsplit* = 5 and *nodesize* = 10. 250 learning sets and one external validation set were generated for each scenario. Higher values indicate higher performance.
$$ Score(t, \hat{M}) = \frac{1}{B} \sum_{b=1}^{B} Score_{b}(t, \mu_{b}, \mu_{b,0}) $$

A higher score is associated with a better performance.

### Integrated counterparts

Above MSE and derived score are time-dependent metrics. While they provide valuable insight of the performance for each time t, there is a need for the estimation of the expectation of single-time MSE and derived score over time (shaded areas in Figure 2). As demonstrated in Bouaziz (2023), above MSE reduces to the Brier score when individuals experience one event at most (21). In the spirit of the integrated version of the Brier score between two time points $\tau_1$ and $\tau_2$, we integrate the MSE and the score:

$$ \left\{ \begin{array}{l} \text{IMSE}(\tau_1, \tau_2, \hat{M}) = \frac{1}{\tau_2 - \tau_1} \int_{\tau_1}^{\tau_2} \overline{\text{MSE}}(t, \hat{M}) dt, \\ \text{IScore}(\tau_1, \tau_2, \hat{M}) = \frac{1}{\tau_2 - \tau_1} \int_{\tau_1}^{\tau_2} Score(t, \hat{M}) dt, \end{array} \right. $$

with $\tau_1 = 0$ and $\tau_2$ the maximum event time on the original sample.

### 3.3.2. OOB errors

OOB errors are used for tuning hyperparameters and evaluating predictive performances and are computed on OOB samples. They are also particularly useful in the absence of external validation data or when dealing with low-dimensional original samples, where allocating a portion for validation is hardly affordable. OOB predictions are calculated by average predictions from OOB trees, and the error rate is complementary to 1. In this work, we consider the IMSE to assess the OOB error:

$$ \text{OOB error} = \overline{\text{IMSE}}^{\text{OOB}}(t, \hat{M}^{\text{OOB}}) $$

In this way, models exhibiting lower OOB errors are consistently favored. Of note, computing OOB errors is not recommended when the number of trees is low as each one of them may underfit.
<figure><img src="image_5.png" /><figcaption>FIG 3. Event plot for readmission data</figcaption></figure>

4. Illustrative example: the readmission data

Readmission dataset from the frailtypack R package from Rondeau, Marzroui and Gonzalez (2012) is widely used to demonstrate methodological prin-
ciples from recurrent events analysis in presence of a terminal event. The
data consist of multiple rehospitalizations after surgery in 403 patients diag-
nosed with colorectal cancer. Available factors are sex (MF), chemotherapy
treatment (YesNo), Dukes' tumoral stage (with levels A-B, C, and D), and
time-dependent comorbidity Charlson's index (with levels 0, 1-2, and ≥ 3).
In average, there were 1.13 (min. – max. = 0 – 22) hospital readmissions
per patients, with 199 patients with no admission and a total of 106 deaths
(Figure 3).

In absence of an external validation set, performances were assessed with
a 10-fold cross-validation procedure. We consider the following models: four
multivariate Ghosh-Lin models with arbitrary combinations of factors, and
**RecForest**. The reference model is the non-parametric estimator. Hyperpa-
rameters from **RecForest minsplit**, **nodesize** and **mtry** were tuned on the to-
tal sample and the OOB score was minimized for {ntrees = 100, minsplit =
2, nodesize = 1, mtry = 2} (Figure 6 in the Supplementary).

In our analysis (results in Table 6), the non-parametric estimator regis-
ters a C-index = 0.58 (0.05). **RecForest** outperforms with C-index = 0.80
(0.04). All GL models, with one to four covariates for adjustment, maintain
relatively consistent C-indices around 0.45 to 0.53. Comparing IMSE and
IScore metrics, **RecForest** and the non-parametric estimator are not directly
comparable due to construction. Specifically, the non-parametric reference
for the integrated score in **RecForest** is constructed for each bootstrap sam-
ple. Integrated scores for GL models operate on the overall dataset from the
### 3.4. Variable importance

The importance of a variable (VImp) is evaluated by permutation, corresponding to the impact of random perturbations in the sample on the OOB error (42). To quantify the VImp of a covariate, a performance metric, as previously defined, is calculated following the permutation of values associated with this covariate. The VImp is determined as the difference between the original and permuted performance metrics. For covariate j and considering K permutations, VImp(j) writes

$$VImp(j) = \frac{1}{K} \sum_{k=1}^{K} (\hat{\theta} - \hat{\theta}_{k}^{j})$$

With $\hat{\theta} = \{-\overline{IMSE}, \hat{C}\}$ the original performance metric and $\hat{\theta}^j$ the permuted performance metric. High relative values of VImp indicate a loss of performance and lower/null values are interpreted as no importance for such covariates.

## 4. Simulation study

We propose the following simulation settings to illustrate the use of RecForest, inspired by (13,21). Simulation scenarios will cover multiple cases with associated covariates, with or without a terminal event, low- and high-dimensional data, and with or without missing data. 250 learning sets and one external validation set were generated for each scenario with n = 250 individuals and p covariates. Next subsections further detail simulation parameters for each case.

For each scenario, we grow 100 trees for RecForest. To ensure comparability across scenarios and reduce computational burden over 250 replications, we fixed `nodesize = 10` and `minsplit = 5`, based on preliminary tuning that balanced tree depth with computational efficiency. We evaluated RecForest using three different values of `mtry` ∈ {1, $\sqrt{p}$, log(p)} to assess performance sensitivity to feature subsampling.

We compared RecForest with non-parametric estimators, as well as a semi-parametric GL model where possible. Performances were measured on the external validation set using the C-index, MSE, the score and their integrated versions.

### 4.1. Simulation scheme

#### 4.1.1. With and without a terminal event
<table><caption>TABLE 6
Means and standard deviations over the 10-fold cross-validation for readmission dataset</caption><thead><tr><th>Metric\Model</th><th>Np</th><th>GL1</th><th>GL2</th><th>GL3</th><th>GL4</th><th>RecForest</th><th>GL*</th></tr></thead><tbody><tr><td>C-index ↑</td><td>0.58<br>(0.05)</td><td>0.53<br>(0.08)</td><td>0.48<br>(0.08)</td><td>0.48<br>(0.07)</td><td>0.45<br>(0.05)</td><td>0.80<br>(0.04)</td><td>0.60<br>(0.06)</td></tr><tr><td>IMSE ↓</td><td>7 883.50<br>(6 229.47)</td><td>7 843.99<br>(6 106.36)</td><td>8 361.16<br>(6 292.29)</td><td>8 229.08<br>(6 478.35)</td><td>9 981.50<br>(6 064.23)</td><td>706.02<br>(508.96)</td><td>7 934.28<br>(6 606.23)</td></tr><tr><td>IScore ↑</td><td>ref.<br>ref.</td><td>39.41<br>(230.6)</td><td>-477.67<br>(348.48)</td><td>-345.62<br>(432.6)</td><td>-2 098.44<br>(541.59)</td><td>188.22<br>(89.00)</td><td>51.33<br>(142.63)</td></tr></tbody></table>

Np = non-parametric estimator; GL1 = Gosh-Lin model with sex; GL2 = Gosh-Lin with sex and chemotherapy; GL3 = Gosh-Lin model with sex, chemotherapy and Dukes' tumoral stage; GL4 = Gosh-Lin model with sex, chemotherapy and Dukes' tumoral stage and Charlson's index; GL* = Gosh-Lin model with best variables from RecForest.

Arrows indicate whether higher are lower scores lead to best performances.

ongoing fold. Consequently, IScore values from RecForest do not simply reflect the difference between IMSE values from the non-parametric estimator and RecForest, as opposed to IScores from GL models.

IMSE and IScore for RecForest indicate lower margin of errors. Among GL models, GL1 (Gosh-Lin model with sex) exhibits lower IMSE than the non-parametric estimator, resulting in a higher IScore. GL2 to GL4 (GL2 = Gosh-Lin with sex and chemotherapy; GL3 = Gosh-Lin model with sex, chemotherapy and Dukes' tumoral stage; GL4 = Gosh-Lin model with sex, chemotherapy and Dukes' tumoral stage and Charlson's index) yield negative IScore values, indicating high variability among GL models observed in our simulation study. Besides, high variability is observed across all approaches.

Variable importance for RecForest was based on both the C-index and the opposite of the integrated MSE (Figure 4). Most important variable identified by RecForest was the Charlson comorbidity index. Sex and chemotherapy did not seem to have an impact on the predictive performance. Variable selection enabled to reach better performance for GL* model.

Prediction curves for RecForest as the expected number of recurrent events are displayed in Figure 5. Predictions were generated for two patients, one with the highest Charlson comorbidity score (in orange), and the other with the lowest (in blue). We observe for the patient in orange that the model predicted an expected number of three readmissions as the patient dies after two observed readmissions. For the patient in blue, the model predictions are in line with observed events.
For $i = 1, ..., n, p_0$-dimensional covariate vector, $X_i = (X_{i,1}, ..., X_{i,p_0})$ is generated, where the first half of the components $X_{i,1:[\lfloor p_0/2 \rfloor]} \sim \mathcal{B}(0.5)$ are Bernoulli variables. The remaining covariates $X_{i,[\lfloor p_0/2 \rfloor]+1:p_0}$ are generated from a multivariate normal distribution with $X_{i,[\lfloor p_0/2 \rfloor]+1:p_0} \sim \mathcal{N}_{p_0-\lfloor p_0/2 \rfloor}(2, 0.5 \cdot \Sigma(\rho))$, where $\Sigma(\rho)$ is a covariance matrix with an autoregressive structure such that $\Sigma_{jk} = \rho^{|j-k|}$, and $\rho \in [0,1]$ controls the level of correlation among covariates. Such design allows to evaluate robustness under varying levels of multicollinearity among (continuous) covariates.

Recurrent events are simulated from a non-homogenous Poisson process with intensity $\lambda(t|X_i) = \lambda_0(t) \exp(\beta^T X_i)$, where the baseline hazard is Weibull with $\lambda_0(t) = \frac{\alpha}{\gamma} \left(\frac{t}{\gamma}\right)^{\alpha-1}$, $\alpha = 2$ (shape parameter), and $\gamma = 0.39$ (scale parameter). The vector of coefficients is $\beta = (\beta_1, ..., \beta_{p_0})^T$ and $\beta_1 = \log(5)$, $\beta_{2:4} = \log(1.3)$, $\beta_{5:p_0} = \log(0.7)$. The true expected number of events is given by $\mu^*(t|X_i) = \int_0^t \lambda(u|X_i) du = \left(\frac{t}{\gamma}\right)^\alpha \exp(\beta^T X_i)$. The censoring process is then simulated from a uniform distribution $U(0,3)$. With a terminal event, the recurrent event process and covariates are simulated in the same way as above.

The censoring process is simulated based on a uniform distribution $U(0,8)$, and the terminal event is simulated using a Cox model with the same covariates and coefficients $\beta$ with shape parameter 8 and scale parameter 1.8. We set $p_0 = 10$.

### 4.1.2. Low- and high-dimensional scenarios

To define low- and high-dimensional scenarios, we introduce $q$ independent noise covariates randomly drawn from a standard normal distribution and add them to simulated datasets. We set $q = 10$ for low-dimensional scenarios, and $q = 290$ in high-dimensional scenarios. The total number of covariates for each scenario is $p = p_0 + q$.

*Complete.* When we analyze scenarios that involve all $p$ generated covariates, we refer to these as 'complete' datasets analyses. Besides, we set $\rho \in \{0, 0.2, 0.5\}$.

*Missing.* To simulate real-world conditions where datasets may have missing values, we intentionally introduced missing data. Specifically, we randomly set 5% of the covariate $X_1$ to NA across all individuals. This was done in a completely random manner, ensuring that the missing data does not follow any pattern
<figure><img src="image_6.png" /><figcaption>FIG 4. Variable importance of RecForest computed on the C-index and the opposite of the integrated MSE. Charlson refers to Charlson comorbidity index, Dukes refers to tumoral Dukes stage.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG 5. Expected cumulative number of recurrent events with RecForest for two patients, one in orange with the highest Charlson comorbidity score, and the other in blue with the lowest. Data points outside the prediction curves are observed data. Triangle indicates the patient died.</figcaption></figure>
<table><caption>Table 1. Summary of investigated scenarios</caption><thead><tr><th rowspan="2"></th><th colspan="2">Without a terminal event</th><th colspan="2">With a terminal event</th></tr><tr><th>Low dimensional<br>{n = 250, p = 20}</th><th>High dimensional<br>{n = 250, p<br>= 300}</th><th>Low dimensional<br>{n = 250, p = 20}</th><th>High dimensional<br>{n = 250, p<br>= 300}</th></tr></thead><tbody><tr><td>Complete,<br>ρ ∈<br>{0, 0.2, 0.5}</td><td>x</td><td>x</td><td>x</td><td>x</td></tr><tr><td>Missing</td><td>x</td><td>x</td><td>x</td><td>x</td></tr><tr><td>Random</td><td>x</td><td>x</td><td>x</td><td>x</td></tr></tbody></table>

$ρ$ controls the level of correlation among covariates within the covariance matrix $Σ(ρ)$.

and is not dependent on any other variables or the values of $X_1$ itself. The missing data mechanism is completely at random. We set $ρ = 0$. We refer to such scenarios as ‘Missing’.

*Random*. We created scenarios where the covariates are generated independently of the recurrent events to simulate a situation where no underlying factors influence the counting process. In such cases, all $β$ were set to 0. We set $ρ = 0$. We refer to these scenarios as ‘Random’.

Table 1 summarizes investigated scenarios.

## 4.2. Results

For each simulation scenario, we generated 250 independent training datasets and one external validation set. Models were trained independently on each dataset (using fixed hyperparameter values for RecForest) and predictions were evaluated on the common external validation set. Performance metrics such as the C-index and integrated score were aggregated across the 250 replicates. The non-parametric estimator uses no covariates, regardless of the dimensionality by construction. For the GL model, no variable selection was performed, meaning all $p$ covariates were included in the model. This limits the analysis for the GL model to low-dimensional scenarios only.

*Without a terminal event*. On average, 62% of the individuals experienced at least one recurrent event, 46% had at least two recurrent events, 26% had at least five recurrent events, and circa four recurrent events per individual.
# 5. Discussion

We developed **RecForest** by extending the RSF algorithm to handle recurrent events in a survival framework, in potential presence of a terminal event, of longitudinal markers, and of missing data. To do so, the splitting rule at each node was tailormade for recurrent events analysis and mean cumulative number of events served in terminal node estimators. We characterized the performance both with discrimination and calibration by introducing a generalized C-index for recurrent event analysis and applying an innovative MSE.

In a simulation study, we compared **RecForest** with two baseline approaches: a non-parametric estimator and a GL model. Scenarios included variations of the presence of a terminal event, the low- or high-dimensionality of the data, and the inclusion of missing values. In instances where missing values or high-dimensional data were present, greater variability was observed across all scenarios for both metrics. In all explored cases, **RecForest** demonstrated higher performances both in terms of C-index and integrated score values compared with baseline. The impact was quite little when introducing massive noisy data. Besides, **RecForest** emerged as the sole modeling approach capable of handling high-dimensionality scenarios with no prior variable selection. Furthermore, we presented a practical application using well-known open-source data, showcasing how the fine-tuning of **RecForest** hyperparameters leads to a more performant model. Again, **RecForest** exhibited superior predictive performance, achieving a C-index of 0.80 alongside strong calibration metrics (IMSE = 706, IScore = 188). Overall, across both simulated and real-world datasets, **RecForest** consistently emerged as the most effective modeling approach.

In practical applications, high-dimensional problems involving recurrent events are often sidestepped by transforming the recurrent event survival framework into alternative formats, such as an event count, a time-to-first-event endpoint, or a classification problem. However, each of these transformations may lead to the voluntary omission of valuable information. In response to this, **RecForest** aims to bridge a recognized gap in handling such applications, ensuring a more comprehensive and nuanced analysis of recurrent event data. Additionally, our algorithm benefits from random forests features, i.e. the ability of handling missing data or multicollinearity, and reducing overfitting thanks to bagging principle.

Additional settings can be explored to integrate a terminal event within the proposed approach. For instance, Charles-Nelson, Katsahian and Schramm (2019) suggested working with inverse probability of survival weighting (IPSW)
Overall, performances based on C-index values are greater in scenarios with neither missing data, nor high-
dimensionality, both in average and in variability (Table 2). As expected, scenarios with random inputs
lead to randomness with C-index values neighboring 0.50 for each model. The non-parametric estimator
provides an average C-index of 0.55. The GL model seems to suffer from not being well-specified, with
average C-index values ranging from 0.49 to 0.55 where assessable. RecForest consistently outperforms
the Ghosh-Lin model across all scenarios. In the low-dimensional setting (p = 20), RecForest achieves a
higher C-index, particularly in the uncorrelated case (ρ = 0) where it reaches up to 0.71 (random scenarios
are not deemed for comparing performance). Although the C-index slightly decreases as correlation
increases (ρ = 0.2, 0.5), RecForest maintains values around 0.57-0.60, showing robustness to
multicollinearity. In the high-dimensional setting (p = 300), RecForest maintains strong discrimination
(C-index=0.68-0.70) across all levels of correlation. These results confirm that RecForest is effective in
ranking individuals by risk, even under complex, correlated, or high-dimensional covariate structures.
Besides, it is not impacted by the introduction of massive noisy data, as C-index values remain similar
across low- and -high-dimensional scenarios.

Table 3 outlines performances in terms of integrated scores. As checked with C-indices, there is no expectations in the interpretation of the random scenarios, hence there are not displayed. The non-parametric estimator is the reference model in the computation of the score. RecForest demonstrates clear advantages over the Ghosh-Lin model, particularly in the presence of multicollinearity. In the low-dimensional setting, Ghosh-Lin's scores deteriorate with increasing correlation, turning negative at ρ = 0.2 and ρ = 0.5, indicating poor calibration. In contrast, RecForest remains stable and well-calibrated, with scores improving slightly as correlation increases, especially for lower mtry values. In the high-dimensional setting, RecForest continues to deliver high scores even as pp increases, showing minimal degradation and outperforming its performance in the independent case (ρ = 0). These findings confirm RecForest's robustness to correlated predictors and its ability to produce accurate and well-calibrated predictions in both low- and high-dimensional scenarios. Notably, RecForest achieved better calibration under missing data than under complete data, as reflected in higher integrated score values. Our intuition is the regularizing effect of imputation may stabilize predictions in the presence of noisy or weakly informative covariates.
to compute coefficient weights in the Ghosh-Lin model, whereas we used inverse probability of censoring weighting. IPSW is typically recommended when modeling the terminal event is also of interest. Another example would be to use frailty models, either joint or additive as per Rondeau, Marzroui and Gonzalez (2012). Besides, natural extensions of random forests, serving as ensemble methods, have been widely used to improve performance through boosting techniques like Gradient Boosting (Friedman (2001)), Extreme Gradient Boosting (Chen and Guestrin (2016)), or LightGBM (Ke et al. (2017)). Since all these methods are grounded in tree-based structures, they offer seamless extensions to the proposed approach, and would hence provide innovative tools for recurrent event analysis.

Our methodology also suffers from several drawbacks. The primary limitation of random forest-like algorithms lies in the computation time, which grows with the number of trees, the dimensionality of the data and the numbers of variables selected at each tree node. Second, we assumed the proportional hazard assumption of the Gosh model, which may not universally hold in real-world settings. Furthermore, variable importance measures provided do not account for potential correlations. To address this limitation, the implementation of grouped variable importance statistics is a promising avenue for further refinement (Devaux et al. (2023); Gregorutti, Michel and Saint-Pierre (2015)). Nevertheless, signs of associations would still be unavailable. Another potential limitation of random forests is their static usage of features. Dynamic predictions could indeed be included as per Cottin et al. (2022) and Moradian et al. (2022).

On the other hand, the issue of interpretability in machine and deep learning, particularly in digital health has garnered significant attention, as pointed out in Farah et al. (2023). Several explainability methods have been proposed such as Local Interpretable Model-agnostic Explanations (LIME, Ribeiro, Singh and Guestrin (2016)), SHapley Additive exPlanations (SHAP, Lundberg and Lee (2017)), and counterfactual explanations (Guidotti (2022); Bhan et al. (2023)). These interpretability techniques have been recently adapted for survival analysis (Cottin et al. (2024); Kovalev, Utkin and Kasimov (2020)). Moreover, random forest-like methods offer a valuable tool for variable selection, especially in addressing high-dimensionality or obtaining hazard ratios that are intrinsically interpretable (Khan and Shaw (2016); Wang and Li (2017)). Approaches such as permutation-based selection, variable hunting, and iterative feature elimination serve as effective means towards this purpose (Genuer, Poggi and Tuleau-Malot (2010); Ishwaran et al. (2010); Pang et al. (2012)).
<table><caption>Table 2. Means and standard deviations of the C-index without a terminal event</caption><thead><tr><th>Scenario \ Method</th><th>Np</th><th>GL</th><th>RecForest<br>mtry = 1</th><th>RecForest<br>mtry = &radic;p</th><th>RecForest<br>mtry = log(p)</th></tr></thead><tbody><tr><td colspan="6" style="text-align:center;">Low dimensional {n = 250, p = 20}</td></tr><tr><td>Complete,<br>&rho; = 0</td><td>0.55 (0.04)</td><td>0.55 (0.12)</td><td>0.68 (0.08)</td><td>0.71 (0.04)</td><td>0.70 (0.05)</td></tr><tr><td>Complete,<br>&rho; = 0.2</td><td></td><td>0.54 (0.12)</td><td>0.59 (0.21)</td><td>0.60 (0.15)</td><td>0.59 (0.17)</td></tr><tr><td>Complete,<br>&rho; = 0.5</td><td></td><td>0.54 (0.15)</td><td>0.58 (0.24)</td><td>0.57 (0.10)</td><td>0.58 (0.16)</td></tr><tr><td>Missing</td><td></td><td>0.52 (0.10)</td><td>0.65 (0.14)</td><td>0.69 (0.15)</td><td>0.67 (0.15)</td></tr><tr><td>Random</td><td></td><td>0.45 (0.09)</td><td>0.52 (0.11)</td><td>0.53 (0.16)</td><td>0.53 (0.12)</td></tr><tr><td colspan="6" style="text-align:center;">High dimensional {n = 250, p = 300}</td></tr><tr><td>Complete,<br>&rho; = 0</td><td></td><td>/</td><td>0.67 (0.21)</td><td>0.70 (0.11)</td><td>0.70 (0.17)</td></tr><tr><td>Complete,<br>&rho; = 0.2</td><td></td><td>/</td><td>0.65 (0.25)</td><td>0.67 (0.13)</td><td>0.70 (0.17)</td></tr><tr><td>Complete,<br>&rho; = 0.5</td><td></td><td>/</td><td>0.64 (0.24)</td><td>0.67 (0.17)</td><td>0.68 (0.20)</td></tr><tr><td>Missing</td><td></td><td>/</td><td>0.60 (0.29)</td><td>0.64 (0.18)</td><td>0.63 (0.24)</td></tr><tr><td>Random</td><td></td><td>/</td><td>0.50 (0.25)</td><td>0.52 (0.29)</td><td>0.50 (0.27)</td></tr></tbody></table>

$\rho$ controls the level of correlation among covariates within the covariance matrix $\Sigma(\rho)$. Np = non-parametric estimator; GL = Ghosh-Lin model with no variable selection. RecForest was trained with fixed values for *minsplit* = 5 and *nodesize* = 10. 250 learning sets and one external validation set were generated for each scenario. Values closer to 1 indicate higher performance.

<table><caption>Table 3. Means and standard deviations of the integrated score without a terminal event</caption><thead><tr><th>Scenario \ Method</th><th>GL</th><th>RecForest<br>mtry = 1</th><th>RecForest<br>mtry = &radic;p</th><th>RecForest<br>mtry = log(p)</th></tr></thead><tbody><tr><td colspan="5" style="text-align:center;">Low dimensional {n = 250, p = 20}</td></tr><tr><td>Complete,<br>&rho; = 0</td><td>50.45 (41.00)</td><td>208.10 (102.42)</td><td>539.49 (451.96)</td><td>161.12 (87.65)</td></tr><tr><td>Complete,<br>&rho; = 0.2</td><td>-29.86 (25.50)</td><td>143.23 (90.05)</td><td>449.37 (400.00)</td><td>170.45 (92.04)</td></tr></tbody></table>
**6. Conclusion**

To conclude, we introduced a new algorithm based on survival theory for re-
current events with or without a terminal event and ensemble-based method-
ology for learning. **RecForest** is readily accessible to adequately answer fur-
ther clinical needs.

Funding

Author Murris J reports a grant from the Association Nationale de la Recherche et de la Technologie, with Pierre Fabre, Convention industrielle de formation par la recherche number 2020/1701.
<table><tbody><tr><td>Complete,<br/>ρ = 0.5</td><td>-54.50 (33.01)</td><td>123.83 (85.78)</td><td>402.10 (375.50)</td><td>139.92 (85.37)</td></tr><tr><td>Missing</td><td>35.78 (17.91)</td><td>325.30 (189.63)</td><td>498.75 (415.28)</td><td>258.90 (112.14)</td></tr><tr><td colspan="5">High dimensional {n = 250, p = 300}</td></tr><tr><td>Complete,<br/>ρ = 0</td><td>/</td><td>309.47 (134.57)</td><td>388.20 (226.95)</td><td>355.65 (117.78)</td></tr><tr><td>Complete,<br/>ρ = 0.2</td><td>/</td><td>256.81 (139.33)</td><td>327.55 (160.00)</td><td>301.55 (141.18)</td></tr><tr><td>Complete,<br/>ρ = 0.5</td><td>/</td><td>239.60 (125.66)</td><td>290.41 (149.59)</td><td>278.73 (134.48)</td></tr><tr><td>Missing</td><td>/</td><td>398.70 (229.57)</td><td>574.50 (318.75)</td><td>475.20 (213.80)</td></tr></tbody></table>

ρ controls the level of correlation among covariates within the covariance matrix Σ(ρ). GL = Ghosh-Lin model with no variable selection. RecForest was trained with fixed values for *minsplit* = 5 and *nodesize* = 10. 250 learning sets and one external validation set were generated for each scenario. Higher values indicate higher performance.

*With a terminal event.* On average, 44% of individuals experienced a terminal event during the observation period.

Overall, performance discrepancies in terms of C-index values (Table 4) were observed when dealing with missing data or randomness, with performance being notably lower compared to complete data, as expected. In the low-dimensional setting (p = 20), RecForest achieves C-index values as high as 0.82 under no correlation (ρ = 0) with *mtry* = *log(p)*, significantly outperforming GL at 0.52. While the C-index decreases slightly with increasing correlation, RecForest remains robust, achieving values above 0.66 even at ρ = 0.5, whereas GL fluctuates around 0.55. RecForest also outperforms GL under missing and random scenarios, maintaining C-indices in the range of 0.71–0.75 under missingness. In the high-dimensional setting (p = 300), RecForest continues to perform well, with C-index values ranging from 0.69 to 0.74 across all ρ levels and *mtry* values, demonstrating strong discrimination even when many features are present and correlated.

Integrated scores for evaluating approaches with a terminal event are displayed in Table 5. RecForest once again shows a clear advantage, particularly in high-dimensional and correlated settings. In the low-dimensional case with a terminal event, GL's score deteriorates significantly with increasing correlation (e.g., from 110.86 at ρ = 0 to -11.43 at ρ = 0.5), indicating poor calibration under multicollinearity.
# Supplementary Material

<figure><img src="image_8.png" /><figcaption>FIG 6. *Hyperparameter optimization on readmission data based on out-of-bag scores*</figcaption></figure>

## Hyperparameter optimization on readmission data based on out-of-bag scores

Hyperparameter optimization was performed using a grid search approach.

## References

* AMORIM, L. D. A. F. and CAI, J. (2015). Modelling recurrent events: a tutorial for analysis in epidemiology. *International Journal of Epidemiology* **44** 324–333.

* ANDERSEN, P. K. and GILL, R. D. (1982). Cox's Regression Model for Counting Processes: A Large Sample Study. *The Annals of Statistics* **10** 1100–1120.

* BHAN, M., VITTAUT, J.-N., CHESNEAU, N. and LESOT, M.-J. (2023). TIGTEC: Token Importance Guided TExt Counterfactuals. In *Machine Learning and Knowledge Discovery in Databases: Research Track: European Conference, ECML PKDD 2023, Turin, Italy, September 18–22, 2023, Proceedings, Part III* 496–512. Springer-Verlag, Berlin, Heidelberg.

* BOUAZIZ, O. (2024). Assessing model prediction performance for the expected cumulative number of recurrent events. *Lifetime Data Analysis* **30** 262–289.

* BREIMAN, L. (2001). Random Forests. *Machine Learning* **45** 5–32.

* CHARLES-NELSON, A., KATSAHIAN, S. and SCHRAMM, C. (2019). How to analyze and interpret recurrent events data in the presence of a terminal event: An application on readmission after colorectal cancer surgery. *Statistics in Medicine* sim.8168.
RecForest, on the other hand, maintains positive and high scores across all scenarios. Even with higher $\rho$, scores remain above 250 with $mtry = 1$, and even higher with larger $mtry$. In the high-dimensional case, RecForest performs strongly across the board, with scores increasing under lower correlation and remaining stable even when $\rho = 0.5$. While performance under missing data shows a slight drop, RecForest still maintains robust scores, demonstrating its ability to deliver well-calibrated scores even in challenging data conditions.

<table><caption>Table 4. Means and standard deviations of the C-index with a terminal event</caption><thead><tr><th>Scenario \ Method</th><th>Np</th><th>GL</th><th>RecForest<br>mtry = 1</th><th>RecForest<br>mtry = &radic;p</th><th>RecForest<br>mtry = log(p)</th></tr></thead><tbody><tr><td colspan="6" style="text-align: center;">Low dimensional {n = 250, p = 20}</td></tr><tr><td>Complete,<br>ρ = 0</td><td>0.52 (0.03)</td><td>0.57 (0.07)</td><td>0.79 (0.05)</td><td>0.80 (0.04)</td><td>0.82 (0.04)</td></tr><tr><td>Complete,<br>ρ = 0.2</td><td></td><td>0.57 (0.12)</td><td>0.69 (0.10)</td><td>0.76 (0.11)</td><td>0.76 (0.11)</td></tr><tr><td>Complete,<br>ρ = 0.5</td><td></td><td>0.55 (0.19)</td><td>0.61 (0.32)</td><td>0.67 (0.24)</td><td>0.66 (0.29)</td></tr><tr><td>Missing</td><td></td><td>0.51 (0.11)</td><td>0.73 (0.19)</td><td>0.71 (0.13)</td><td>0.75 (0.11)</td></tr><tr><td>Random</td><td></td><td>0.48 (0.10)</td><td>0.51 (0.19)</td><td>0.50 (0.17)</td><td>0.52 (0.15)</td></tr><tr><td colspan="6" style="text-align: center;">High dimensional {n = 250, p = 300}</td></tr><tr><td>Complete,<br>ρ = 0</td><td></td><td>/</td><td>0.71 (0.19)</td><td>0.69 (0.13)</td><td>0.74 (0.10)</td></tr><tr><td>Complete,<br>ρ = 0.2</td><td></td><td>/</td><td>0.71 (0.21)</td><td>0.70 (0.13)</td><td>0.74 (0.10)</td></tr><tr><td>Complete,<br>ρ = 0.5</td><td></td><td>/</td><td>0.69 (0.23)</td><td>0.70 (0.15)</td><td>0.73 (0.17)</td></tr><tr><td>Missing</td><td></td><td>/</td><td>0.64 (0.20)</td><td>0.68 (0.13)</td><td>0.71 (0.11)</td></tr><tr><td>Random</td><td></td><td>/</td><td>0.49 (0.22)</td><td>0.52 (0.17)</td><td>0.53 (0.20)</td></tr></tbody></table>

$\rho$ controls the level of correlation among covariates within the covariance matrix $\Sigma(\rho)$. Np = non-parametric estimator; GL = Ghosh-Lin model with no variable selection. RecForest was trained with fixed values for $minsplit$ = 5 and $nodesize$ = 10. 250 learning sets and one external validation set were generated for each scenario.
CHEN, T. and GUESTRIN, C. (2016). XGBoost: A Scalable Tree Boosting System. In *Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining*. KDD '16 785-794. Association for Computing Machinery, New York, NY, USA.

CHEN, S. and XU, C. (2023). Handling high-dimensional data with missing values by modern machine learning techniques. *Journal of Applied Statistics* **50** 786-804.

COOK, R. J. and LAWLESS, J. F. (1997). Marginal analysis of recurrent events and a terminating event. *Statistics in Medicine* **16** 911-924.

COOK, R. J., LAWLESS, J. F. and LEE, K.-A. (2010). A copula-based mixed Poisson model for bivariate recurrent events under event-dependent censoring. *Statistics in Medicine* **29** 694-707.

COTTIN, A., PECUCHET, N., ZULIAN, M., GUILLOUX, A. and KATSAHIAN, S. (2022). IDNetwork: A deep illness-death network based on multi-state event history process for disease prognostication. *Statistics in Medicine* **41** 1573-1598.

COTTIN, A., ZULIAN, M., PÉCUCHET, N., GUILLOUX, A. and KATSAHIAN, S. (2024). MS-CPFI: A model-agnostic Counterfactual Perturbation Feature Importance algorithm for interpreting black-box Multi-State models. *Artificial Intelligence in Medicine* **147** 1027-1041.

COX, D. R. (1972). Regression Models and Life-Tables. *Journal of the Royal Statistical Society: Series B (Methodological)* **34** 187-202.

DEVAUX, A., HELMER, C., GENUER, R. and PROUST-LIMA, C. (2023). Random survival forests with multivariate longitudinal endogenous covariates. *Statistical Methods in Medical Research* **32** 2331-2346.

FARAH, L., MURRIS, J. M., BORGET, I., GUILLOUX, A., MARTELLI, N. M. and KATSAHIAN, S. I. M. (2023). Assessment of Performance, Interpretability, and Explainability in Artificial Intelligence-Based Health Technologies: What Healthcare Stakeholders Need to Know. *Mayo Clinic Proceedings: Digital Health* **1** 120-138.

FRIEDMAN, J. H. (2001). Greedy function approximation: A gradient boosting machine. *The Annals of Statistics* **29** 1189-1232.

GENUER, R., POGGI, J.-M. and TULEAU-MALOT, C. (2010). Variable selection using random forests. *Pattern Recognition Letters* **31** 2225-2236.

GHOSH, D. and LIN, D. Y. (2000). Nonparametric Analysis of Recurrent Events and Death. *Biometrics* **56** 554-562.

GHOSH, D. and LIN, D. Y. (2002). Marginal Regression Models for Recurrent and Terminal Events. *Statistica Sinica* **12** 663-688.

GRAF, E., SCHMOOR, C., SAUERBREI, W. and SCHUMACHER, M. (1999). Assessment and comparison of prognostic classification schemes for sur-
<table><caption>Table 5. Means and standard deviations of the integrated score with a terminal event</caption><thead><tr><th>Scenario \<br>Method</th><th>GL</th><th>RecForest<br>mtry = 1</th><th>RecForest<br>mtry = &radic;p</th><th>RecForest<br>mtry = log(p)</th></tr></thead><tbody><tr><td colspan="5" style="text-align: center;">Low dimensional {n = 250, p = 20}</td></tr><tr><td>Complete,<br>p = 0</td><td>110.86 (75.14)</td><td>315.76 (119.41)</td><td>446.14 (410.88)</td><td>410.90 (115.54)</td></tr><tr><td>Complete,<br>p = 0.2</td><td>54.44 (12.80)</td><td>299.23 (140.90)</td><td>400.01 (310.15)</td><td>387.58 (306.06)</td></tr><tr><td>Complete,<br>p = 0.5</td><td>-11.43 (18.05)</td><td>250.39 (133.32)</td><td>385.95 (298.55)</td><td>362.05 (292.95)</td></tr><tr><td>Missing</td><td>112.81 (77.64)</td><td>368.62 (211.11)</td><td>406.20 (275.57)</td><td>392.85 (138.23)</td></tr><tr><td colspan="5" style="text-align: center;">High dimensional {n = 250, p = 300}</td></tr><tr><td>Complete,<br>p = 0</td><td>/</td><td>547.89 (229.37)</td><td>589.14 (472.33)</td><td>628.67 (122.41)</td></tr><tr><td>Complete,<br>p = 0.2</td><td>/</td><td>513.80 (213.98)</td><td>582.45 (436.72)</td><td>593.15 (228.43)</td></tr><tr><td>Complete,<br>p = 0.5</td><td>/</td><td>458.22 (401.37)</td><td>530.57 (387.55)</td><td>601.08 (334.19)</td></tr><tr><td>Missing</td><td>/</td><td>392.34 (39.16)</td><td>578.52 (336.71)</td><td>512.85 (441.29)</td></tr></tbody></table>

$\rho$ controls the level of correlation among covariates within the covariance matrix $\Sigma(\rho)$. GL = Ghosh-Lin model with no variable selection. RecForest was trained with fixed values for *minsplit* = 5 and *nodesize* = 10. 250 learning sets and one external validation set were generated for each scenario.

In summary of the simulation study, our findings illustrate RecForest superior performance across all examined scenarios. Unlike the comparator GL model, RecForest effectively addresses both missing data and high-dimensionality. Furthermore, in random scenarios, RecForest outputs randomness, implying its reliability when the input lacks discernible patterns.
vival data. *Statistics in Medicine* **18** 2529–2545.

GREGORUTTI, B., MICHEL, B. and SAINT-PIERRE, P. (2015). Grouped variable importance with random forests and application to multiple functional data analysis. *Computational Statistics & Data Analysis* **90** 15–35.

GUIDOTTI, R. (2022). Counterfactual explanations and how to find them: literature review and benchmarking. *Data Mining and Knowledge Discovery*.

HARRELL, F. E., LEE, K. L. and MARK, D. B. (1996). MULTIVARIABLE PROGNOSTIC MODELS: ISSUES IN DEVELOPING MODELS, EVALUATING ASSUMPTIONS AND ADEQUACY, AND MEASURING AND REDUCING ERRORS. *Statistics in Medicine* **15** 361–387.

HUANG, Y., LI, J., LI, M. and APARASU, R. R. (2023). Application of machine learning in predicting survival outcomes involving real-world data: a scoping review. *BMC medical research methodology* **23** 268.

ISHWARAN, H., KOGALUR, U. B., BLACKSTONE, E. H. and LAUER, M. S. (2008). Random survival forests. *The Annals of Applied Statistics* **2** 841–860.

ISHWARAN, H., KOGALUR, U. B., GORODESKI, E. Z., MINN, A. J. and LAUER, M. S. (2010). High-Dimensional Variable Selection for Survival Data. *Journal of the American Statistical Association* **105** 205–217.

ISHWARAN, H., GERDS, T. A., KOGALUR, U. B., MOORE, R. D., GANGE, S. J. and LAU, B. M. (2014). Random survival forests for competing risks. *Biostatistics (Oxford, England)* **15** 757–773.

KAPLAN, E. L. and MEIER, P. (1958). Nonparametric Estimation from Incomplete Observations. *Journal of the American Statistical Association* **53** 457–481.

KE, G., MENG, Q., FINLEY, T., WANG, T., CHEN, W., MA, W., YE, Q. and LIU, T.-Y. (2017). LightGBM: A Highly Efficient Gradient Boosting Decision Tree. In *Advances in Neural Information Processing Systems* **30**. Curran Associates, Inc.

KHAN, M. H. R. and SHAW, J. E. H. (2016). Variable selection for survival data with a class of adaptive elastic net techniques. *Statistics and Computing* **26** 725–741.

KIM, S., SCHAUBEL, D. E. and McCULLOUGH, K. P. (2018). A C-index for recurrent event data: Application to hospitalizations among dialysis patients. *Biometrics* **74** 734–743.

KOVALEV, M. S., UTKIN, L. V. and KASIMOV, E. M. (2020). SurvLIME: A method for explaining machine learning survival models. *Knowledge-Based Systems* **203** 106164.

LAWLESS, J. F. and NADEAU, C. (1995). Some Simple Robust Methods for
5. Illustrative example: the readmission data

Readmission dataset from the frailtypack R package is widely used to demonstrate methodological
principles from recurrent events analysis in presence of a terminal event (27,43). The data consist of
multiple rehospitalizations after surgery in 403 patients diagnosed with colorectal cancer. Available factors
are sex (M/F), chemotherapy treatment (Yes/No), Dukes' tumoral stage (with levels A-B, C, and D), and
Charlson comorbidity index (with levels 0, 1-2, and ≥ 3). In average, there were 1.13 (min. – max. = 0 –
22) hospital readmissions per patients, with 199 patients with no admission and a total of 106 deaths (Figure
3).

<figure><img src="image_4.png" /><figcaption>Figure 3. Event plot for readmission data</figcaption></figure>

In absence of an external validation set, performances were assessed with a 10-fold cross-validation procedure. We consider the following models: four multivariate Ghosh-Lin models with arbitrary combinations of factors, and RecForest. The Charlson comorbidity index was treated as a time-dependent covariate. The reference model is the non-parametric estimator.

Hyperparameter tuning was performed. explored multiple values of *mtry* ∈ {1, 2, 6} (also derived from *p* the number of covariates), along with *nodesize* ∈ {5, 10, 20} and *minsplit* ∈ {2, 5, 10}. The best-performing model, identified through hyperparameter tuning, was obtained with *mtry* = 4, *minsplit* = 2, and *nodesize* = 5. Further details are given in Appendix.

Table 6 summarizes model performance across a 10-fold cross-validation using the concordance index (C-index) and the integrated score (IScore). These metrics assess, respectively, the models' discriminative ability and calibration relative to a non-parametric reference.
the Analysis of Recurrent Events. *Technometrics* **37** 158–168.

LUNDBERG, S. M. and LEE, S.-I. (2017). A unified approach to interpreting model predictions. In *Proceedings of the 31st International Conference on Neural Information Processing Systems. NIPS'17* 4768–4777. Curran Associates Inc., Red Hook, NY, USA.

MORADIAN, H., YAO, W., LAROCQUE, D., SIMONOFF, J. S. and FRYD-MAN, H. (2022). Dynamic estimation with random forests for discrete-time survival data. *Canadian Journal of Statistics* **50** 533–548.

MURRIS, J., CHARLES-NELSON, A., TADMOURI SELLIER, A., LAVENU, A. and KATSAHIAN, S. (2023). Towards filling the gaps around recurrent events in high dimensional framework: a systematic literature review and application*. *Biostatistics & Epidemiology* **7** e2283650.

OZGA, A.-K., KIESER, M. and RAUCH, G. (2018). A systematic comparison of recurrent event models for application to composite endpoints. *BMC medical research methodology* **18** 2.

PANG, H., GEORGE, S. L., HUI, K. and TONG, T. (2012). Gene selection using iterative feature elimination random forests for survival outcomes. *IEEE/ACM transactions on computational biology and bioinformatics* / *IEEE, ACM* **9** 1422–1431.

PRENTICE, R. L., WILLIAMS, B. J. and PETERSON, A. V. (1981). On the regression analysis of multivariate failure time data. *Biometrika* **68** 373–379.

RIBEIRO, M. T., SINGH, S. and GUESTRIN, C. (2016). "Why Should I Trust You?": Explaining the Predictions of Any Classifier. In *Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. KDD '16* 1135–1144. Association for Computing Machinery, New York, NY, USA.

RONDEAU, V., MARZROUI, Y. and GONZALEZ, J. R. (2012). frailtypack: An R Package for the Analysis of Correlated Survival Data with Frailty Models Using Penalized Likelihood Estimation or Parametrical Estimation. *Journal of Statistical Software* **47** 1–28.

TIBSHIRANI, R. (1997). The lasso method for variable selection in the Cox model. *Statistics in Medicine* **16** 385–395.

UNO, H., CAI, T., PENCINA, M. J., D'AGOSTINO, R. B. and WEI, L. J. (2011). On the C-statistics for evaluating overall adequacy of risk prediction procedures with censored survival data. *Statistics in Medicine* **30** 1105–1117.

VAN BELLE, V., PELCKMANS, K., VAN HUFFEL, S. and SUYKENS, J. A. K. (2011). Support vector methods for survival analysis: a comparison between ranking and regression approaches. *Artificial Intelligence in*
RecForest achieved the highest overall performance, with a C-index of 0.78 (SD = 0.07) and a positive IScore of 252.10 (14.69), indicating strong discrimination and improved calibration relative to the benchmark. In contrast, the best-performing parametric model (GL*), reached a C-index of 0.60 (0.06) and an IScore of 51.33 (142.63), outperforming other GL variants but still falling short of RecForest. The non-parametric estimator used as a reference model showed a C-index of 0.58 (0.05), which notably exceeded that of most GL models. This suggests that the added complexity of GL2 – GL4 did not translate into improved predictive performance and may reflect an underfitting tendency or inability to capture non-linear patterns in the data.

Comparing IScore metrics, RecForest and the non-parametric estimator are not directly comparable due to construction. The IScore values reported represent relative calibration performance with respect to a reference estimator. This makes them suitable for comparison across models, since all methods are evaluated against the same benchmark. A higher IScore thus indicates better calibration with regards to the reference model (here the non-parametric estimator). Specifically, the non-parametric reference for the integrated score in RecForest is constructed for each bootstrap sample. Integrated scores for GL models operate on the overall dataset from the ongoing fold. Consequently, IScore values from RecForest do not simply reflect the difference between error values from the non-parametric estimator and RecForest, as opposed to IScores from GL models. Besides, IScore values for RecForest indicate lower margin of errors. Among GL models, GL1 (Ghosh-Lin model with sex) exhibits a higher IScore, compared to GL2 to GL4 (GL2 = Ghosh-Lin with sex and chemotherapy; GL3 = Ghosh-Lin model with sex, chemotherapy and Dukes' tumoral stage; GL4 = Ghosh-Lin model with sex, chemotherapy and Dukes' tumoral stage and Charlson's index) yield negative IScore values, indicating a worse calibration than a non-parametric estimator. RecForest reports not only positive IScore indicating better calibration than reference model, but also higher raw value. Besides, high variability is observed across all approaches.

Variable importance for RecForest was based on both the C-index and the opposite of the integrated MSE (Figure 4). Most important variable identified by RecForest was the Charlson comorbidity index. Sex and chemotherapy did not seem to have an impact on the predictive performance. Variable selection enabled to reach better performance for GL* model.

Prediction curves for RecForest as the expected number of recurrent events are displayed in Figure 5. Predictions were generated for two patients, one with the highest Charlson comorbidity score (in orange), and the other with the lowest (in blue). We observe for the patient in orange that the model overestimated
*Medicine* **53** 107–118.

WANG, H. and LI, G. (2017). A Selective Review on Random Survival Forests for High Dimensional Data. *Quantitative bio-science* **36** 85–96.

WEI, L. J., LIN, D. Y. and WEISSFELD, L. (1989). Regression Analysis of Multivariate Incomplete Failure Time Data by Modeling Marginal Distributions. *Journal of the American Statistical Association* **84** 1065–1073.
the expected number of three readmissions as the patient dies after two observed readmissions. For the patient in blue, the model predictions are in line with observed events.

In terms of computational time, RecForest remains tractable for practical use. Fitting a forest with 100 trees takes approximately one minute on a standard machine (M3 Pro, 18GB RAM). In comparison, the traditional Ghosh-Lin model fits in under 2 seconds. While RecForest involves greater computational cost due to tree construction and recurrent event-specific splitting, this remains manageable and is offset by the model's flexibility and ability to capture complex nonlinear effects.

<table><caption>**Table 6.** Means and standard deviations over the 10-fold cross-validation</caption><thead><tr><th>Metric</th><th>Np</th><th>GL1</th><th>GL2</th><th>GL3</th><th>GL4</th><th>RecForest</th><th>GL*</th></tr></thead><tbody><tr><td>C-index ↑</td><td>0.58<br>(0.05)</td><td>0.53<br>(0.08)</td><td>0.48<br>(0.08)</td><td>0.48<br>(0.07)</td><td>0.45<br>(0.05)</td><td>0.78<br>(0.07)</td><td>0.60<br>(0.06)</td></tr><tr><td>IScore ↑</td><td>ref.</td><td>39.41<br>(230.6)</td><td>-477.67<br>(348.48)</td><td>-345.62<br>(432.6)</td><td>-2 098.44<br>(541.59)</td><td>252.10<br>(14.69)</td><td>51.33<br>(142.63)</td></tr></tbody></table>

Np = non-parametric estimator; GL1 = Ghosh-Lin model with sex; GL2 = Ghosh-Lin with sex and chemotherapy; GL3 = Ghosh-Lin model with sex, chemotherapy and Dukes' tumoral stage; GL4 = Ghosh-Lin model with sex, chemotherapy and Dukes' tumoral stage and Charlson's index; GL* = Ghosh-Lin model with best variables from RecForest. Arrows indicate whether higher are lower scores lead to best performances.
<figure><img src="image_5.png" /><figcaption>Figure 4. Variable importance of RecForest computed on the C-index and the opposite of the integrated MSE. Charlson refers to Charlson comorbidity index, Dukes refers to tumoral Dukes stage.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Expected cumulative number of recurrent events with RecForest for two patients, one in orange with the highest Charlson comorbidity score, and the other in blue with the lowest. Data points outside the prediction curves are observed data. Triangle indicates the patient died.</figcaption></figure>
## 6. Discussion

We developed *RecForest* by extending the RSF algorithm to handle recurrent events in a survival framework, in potential presence of a terminal event, of longitudinal markers, and of missing data. To do so, the splitting rule at each node was tailormade for recurrent events analysis and mean cumulative number of events served in terminal node estimators. We characterized the performance both with discrimination and calibration by introducing a generalized C-index for recurrent event analysis and applying an innovative MSE.

In a simulation study, we compared *RecForest* with two baseline approaches: a non-parametric estimator and a GL model. Scenarios included variations of the presence of a terminal event, the low- or high-dimensionality of the data, multicollinearity amongst covariates, and the inclusion of missing values. In instances where missing values or high-dimensional data were present, greater variability was observed across all scenarios for both metrics. In all explored cases, *RecForest* demonstrated higher performances in terms of discrimination with the C-index with baseline GL models, and in terms of calibration with the IScore compared with the non-parametric estimator. The impact was quite little when introducing massive noisy data. Besides, *RecForest* emerged as the sole modeling approach capable of handling high-dimensionality scenarios with no prior variable selection. Furthermore, we presented a practical application using well-known open-source data, showcasing how the fine-tuning of *RecForest* hyperparameters leads to a more performant model. Again, *RecForest* exhibited superior predictive performance, achieving a C-index of 0.78 (0.07), alongside strong calibration (IScore = 252.10 (14.69)). Overall, across both simulated and real-world datasets, *RecForest* consistently emerged as the most effective modeling approach.

In practical applications, high-dimensional problems involving recurrent events are often sidestepped by transforming the recurrent event survival framework into alternative formats, such as an event count, a time-to-first-event endpoint, or a classification problem. However, each of these transformations may lead to the voluntary omission of valuable information. In response to this, *RecForest* aims to bridge a recognized gap in handling such applications, ensuring a more comprehensive and nuanced analysis of recurrent event data. Additionally, our algorithm benefits from random forests features, i.e. the ability of handling missing data or multicollinearity, and reducing overfitting thanks to bagging principle.

Additional settings can be explored to integrate a terminal event within the proposed approach. For instance, Charles-Nelson, et al. (2019) suggested working with inverse probability of survival weighting (IPSW) to compute coefficient weights in the Ghosh-Lin model, whereas we used inverse probability of censoring
weighting (7). IPSW is typically recommended when analyzing the terminal event is also of interest. Another example would be to use frailty models, either joint or additive as per Rondeau, et al. (2012) (27). Moreover, natural extensions of random forests, serving as ensemble methods, have been widely used to improve performance through boosting techniques like Gradient Boosting, Extreme Gradient Boosting, or LightGBM (44–46). Since all these methods are grounded in tree-based structures, they offer seamless extensions to the proposed approach, and would hence provide innovative tools for recurrent event analysis.

With regards to the MCAR mechanism adopted in our work, we acknowledge it provides a simplified benchmark to assess model robustness under incomplete data. While this assumption is strong, it allows for controlled interpretation of the missing data effects. Extending this framework to more realistic missing data mechanisms, such as Missing At Random (MAR), represents a promising direction for future work. Recent applied work in environmental science has demonstrated the practical relevance of random forest-based methods for handling complex missingness structures under MAR assumptions, including the use of multiple imputation strategies in air quality monitoring data (47).

Besides, in our work, each RecForest tree was trained on a bootstrapped sample drawn with replacement from the original dataset, with sample size approximately 2n/3, following the classical random forest framework (42). A valid alternative is subsampling without replacement, which is supported in many implementations, for instance in scikit-learn (48). Using subsampling can be advantageous in several scenarios, such as when working with very large datasets (where sufficient tree diversity can be achieved without replacement), when maintaining a fixed class balance is important, or when addressing imbalanced datasets where more control over class representation is desired. While we did not apply subsampling in this study, we acknowledge its potential to improve computational efficiency and consider it a promising direction for future optimization of the RecForest algorithm.

Our methodology also suffers from several drawbacks. The primary limitation of random forest-like algorithms lies in the computation time, which grows with the number of trees, the dimensionality of the data and the numbers of variables selected at each tree node. Typically, our use of 250 replicated simulations with a shared external validation set, along with reported standard deviations for key performance metrics, offers insights into the variability and robustness of RecForest across scenarios. Incorporating hyperparameter tuning within each replicate could further improve performance, albeit at a significantly increased computational cost.
Second, we assumed the proportional hazard assumption of the Ghosh model, which may not universally hold in real-world settings. In our analysis, the use of a log-rank-akin test for splitting in survival trees with no terminal event presents potential limitations. As highlighted by Cui et al. (2022), traditional log-rank-based splitting rules are prone to bias when censoring is dependent on covariates, especially under non-i.i.d. conditions (49). While in our work we make the strong assumption that the censoring process and the covariates are independent, this bias stems from the entanglement between failure and censoring distributions, leading to the selection of variables associated with censoring rather than the event process itself. The issue becomes even more pronounced in high-dimensional settings, where noise variables may carry marginal associations with the censoring distribution, thus misleading the splitting rule and compromising model consistency. Cui et al. (2023) further emphasize that such splitting bias can impair the estimation of treatment effects and highlight the importance of robust adjustments to account for censoring (50). To mitigate this, we recommend the use of *RecForest* with the Ghosh-Lin splitting even without a terminal event, which involves censoring weighting. This adjustment corrects for the splitting bias by more accurately estimating within-node failure distributions, even in the absence of terminal events, and enables convergence rates to depend only on the number of true failure-related variables rather than the full covariate space.

Furthermore, variable importance measures provided do not account for potential correlations. To address this limitation, the implementation of grouped variable importance statistics is a promising avenue for further refinement (14,51). Nevertheless, signs of associations would still be unavailable. Another potential limitation of random forests is their static usage of features. Dynamic predictions could indeed be included as per (52,53).

On the other hand, the issue of interpretability in machine and deep learning, particularly in digital health has garnered significant attention (54). Several explainability methods have been proposed such as Local Interpretable Model-agnostic Explanations (LIME), SHapley Additive exPlanations (SHAP), and counterfactual explanations (55–58). These interpretability techniques have been recently adapted for survival analysis (59–61). Moreover, random forest-like methods offer a valuable tool for variable selection, especially in addressing high-dimensionality or obtaining hazard ratios that are intrinsically interpretable (54,62,63). Approaches such as permutation-based selection, variable hunting, and iterative feature elimination serve as effective means towards this purpose (35,36,64).
## 7. Software availability

Although RSF and its derivatives have mostly been implemented in both R and Python, programs for the analysis of recurrent events in a survival framework mainly rely on packages in R. Consequently, we also developed our code in R to make it available to the wider community.

The proposed method is implemented in the open-source R package `recforest`, available on CRAN: https://cran.r-project.org/package=recforest. This package includes functions for training RecForest models, evaluating their performance both in terms of discrimination and calibration, and generating predictions on new data. Detailed documentation and examples are provided in the package vignettes.

## 8. Conclusion

To conclude, we introduced a new algorithm based on survival theory for recurrent events with or without a terminal event and ensemble-based methodology for learning. RecForest is readily accessible to adequately answer further clinical needs.

### List of abbreviations

C-index = Concordance index; GL = Ghosh-Lin; MSE = Mean Square Error; Np = Non-parametric; LASSO = Least Absolute Shrinkage and Selection Operator; OOB = Out-Of-Bag; RSF = Random Survival Forests; VImp = Variable Importance.

### Declarations

#### Ethics approval and consent to participate

Not applicable.

#### Consent for publication

Not applicable.

#### Availability of data and materials

Data are opensource under the frailtypack R package, fully available at: https://cran.r-project.org/web/packages/frailtypack/index.html.
## Conflicting Interests

The authors declare no conflict of interest.

## Funding

JM reports a grant from the Association Nationale de la Recherche et de la Technologie, with Pierre Fabre, Convention industrielle de formation par la recherche number 2020/1701.

## Appendix

To optimize the performance of RecForest, we conducted a grid search over key hyperparameters using the train_forest() function from the recforest R package. The goal was to identify the configuration that maximized the model's calibration, measured using the integrated score (mse_iscore). The hyperparameter search space was defined as follows:

* mtry ∈ {1, 2, 4}, (where p = 4, the number of covariates; corresponding to 1, $\sqrt{p}$, and p);

* minsplit ∈ {2, 5, 10};

* nodesize ∈ {5, 10, 20}.

For each hyperparameter combination, 100 trees were grown and the Ghosh-Lin splitting rule was used to take into account both a terminal event and the Charlson comorbidity index as a time-varying covariate.

A heatmap (Figure S1) was used to visualize the variation of the integrated score across the tuning grid. Each panel represents a different *mtry* value. Within each panel, the x-axis corresponds to *minsplit* and the y-axis to *nodesize*. Color gradients reflect the error, where lighter colors indicate better performance (higher score and lower error).
Figure S2 displays line plots of integrated scores mse_iscore by mtry, with separate panels for minsplit values and colored lines for each nodesize. This visualization further illustrates how model performance improves with increased tree complexity (*mtry* = 4), particularly under low minsplit and nodesize.

The grid search and visual analysis demonstrate that increased tree complexity and depth lead to improved model performance for recurrent event modeling in the readmission dataset. The final chosen model (*mtry* = 4, *minsplit* = 2, *nodesize* = 5) is used in subsequent analyses.

<figure><img src="image_7.png" /><figcaption>Figure S6. Heatmap of integrated scores across hyperparameter grid</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure S7. Line plot of integrated scores across hyperparameter grid</figcaption></figure>

References

1. Andersen PK, Gill RD. Cox's Regression Model for Counting Processes: A Large Sample Study. The Annals of Statistics. 1982;10(4):1100-20.

2. Prentice RL, Williams BJ, Peterson AV. On the regression analysis of multivariate failure time data. Biometrika. 1981;68(2):373-9.

3. Wei LJ, Lin DY, Weissfeld L. Regression Analysis of Multivariate Incomplete Failure Time Data by Modeling Marginal Distributions. Journal of the American Statistical Association. 1989;84(408):1065-73.

4. Cook RJ, Lawless JF, Lee KA. A copula-based mixed Poisson model for bivariate recurrent events under event-dependent censoring. Statistics in Medicine. 2010;29(6):694-707.

5. Amorim LDAF, Cai J. Modelling recurrent events: a tutorial for analysis in epidemiology. Int J Epidemiol. 2015;44(1):324-33.

6. Ozga AK, Kieser M, Rauch G. A systematic comparison of recurrent event models for application to composite endpoints. BMC Med Res Methodol. 2018;18(1):2.

7. Charles-Nelson A, Katsahian S, Schramm C. How to analyze and interpret recurrent events data in the presence of a terminal event: An application on readmission after colorectal cancer surgery. Statistics in Medicine. 2019;sim.8168.
8. Huang Y, Li J, Li M, Aparasu RR. Application of machine learning in predicting survival outcomes involving real-world data: a scoping review. BMC Med Res Methodol. 2023 Nov 13;23(1):268.

9. Cox DR. Regression Models and Life-Tables. Journal of the Royal Statistical Society: Series B (Methodological). 1972;34(2):187-202.

10. Tibshirani R. The lasso method for variable selection in the Cox model. Stat Med. 1997;16(4):385-95.

11. Van Belle V, Pelckmans K, Van Huffel S, Suykens JAK. Support vector methods for survival analysis: a comparison between ranking and regression approaches. Artif Intell Med. 2011 Oct;53(2):107-18.

12. Ishwaran H, Kogalur UB, Blackstone EH, Lauer MS. Random survival forests. The Annals of Applied Statistics. 2008 Sep;2(3):841-60.

13. Ishwaran H, Gerds TA, Kogalur UB, Moore RD, Gange SJ, Lau BM. Random survival forests for competing risks. Biostatistics. 2014 Oct;15(4):757-73.

14. Devaux A, Helmer C, Genuer R, Proust-Lima C. Random survival forests with multivariate longitudinal endogenous covariates. Stat Methods Med Res. 2023 Dec 1;32(12):2331-46.

15. Murris J, Charles-Nelson A, Tadmouri Sellier A, Lavenu A, Katsahian S. Towards filling the gaps around recurrent events in high dimensional framework: a systematic literature review and application*. Biostatistics & Epidemiology. 2023 Jan 2;7(1):e2283650.

16. Ghosh D, Lin DY. Nonparametric Analysis of Recurrent Events and Death. Biometrics. 2000;56(2):554-62.

17. Cook RJ, Lawless JF. MARGINAL ANALYSIS OF RECURRENT EVENTS AND A TERMINATING EVENT. Statist Med. 1997 Apr 30;16(8):911-24.

18. Harrell FE, Lee KL, Mark DB. Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Stat Med. 1996 Feb 28;15(4):361-87.

19. Uno H, Cai T, Pencina MJ, D'Agostino RB, Wei LJ. On the C-statistics for evaluating overall adequacy of risk prediction procedures with censored survival data. Statistics in Medicine. 2011;30(10):1105-17.

20. Kim S, Schaubel DE, McCullough KP. A C-index for recurrent event data: Application to hospitalizations among dialysis patients. Biometrics. 2018;74(2):734-43.

21. Bouaziz O. Assessing model prediction performance for the expected cumulative number of recurrent events. Lifetime Data Anal. 2024 Jan;30(1):262-89.

22. Cook RJ, Lawless JF. The statistical analysis of recurrent events. New York: Springer; 2007. (Statistics for biology and health).
23. Andersen PK. Statistical Models Based on Counting Processes. 1st ed. New York, NY: Springer New York; 1993. 1 p. (Springer Series in Statistics Ser).

24. Zhang Y, Wong G, Mann G, Muller S, Yang JYH. SurvBenchmark: comprehensive benchmarking study of survival analysis methods using both omics data and clinical data. GigaScience. 2022 Jul 30;11:giac071.

25. Andersen PK, Gill RD. Cox's Regression Model for Counting Processes: A Large Sample Study. Ann Statist [Internet]. 1982 Dec 1 [cited 2025 Apr 4];10(4). Available from: https://projecteuclid.org/journals/annals-of-statistics/volume-10/issue-4/Coxs-Regression-Model-for-Counting-Processes--A-Large-Sample/10.1214/ aos/1176345976 full

26. Prentice RL, Williams BJ, Peterson AV. On the regression analysis of multivariate failure time data. Biometrika. 1981;68(2):373-9.

27. Rondeau V, Marzroui Y, Gonzalez JR. frailtypack: An R Package for the Analysis of Correlated Survival Data with Frailty Models Using Penalized Likelihood Estimation or Parametrical Estimation. Journal of Statistical Software. 2012 Apr 17;47:1-28.

28. Wei LJ, Lin DY, Weissfeld L. Regression Analysis of Multivariate Incomplete Failure Time Data by Modeling Marginal Distributions. Journal of the American Statistical Association. 1989 Dec;84(408):1065-73.

29. Schmidli H, Roger JH, Akacha M. Estimands for Recurrent Event Endpoints in the Presence of a Terminal Event. Statistics in Biopharmaceutical Research. 2023 Apr 3;15(2):238-48.

30. Austin PC, Lee DS, Fine JP. Introduction to the Analysis of Survival Data in the Presence of Competing Risks. Circulation. 2016 Feb 9;133(6):601-9.

31. Andersen PK, Geskus RB, De Witte T, Putter H. Competing risks in epidemiology: possibilities and pitfalls. International Journal of Epidemiology. 2012 Jun;41(3):861-70.

32. Andersen PK, Keiding N. Multi-state models for event history analysis. Stat Methods Med Res. 2002 Apr 1;11(2):91-115.

33. Pickett KL, Suresh K, Campbell KR, Davis S, Juarez-Colunga E. Random survival forests for dynamic predictions of a time-to-event outcome using a longitudinal biomarker. BMC Med Res Methodol. 2021 Oct 17;21(1):216.

34. Chipman HA, George EI, McCulloch RE. BART: Bayesian additive regression trees. Ann Appl Stat [Internet]. 2010 Mar 1 [cited 2025 Apr 4];4(1). Available from: http://arxiv.org/abs/0806.3286

35. Ishwaran H, Kogalur UB, Gorodeski EZ, Minn AJ, Lauer MS. High-Dimensional Variable Selection for Survival Data. Journal of the American Statistical Association. 2010 Mar 1;105(489):205-17.

36. Genuer R, Poggi JM, Tuleau-Malot C. Variable selection using random forests. Pattern Recognition Letters. 2010 Oct 15;31(14):2225-36.
37. Pang H, George SL, Hui K, Tiejun Tong. Gene Selection Using Iterative Feature Elimination Random Forests for Survival Outcomes. IEEE/ACM Trans Comput Biol and Bioinf. 2012 Sep;9(5):1422-31.

38. Lawless JF, Nadeau C. Some Simple Robust Methods for the Analysis of Recurrent Events. Technometrics. 1995;37(2):158-68.

39. Ghosh D, Lin DY. Marginal Regression Models for Recurrent and Terminal Events. Statistica Sinica. 2002;12(3):663-88.

40. Chen S, Xu C. Handling high-dimensional data with missing values by modern machine learning techniques. Journal of Applied Statistics. 2023 Feb 17;50(3):786-804.

41. Graf E, Schmoor C, Sauerbrei W, Schumacher M. Assessment and comparison of prognostic classification schemes for survival data. Stat Med. 1999 Sep 15;18(17-18):2529-45.

42. Breiman L. Random Forests. Machine Learning. 2001 Oct 1;45(1):5-32.

43. Gonzalez JR, Fernandez E, Moreno V, Ribes J, Peris M, Navarro M, et al. Sex differences in hospital readmission among colorectal cancer patients. J Epidemiol Community Health. 2005 Jun;59(6):506-11.

44. Friedman JH. Greedy function approximation: A gradient boosting machine. The Annals of Statistics. 2001 Oct;29(5):1189-232.

45. Chen T, Guestrin C. XGBoost: A Scalable Tree Boosting System. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining [Internet]. New York, NY, USA: Association for Computing Machinery; 2016 [cited 2024 Feb 16]. p. 785-94. (KDD '16). Available from: https://dl.acm.org/doi/10.1145/2939672.2939785

46. Ke G, Meng Q, Finley T, Wang T, Chen W, Ma W, et al. LightGBM: A Highly Efficient Gradient Boosting Decision Tree. In: Advances in Neural Information Processing Systems [Internet]. Curran Associates, Inc.; 2017 [cited 2024 Feb 16]. Available from: https://papers.nips.cc/paper_files/paper/2017/hash/6449f44a102fde848669bdd9eb6b76fa-Abstract.html

47. Alsaber AR, Pan J, Al-Hurban A. Handling Complex Missing Data Using Random Forest Approach for an Air Quality Monitoring Dataset: A Case Study of Kuwait Environmental Data (2012 to 2018). International Journal of Environmental Research and Public Health. 2021 Jan;18(3):1333.

48. Pedregosa F, Pedregosa F, Varoquaux G, Varoquaux G, Org N, Gramfort A, et al. Scikit-learn: Machine Learning in Python. MACHINE LEARNING IN PYTHON.

49. Cui Y, Zhu R, Zhou M, Kosorok M. Consistency of Survival Tree and Forest Models: Splitting Bias and Correction. Statistica Sinica. 2022;32(3):1245-67.
50. Cui Y, Kosorok MR, Sverdrup E, Wager S, Zhu R. Estimating heterogeneous treatment effects with right-censored data via causal survival forests. Journal of the Royal Statistical Society Series B: Statistical Methodology. 2023 May 4;85(2):179-211.

51. Gregorutti B, Michel B, Saint-Pierre P. Grouped variable importance with random forests and application to multiple functional data analysis. Computational Statistics & Data Analysis. 2015 Oct 1;90:15-35.

52. Cottin A, Pécuchet N, Zulian M, Guilloux A, Katsahian S. IDNetwork: A deep illness-death network based on multi-state event history process for disease prognostication. Statistics in Medicine. 2022;41(9):1573-98.

53. Moradian H, Yao W, Larocque D, Simonoff JS, Frydman H. Dynamic estimation with random forests for discrete-time survival data. Canadian Journal of Statistics. 2022;50(2):533-48.

54. Farah L, Murris JM, Borget I, Guilloux A, Martelli NM, Katsahian SIM. Assessment of Performance, Interpretability, and Explainability in Artificial Intelligence-Based Health Technologies: What Healthcare Stakeholders Need to Know. Mayo Clinic Proceedings: Digital Health. 2023 Jun 1;1(2):120-38.

55. Ribeiro MT, Singh S, Guestrin C. "Why Should I Trust You?": Explaining the Predictions of Any Classifier. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining [Internet]. New York, NY, USA: Association for Computing Machinery; 2016 [cited 2022 Jun 30]. p. 1135-44. (KDD '16). Available from: https://doi.org/10.1145/2939672.2939778

56. Lundberg SM, Lee SI. A unified approach to interpreting model predictions. In: Proceedings of the 31st International Conference on Neural Information Processing Systems. Red Hook, NY, USA: Curran Associates Inc.; 2017. p. 4768-77. (NIPS'17).

57. Guidotti R. Counterfactual explanations and how to find them: literature review and benchmarking. Data Min Knowl Disc [Internet]. 2022 Apr 28 [cited 2024 Feb 16]; Available from: https://doi.org/10.1007/s10618-022-00831-6

58. Bhan M, Vittaut JN, Chesneau N, Lesot MJ. TIGTEC: Token Importance Guided TExt Counterfactuals. In: Machine Learning and Knowledge Discovery in Databases: Research Track: European Conference, ECML PKDD 2023, Turin, Italy, September 18-22, 2023, Proceedings, Part III [Internet]. Berlin, Heidelberg: Springer-Verlag; 2023 [cited 2024 Feb 16]. p. 496-512. Available from: https://doi.org/10.1007/978-3-031-43418-1_30

59. Spytek M, Krzyziński M, Langbein SH, Baniecki H, Wright MN, Biecek P. survex: an R package for explaining machine learning survival models. Bioinformatics. 2023 Dec 1;39(12):btad723.

60. Cottin A, Zulian M, Pécuchet N, Guilloux A, Katsahian S. MS-CPFI: A model-agnostic Counterfactual Perturbation Feature Importance algorithm for interpreting black-box Multi-State models. Artif Intell Med. 2024 Jan;147:102741.
61. Kovalev MS, Utkin LV, Kasimov EM. SurvLIME: A method for explaining machine learning survival models. Knowledge-Based Systems. 2020 Sep 5;203:106164.

62. Khan MHR, Shaw JEH. Variable selection for survival data with a class of adaptive elastic net techniques. Stat Comput. 2016 May 1;26(3):725-41.

63. Wang H, Li G. A Selective Review on Random Survival Forests for High Dimensional Data. Quant Biosci. 2017;36(2):85-96.

64. Pang H, George SL, Hui K, Tong T. Gene selection using iterative feature elimination random forests for survival outcomes. IEEE/ACM Trans Comput Biol Bioinform. 2012;9(5):1422-31.