RecSBART: Bayesian Machine Learning for Recurrent Events

• RecSBART is a Bayesian nonparametric method that models recurrent event data by combining subject-specific frailty with soft Bayesian additive regression trees.
• It decomposes event intensity into a time-constant baseline, a random frailty term, and a smooth regression component capturing nonlinear covariate and time interactions.
• Simulation and empirical studies demonstrate RecSBART’s competitive accuracy and reduced overfitting compared to traditional proportional intensity models and RecForest.

RecSBART is a Bayesian nonparametric method for recurrent event data in which each subject may experience multiple events over time, such as repeated hospitalizations. It models the subject-specific recurrence process under a non-homogeneous Poisson process (NHPP) framework whose conditional intensity is decomposed into a time-constant baseline intensity, a subject-specific frailty random effect, and a soft Bayesian additive regression tree (SBART) component that captures nonlinear covariate effects together with covariate-by-time interactions. The method was introduced as “Bayesian machine learning approach for recurrent events studies using Soft Bayesian Additive Regression Trees (SBART)” (Chen et al., 10 Jun 2026).

1. Problem setting and modeling objective

RecSBART is formulated for recurrent event studies in which, for subject $i=1,\dots,n$, $N_i(t)$ denotes the number of recurrent events in $(0,t]$, observed up to a terminal time $a_i$. The observed data are

$$\mathcal{D}=\{ 0<y_{i1}<\ldots<y_{in_i}< a_i;\mathbf{x}_i \in \mathbb{R}^p\}_{i=1}^n,$$

where $y_{ij}$ are recurrence times and $\mathbf{x}_i$ are fixed covariates. The motivating problem is that recurrent event studies often exhibit nonlinear, interactive, and time-varying covariate effects together with subject-level heterogeneity. Standard proportional intensity models are described as restrictive because they assume that the ratio of intensities between subjects does not depend on time, whereas marginal approaches typically do not model subject-specific trajectories or latent heterogeneity.

Within this setting, RecSBART is designed to combine frailty modeling with SBART’s flexible function estimation. This places it between classical semiparametric proportional intensity frailty models and purely marginal tree-based approaches such as RecForest. A central distinction is that RecSBART targets subject-specific recurrence processes rather than only marginal cumulative intensity (Chen et al., 10 Jun 2026).

2. Intensity decomposition and recurrent-event formulation

RecSBART assumes that, conditional on frailty $W_i$ and covariates $\mathbf{x}_i$, the recurrent event process follows an NHPP,

$$N_i(t) \sim \text{NHPP}(\lambda_i(t \mid W_i, \mathbf{x}_i)).$$

The conditional intensity is defined as

$N_i(t)$0

where $N_i(t)$1 is a baseline intensity, $N_i(t)$2 is a subject-specific frailty, $N_i(t)$3 is the standard normal CDF, and $N_i(t)$4 is an unknown nonparametric function. The factor $N_i(t)$5 guarantees positivity and boundedness of the regression component. For computation, the paper notes that $N_i(t)$6 can be reduced to a time-constant baseline $N_i(t)$7, yielding the working model

$N_i(t)$8

The three components have distinct roles. The baseline $N_i(t)$9 captures overall event frequency common to all subjects. The frailty $(0,t]$0 is a multiplicative random effect that accounts for unobserved heterogeneity and within-subject dependence. The SBART term $(0,t]$1 provides a smooth, flexible representation of nonlinear covariate effects, interactions among covariates, and covariate-time interactions. This structure departs from proportional intensity frailty models, in which covariate effects are typically log-linear and time-invariant (Chen et al., 10 Jun 2026).

3. Frailty prior and the SBART component

The frailties are assumed i.i.d. Gamma with mean 1,

$(0,t]$2

with density

$(0,t]$3

so that

$(0,t]$4

Fixing the mean at 1 serves identifiability. The model parameters are

$(0,t]$5

