Bayesian machine learning approach for recurrent events studies using Soft Bayesian Additive Regression Trees (SBART)

Abstract: Recurrent event data frequently arise in biomedical studies, where individuals may experience multiple recurrences of the same type of events, such as recurrent hospitalizations. This article introduces a nonparametric method for recurrent events under a Bayesian ensemble learning framework, called Soft Bayesian Additive Regression Trees (SBART), which combines multiple soft decision trees to achieve high predictive accuracy and a smooth estimator of the underlying intensity of the recurrent events. The proposed model represents the conditional intensity function of the non-homogeneous Poisson process as the product of a time-constant baseline, a subject-specific frailty random effect, and a nonparametric component capturing potentially nonlinear covariate effects and unknown interactions among covariates and time. A two-layer data augmentation scheme is employed to efficiently incorporate the SBART component within our computational algorithm. Simulation studies demonstrate that our method, called RecSBART in short, achieves superior accuracy in estimating cumulative intensity compared to existing approaches, even when our modeling assumptions are not true. With the Bayesian analysis of a study of recurrent hospitalizations of colorectal cancer patients, we further demonstrate our RecSBART method's ability to reveal and interpret the underlying complex relationships among covariates in a recurrent events study.

• The paper presents RecSBART, a fully Bayesian approach using soft decision trees for nonparametric modeling of subject-specific recurrent event intensities.
• It employs a two-level data augmentation strategy to efficiently address intractable integrals, achieving lower mean squared error than competing methods.
• Application to colorectal cancer hospitalization data demonstrates RecSBART's capability to capture nonlinear covariate-time interactions and deliver superior generalization.

Bayesian Machine Learning for Recurrent Events via Soft Bayesian Additive Regression Trees

Introduction

This paper presents a novel Bayesian machine learning methodology for the analysis of recurrent event data, leveraging Soft Bayesian Additive Regression Trees (SBART). The methodology, termed RecSBART, addresses major limitations in classical and machine learning-based recurrent event models by enabling nonparametric estimation of the subject-specific conditional event intensity under a non-homogeneous Poisson process (NHPP) framework, with explicit modeling of complex, nonlinear covariate-time interactions and unobserved heterogeneity through frailty.

Methodological Framework

RecSBART models the conditional intensity of recurrent events as

$$\lambda_i(t \mid W_i, \mathbf{x}_i) = \lambda_0 W_i \Phi(b(t, \mathbf{x}_i)),$$

where $\lambda_0$ is a baseline intensity parameter, $W_i$ is a subject-specific frailty (modeled as Gamma distributed), and $b(t, \mathbf{x}_i)$ is an unknown function of time and covariates, modeled nonparametrically via an ensemble of soft decision trees (SBART). The function $\Phi$ denotes the standard normal CDF, ensuring stability and upper-bounding the intensity.

The SBART extension over classical BART employs soft splits (using logistic weight functions) to induce smoothness and mitigate piecewise-constant artifacts inherent to BART. RecSBART supports flexible adaptation to unknown smoothness and automatically models high-order interactions without parametric commitment.

Integration of RecSBART into Bayesian inference for NHPP with frailty presents computational challenges due to the presence of intractable integrals and high-dimensional latent structure. These are efficiently addressed via a two-level data augmentation strategy: (1) latent event times for the thinned Poisson process, and (2) auxiliary variables for Gibbs updates of the SBART parameterization, leveraging the Poisson thinning framework and probit data augmentation, respectively.

Simulation Studies

The authors perform extensive simulation studies comparing RecSBART with RecForest (a frequentist random forest for recurrent events) and Bayesian proportional intensity frailty models. Simulated settings encompass correct and misspecified intensities, including both homogeneous and non-homogeneous processes and correct/incorrect frailty distributions.

RecSBART consistently yields the lowest average mean squared error (AMSE) in estimation of the cumulative intensity $\Lambda(t \mid \mathbf{x})$ across all tested scenarios, including when frailty is misspecified. For instance, in a challenging non-homogeneous scenario with frailty misspecification, RecSBART attains AMSE = 0.032 versus RecForest's 0.036, and the proportional intensity model's 1.035. These results underscore the robustness of RecSBART to both intensity and frailty model misspecification.

Assessment of subject-level frailty recovery via AMSE further demonstrates superiority of RecSBART over the parametric approach, with lower error in all tested scenarios.

Application: Recurrent Hospitalizations in Colorectal Cancer

RecSBART is used to analyze a well-studied biomedical recurrent event dataset of rehospitalizations following colorectal cancer surgery (403 patients, 458 events). Covariates include gender, chemotherapy status, Dukes' stage, and Charlson's index. Competing methods are again RecForest and the Bayesian proportional intensity frailty model.

Model fit is assessed via the empirical distribution of martingale residuals. RecSBART achieves a residual distribution that is most tightly centered around zero and exhibits the lightest tails, indicating superior fit over both RecForest and the proportional intensity model.

Figure 1: Empirical density of martingale residuals for RecSBART, RecForest, and the Bayesian proportional intensity model demonstrating RecSBART's lower bias and variance in residuals.

Estimation of the cumulative intensity function further reveals interpretable, nonlinear interactions among covariates and time, such as time-varying differential risk by gender and chemotherapy status not captured by standard models. Marginal effects and conditional covariate effect plots detect strong evidence against proportional hazards, as the effects are non-constant in time—a result the proportional intensity model cannot accommodate.

Model Generalizability and Overfitting Analysis

Out-of-sample predictive performance is assessed using 5-fold cross-validated mean squared martingale residuals (MSMRs). RecSBART maintains parity between training and test error (relative increase 8.7%), while RecForest displays a 28.4% test-train gap indicative of stronger overfitting and less stable generalization.

Figure 2: Five-fold cross-validated MSMRs for RecForest and RecSBART, illustrating superior generalization and reduced overfitting in RecSBART.

Practical and Theoretical Implications

RecSBART provides a flexible, computationally scalable solution to recurrent event modeling, supporting nonlinear and interactive covariate effects and explicit, interpretable frailty adjustment, all while delivering exact posterior Bayesian inference. Empirical robustness to model misspecification and improved generalization recommend RecSBART in biomedical and reliability studies where complex dependency structures are both plausible and consequential. The model's applicability is limited to vectorial covariates and outcomes but is readily extensible to ultrahigh-dimensional settings via sparsity-inducing priors and automated variable selection—a direction prompted by recent advances in high-dimensional SBART (2606.12701).

The demonstrated performance in recovering intricate covariate-time interactions and resilience under frailty and process misspecification positions RecSBART as a strong candidate for recurrent event analysis in settings where model fidelity is a critical design criterion.

Conclusion

RecSBART introduces a fully Bayesian, semiparametric approach for recurrent events analysis with robust and flexible modeling of subject-specific intensities and frailties. The model offers substantial improvements over classic proportional intensity models and tree-based machine learning competitors, in both simulation and real-world biomedical applications. Extension to structured (e.g., tensor) covariates and development of scalable, variable-selection regimes constitute promising avenues for future work.