where $(0,t]$6 are the SBART tree parameters, and the prior factorizes as

$(0,t]$7

The paper specifies $(0,t]$8, $(0,t]$9, the SBART tree prior $a_i$0, and the leaf-parameter prior $a_i$1. The default SBART regularization settings are $a_i$2, $a_i$3, and $a_i$4, with split bandwidths satisfying

$a_i$5

shared across branches within a tree.

RecSBART inherits SBART’s modification of standard BART. Standard BART expresses a regression function as

$a_i$6

with each tree contributing a piecewise constant function,

$a_i$7

where $a_i$8 is an indicator selecting a terminal node. SBART replaces hard splits with soft probabilistic splits, so the tree contribution becomes a smooth weighted sum. The weights are products of soft decision probabilities along the path to each leaf:

$a_i$9

where $$\mathcal{D}=\{ 0<y_{i1}<\ldots<y_{in_i}< a_i;\mathbf{x}_i \in \mathbb{R}^p\}_{i=1}^n,$$0 is the set of ancestor nodes of leaf $$\mathcal{D}=\{ 0<y_{i1}<\ldots<y_{in_i}< a_i;\mathbf{x}_i \in \mathbb{R}^p\}_{i=1}^n,$$1, $$\mathcal{D}=\{ 0<y_{i1}<\ldots<y_{in_i}< a_i;\mathbf{x}_i \in \mathbb{R}^p\}_{i=1}^n,$$2 indicates whether the path goes right at node $$\mathcal{D}=\{ 0<y_{i1}<\ldots<y_{in_i}< a_i;\mathbf{x}_i \in \mathbb{R}^p\}_{i=1}^n,$$3, $$\mathcal{D}=\{ 0<y_{i1}<\ldots<y_{in_i}< a_i;\mathbf{x}_i \in \mathbb{R}^p\}_{i=1}^n,$$4 is the split point, and $$\mathcal{D}=\{ 0<y_{i1}<\ldots<y_{in_i}< a_i;\mathbf{x}_i \in \mathbb{R}^p\}_{i=1}^n,$$5 controls smoothness. The soft-splitting function is the inverse-logit CDF,

$$\mathcal{D}=\{ 0<y_{i1}<\ldots<y_{in_i}< a_i;\mathbf{x}_i \in \mathbb{R}^p\}_{i=1}^n,$$6

As $$\mathcal{D}=\{ 0<y_{i1}<\ldots<y_{in_i}< a_i;\mathbf{x}_i \in \mathbb{R}^p\}_{i=1}^n,$$7, SBART approaches ordinary BART. In RecSBART, this machinery is used to model the latent function $$\mathcal{D}=\{ 0<y_{i1}<\ldots<y_{in_i}< a_i;\mathbf{x}_i \in \mathbb{R}^p\}_{i=1}^n,$$8, thereby producing a smooth estimator for the intensity component $$\mathcal{D}=\{ 0<y_{i1}<\ldots<y_{in_i}< a_i;\mathbf{x}_i \in \mathbb{R}^p\}_{i=1}^n,$$9.

4. Likelihood, data augmentation, and posterior computation

Given the NHPP intensity, the likelihood is

$y_{ij}$0

The main computational difficulty is the integral

$y_{ij}$1

which must be repeatedly evaluated within MCMC. The paper’s principal computational device is a two-layer data augmentation scheme that converts the problem into a tractable binary-regression-like representation (Chen et al., 10 Jun 2026).

The first augmentation introduces latent thinned event times

$y_{ij}$2

generated from an NHPP with intensity

$y_{ij}$3

These latent points are independent of the observed event times. Conditional on them, the complete-data likelihood becomes

$y_{ij}$4

This step removes the difficult integral.

The second augmentation applies Albert–Chib latent Gaussian variables,

$y_{ij}$5

leading to the augmented likelihood

$y_{ij}$6

with truncation rules $y_{ij}$7 for thinned times and $y_{ij}$8 for observed event times. This makes the SBART update resemble standard Gaussian regression.

Posterior computation proceeds by MCMC. Several full conditionals are conjugate:

$y_{ij}$9

and

$\mathbf{x}_i$0

The frailty precision $\mathbf{x}_i$1 has a nonstandard conditional posterior and is updated via slice sampling. Given $\mathbf{x}_i$2, the SBART component is updated by standard tree-ensemble MCMC moves on tree structure and leaf parameters. The full algorithm iterates through latent thinned times, truncated Gaussian variables, $\mathbf{x}_i$3, $\mathbf{x}_i$4, $\mathbf{x}_i$5, and $\mathbf{x}_i$6, after which posterior inference can be drawn for cumulative intensity, frailties, covariate effects, time-varying effects, and uncertainty intervals.

5. Simulation study and comparative performance

The simulation study compares RecSBART with the semiparametric Bayesian proportional intensity frailty model of Sinha (1993) and with RecForest. The design uses $\mathbf{x}_i$7 subjects, $\mathbf{x}_i$8 covariates independently distributed as $\mathbf{x}_i$9, 20 replicated datasets per scenario, 50 trees for tree-ensemble methods, and 2500 burn-in plus 2500 sampling iterations for Bayesian methods. The covariate signal is

$W_i$0

Three scenarios are examined. Simulation A is a homogeneous Poisson process with model misspecification,

$W_i$1

Simulation B is a non-homogeneous Poisson process with model misspecification,

$W_i$2

Simulation C is an NHPP correctly specified relative to the RecSBART framework,

$W_i$3

Performance is assessed by the average mean squared error (AMSE) of the estimated cumulative intensity $W_i$4,

$W_i$5

with

$W_i$6

where $W_i$7. RecSBART attains the smallest AMSE in all three scenarios: $W_i$8 in Simulation A, $W_i$9 in Simulation B, and $\mathbf{x}_i$0 in Simulation C. The corresponding values are $\mathbf{x}_i$1, $\mathbf{x}_i$2, and $\mathbf{x}_i$3 for RecForest, and $\mathbf{x}_i$4, $\mathbf{x}_i$5, and $\mathbf{x}_i$6 for the proportional-intensity model (Chen et al., 10 Jun 2026).

The paper also evaluates Martingale residuals,

$\mathbf{x}_i$7

reporting that RecSBART and RecForest both produce residuals more concentrated around zero than the proportional intensity model, with RecSBART competitive or slightly better in some settings. For frailty estimation, the AMSE is $\mathbf{x}_i$8 for RecSBART versus $\mathbf{x}_i$9 for the proportional intensity model. Taken together, these results indicate accuracy under correct specification and robustness under misspecification. A plausible implication is that the combination of subject-specific frailty with flexible SBART regression stabilizes estimation even when the assumed frailty distribution or intensity form is imperfect.

6. Colorectal-cancer hospitalization study and Bayesian marginal effects

The empirical application analyzes recurrent hospitalizations among 403 colorectal cancer patients with 458 total recurrent hospitalizations. The covariates are gender ($$N_i(t) \sim \text{NHPP}(\lambda_i(t \mid W_i, \mathbf{x}_i)).$$0, $$N_i(t) \sim \text{NHPP}(\lambda_i(t \mid W_i, \mathbf{x}_i)).$$1 female, $$N_i(t) \sim \text{NHPP}(\lambda_i(t \mid W_i, \mathbf{x}_i)).$$2 male), chemotherapy status ($$N_i(t) \sim \text{NHPP}(\lambda_i(t \mid W_i, \mathbf{x}_i)).$$3, $$N_i(t) \sim \text{NHPP}(\lambda_i(t \mid W_i, \mathbf{x}_i)).$$4 placebo, $$N_i(t) \sim \text{NHPP}(\lambda_i(t \mid W_i, \mathbf{x}_i)).$$5 treatment), Dukes’ stage ($$N_i(t) \sim \text{NHPP}(\lambda_i(t \mid W_i, \mathbf{x}_i)).$$6, continuous), and Charlson’s index ($$N_i(t) \sim \text{NHPP}(\lambda_i(t \mid W_i, \mathbf{x}_i)).$$7, continuous). RecSBART and RecForest are each fit with 50 trees, and the Bayesian methods use 2500 burn-in plus 2500 retained iterations; for RecSBART, the prior $$N_i(t) \sim \text{NHPP}(\lambda_i(t \mid W_i, \mathbf{x}_i)).$$8 is used (Chen et al., 10 Jun 2026).

Goodness of fit is evaluated through Martingale residuals. The mean squared Martingale residuals on the full dataset are $$N_i(t) \sim \text{NHPP}(\lambda_i(t \mid W_i, \mathbf{x}_i)).$$9 for RecSBART, $N_i(t)$00 for RecForest, and $N_i(t)$01 for the proportional intensity model. In 5-fold subject-level cross-validation, the relative gap is $N_i(t)$02 for RecForest and $N_i(t)$03 for RecSBART. The paper interprets these results as indicating less overfitting for RecSBART.

The estimated cumulative intensity curves show that males have higher recurrence risk than females at the same Dukes’ stage, that higher Dukes’ stage is associated with higher risk, that placebo recipients tend to have higher rehospitalization risk than treated patients, and that these effects interact in nontrivial ways. This directly supports the methodological premise that the effect of one covariate may depend on the values of others.

To summarize fitted effects, the paper introduces Bayesian Marginal Effects (BME), adapted from local marginal analysis. For a subset $N_i(t)$04,

$N_i(t)$05

This quantity measures the average change in predicted cumulative intensity when selected covariates are changed while averaging over the observed distribution of the remaining covariates. The BME analysis suggests that the effect of gender increases over time, the effect of chemotherapy also changes over time, Dukes’ stage has a strong positive effect on cumulative intensity, and interactions are especially important among gender, chemotherapy, and Dukes’ stage. More specifically, at higher Dukes’ stage the gender-by-chemotherapy interaction is weaker, at lower Dukes’ stage the interaction is stronger, and among women the effect of Dukes’ stage varies more with chemotherapy status than among men.

The paper also examines BMEs on the log cumulative intensity as a diagnostic for proportionality. Because proportional intensity would imply time-constant BME on the log scale, the observed time variation is interpreted as evidence that a proportional intensity model is likely inappropriate.

7. Position relative to related approaches

RecSBART is distinguished along several axes. Relative to proportional intensity frailty models, it does not require linear covariate effects or proportionality over time. Relative to marginal recurrent-event methods, it provides subject-specific inference through frailties. Relative to RecForest, it models individual recurrence trajectories and latent heterogeneity rather than only marginal cumulative intensity. Relative to standard BART, it uses SBART’s soft splitting to obtain smoother estimation of continuous time-varying effects. Relative to survival-tree methods for first events, it is built directly for recurrent event processes rather than only first-event survival (Chen et al., 10 Jun 2026).

A common misconception would be to treat RecSBART as merely a recurrent-event version of a generic tree ensemble. The formulation in fact embeds the tree ensemble inside a frailty-based NHPP intensity model, so its target of inference is not only prediction of event counts but also the subject-specific stochastic intensity. Another misconception would be to view the method as enforcing proportionality through its multiplicative decomposition; the BME analysis on the log cumulative intensity is presented precisely to diagnose departures from proportional intensity.

In substantive terms, RecSBART is positioned for recurrent-event settings in which nonlinear covariate effects, interactions, time variation, and latent heterogeneity are all plausible. This suggests that its principal contribution lies not in replacing frailty models with tree ensembles, but in integrating the two within a coherent Bayesian framework for recurrent event analysis.