📊 statistics

Principal stratification with recurrent events truncated by a terminal event: A nested Bayesian nonparametric approach

This paper proposes a nested Bayesian nonparametric framework using an enriched dependent Dirichlet process to address selection bias from death-truncated recurrent events, enabling transparent sensitivity analysis and demonstrating superior performance in estimating causal effects, such as the impact of exercise on rehospitalizations.

Original authors: Yuki Ohnishi, Michael O. Harhay, Guangyu Tong, Fan Li

Published 2026-03-18

📖 6 min read🧠 Deep dive

CC BY 4.0

Original authors: Yuki Ohnishi, Michael O. Harhay, Guangyu Tong, Fan Li

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are a doctor trying to figure out if a new exercise program helps heart failure patients. You want to know: Does exercise stop patients from being re-hospitalized?

But there's a catch. Some patients pass away during the study. When a patient dies, they stop coming to the hospital, so we stop counting their hospital visits. This creates a tricky problem:

If the exercise program helps people live longer, they have more time to get sick and get re-hospitalized.
If the exercise program doesn't work, people might die sooner, and we stop counting their hospital visits early.

If you just compare the total number of hospital visits between the exercise group and the control group, you might get a misleading answer. You aren't just measuring the effect of exercise on sickness; you are also measuring the effect of exercise on survival, which changes the rules of the game.

This paper proposes a clever new way to solve this puzzle using a method called Principal Stratification and a fancy statistical tool called Enriched Dependent Dirichlet Process (EDDP).

Here is the breakdown in simple terms:

1. The Problem: The "Survivor Bias" Trap

Imagine you are counting how many times people fall off their bikes.

Group A rides on smooth pavement.
Group B rides on rough, dangerous rocks.

If Group B crashes so hard that their bikes break and they stop riding entirely, you might think, "Wow, Group B had fewer falls!" But that's only because they stopped riding. They didn't have fewer falls; they just stopped having the opportunity to fall.

In the heart failure study, death is like the bike breaking. If we only look at the people still alive at the end, we are looking at a group that has been "filtered" by death. We need to compare apples to apples: What would have happened to the people who would have survived regardless of which group they were in?

2. The Solution: The "Always-Survivor" Club

The authors create a special imaginary club called the "Always-Survivor Stratum."

This club includes only the patients who would have lived past a certain date (say, 4 years) whether they did the exercise or not.
They exclude the "fragile" patients who would have died early no matter what.
By focusing only on this club, they can fairly ask: "For the people who were going to survive anyway, did exercise reduce their hospital visits?"

3. The Innovation: Two Clocks Instead of One

Previous methods tried to answer this question using just one clock. They asked: "How many hospital visits did the survivors have by Year 4?"

The Flaw: As time goes on, the "survivor" group changes. The people surviving at Year 1 are different from the people surviving at Year 4. It's like trying to measure the growth of a tree by looking at a different tree every year.

The Paper's New Idea: Two Clocks (Double-Indexing)
The authors use two clocks:

Clock R (The Survival Clock): Sets a fixed date (e.g., 4 years) to define who is in the "Always-Survivor" club. This group stays the same.
Clock T (The Event Clock): Counts the hospital visits up to a specific time (e.g., Year 1, Year 2, Year 3) within that fixed group.

The Analogy: Imagine you want to see how fast a car accelerates.

Old Way: You pick a car that is still running at the finish line and measure its speed at every mile. But the cars that broke down earlier are gone, so your data is messy.
New Way: You pick a specific group of cars that you know will finish the race (the "Always-Survivors"). Then, you watch that same group of cars at Mile 1, Mile 2, and Mile 3. You can now see exactly how their speed changes over time without the group changing underneath you.

4. The Secret Sauce: The "Enriched" Statistical Tool

To do this math, the authors had to invent a new statistical tool.

The Old Tool (Standard DP): Imagine you are sorting a pile of mixed nuts (almonds, walnuts, cashews) into buckets. The old tool was so obsessed with sorting the almonds (the frequent hospital visits) that it ignored the walnuts (the deaths). It ended up putting people with different death risks into the same bucket, which made the predictions about death very inaccurate.
The New Tool (EDDP): This is a "nested" sorter. It first sorts the nuts into big buckets based on the walnuts (death risk). Then, inside each big bucket, it sorts the almonds (hospital visits) into smaller sub-buckets.
- This ensures that people who are at high risk of dying are grouped together first, and then their hospital patterns are analyzed within that safe group. It's like sorting a library first by genre (Death Risk), and then by author (Hospital Pattern), rather than just by author.

5. The "What If" Game (Sensitivity Analysis)

There is one thing we can never know for sure: How are the "what if" scenarios related?

Question: If a patient would have survived 10 years with exercise, would they have also survived 10 years without it?
The authors admit they don't know the answer. So, they play a "What If" game. They run the math assuming different levels of connection between these two scenarios (from "totally unrelated" to "perfectly linked").
The Result: They found that for most reasonable guesses, the exercise program did reduce hospital visits for the "Always-Survivors." However, if you assume the two scenarios are perfectly linked (a very extreme guess), the benefit looks smaller. This transparency helps doctors understand how strong the evidence really is.

The Bottom Line

The authors applied this to the real HF-ACTION trial (a famous heart failure study).

What they found: The exercise program didn't just help people live longer; it actually reduced the total burden of hospital stays for the patients who were going to survive anyway.
The Takeaway: The benefit of exercise wasn't immediate; it grew over time. By using their new "Two Clock" method, they could see this slow, steady improvement that older methods missed because they were confused by the changing group of survivors.

In short: They built a better statistical microscope that lets us look at the "survivors" without getting confused by the people who died, revealing that exercise truly helps reduce the stress of repeated hospital stays for heart failure patients.

1. Problem Statement

The paper addresses a critical challenge in causal inference for clinical trials involving recurrent events (e.g., hospitalizations) that are truncated by a terminal event (e.g., death).

The Core Issue: In chronic diseases like heart failure, patients may experience multiple hospitalizations. However, death acts as an absorbing state that stops the observation of future hospitalizations. This creates selection bias: comparing observed hospitalization rates between treatment arms is confounded because the treatment may affect survival. Patients who survive longer have more "time at risk" for hospitalization, while those who die early have their recurrent processes truncated.
Limitations of Existing Methods:
- Single-Index Estimands: Previous approaches (e.g., Lyu et al., 2023) often define the "always-survivor" subpopulation at a single time point $t$ . This creates a "snapshot" effect where the target subpopulation changes as $t$ increases, making it difficult to distinguish between genuine treatment effects on disease burden and shifts in the composition of the survivor cohort.
- Parametric Assumptions: Existing methods often rely on strict parametric models (e.g., proportional hazards or specific frailty distributions), which are prone to misspecification bias in complex data-generating processes.
- Standard Dirichlet Process (DP) Pitfalls: Standard Bayesian nonparametric (BNP) approaches using DP priors for joint modeling often induce random partitions dominated by the frequent recurrent events. This leads to poor predictive performance for the terminal event (death) because the clustering is driven by the high-volume recurrent data rather than the sparse terminal events.

2. Methodology

The authors propose a Bayesian Nonparametric (BNP) framework within the Principal Stratification setting.

A. Causal Estimands: Double-Index Approach

Instead of a single time index, the authors define a Survivor-Average Number of Recurrence (SANR) using a double-index $(t, r)$ :

$r$ (Survivorship Horizon): Defines the principal stratum $AS(r) = \{i : D_i^0 > r, D_i^1 > r\}$ , the subpopulation of patients who would survive beyond time $r$ under both treatment arms. This fixes the target subpopulation.
$t$ (Event Horizon): Defines the time up to which recurrent events are counted ( $t < r$ ).
Estimand: $\text{SANR}(t; r) = g(\mu_1(t; r), \mu_0(t; r))$ , where $\mu_z(t; r)$ is the expected number of recurrent events up to time $t$ for the $AS(r)$ stratum under treatment $z$ .
Benefit: This decouples the definition of the subpopulation from the observation window, allowing researchers to track how treatment effects accumulate over time for a fixed latent subgroup.

B. Identification via Dual-Frailty and Sensitivity Analysis

Since potential outcomes under different treatments are never jointly observed, the cross-world dependence is non-identifiable.

Dual-Frailty Framework: The authors introduce a bivariate frailty vector $\gamma_i = (\gamma_i^0, \gamma_i^1)$ $γ_{i} = (γ_{i}^{0}, γ_{i}^{1})$ for each unit.
- $\gamma_i^z$ governs the within-world dependence between the terminal event and recurrent events under treatment $z$ .
- The correlation $\rho = \text{Cor}(\gamma_i^0, \gamma_i^1)$ captures the cross-world dependence.
Partial Identification: The distribution of frailties is estimable from data, but $\rho$ is not. The authors treat $\rho$ as a sensitivity parameter, performing transparent sensitivity analysis across a range of plausible values to assess the robustness of causal conclusions.

C. Bayesian Nonparametric Model: Enriched Dependent Dirichlet Process (EDDP)

To overcome the limitations of standard DP priors, the authors propose the Enriched Dependent Dirichlet Process (EDDP):

Nested Structure: The model uses a two-level partition:
1. Top-level ( $\phi$ -clusters): Clusters subjects based on the joint distribution of survival and recurrent events. This level is driven by the terminal event to ensure adequate learning about survival.
2. Nested level ( $\theta$ -subclusters): Within each $\phi$ -cluster, subjects are further subdivided based on finer details of the recurrent event process.
Covariate Dependence: The model incorporates covariates ( $X$ ) and time-varying covariates ( $V$ ) via a Dependent Dirichlet Process (DDP) structure, relaxing the exchangeability assumption.
Frailty Propagation: The subject-specific frailty $\gamma_i$ is sampled at the top level (governing survival) and propagated down to the recurrent event kernel. This ensures that the clustering for survival is not overwhelmed by the volume of recurrent events.
Inference: A fully tractable Gibbs sampling algorithm is developed, utilizing data augmentation to impute censored outcomes and latent cluster indicators.

3. Key Contributions

New Causal Estimand: Introduction of the double-index SANR estimand, which resolves the "snapshot" limitation of single-index principal stratification, allowing for the analysis of treatment effect accumulation over time within a fixed subpopulation.
Partial Identification Strategy: A dual-frailty framework that cleanly separates identifiable within-world dependence from non-identifiable cross-world dependence, enabling transparent sensitivity analysis via the parameter $\rho$ .
Novel BNP Prior (EDDP): Development of the Enriched Dependent Dirichlet Process, which specifically addresses the "recurrent-event dominance" problem in standard DP priors. It ensures that the terminal event (death) drives the primary clustering, leading to better predictive performance for survival.
Algorithmic Efficiency: A computationally efficient blocked Gibbs sampler for posterior inference in this complex hierarchical model.

4. Results

Simulation Studies

Setup: 300 datasets ( $n=1000$ ) generated from complex mixture models with latent heterogeneity.
Comparison: The proposed EDDPM was compared against Linear Models (LM), standard Dirichlet Process Mixtures (DPM), and Dependent Dirichlet Process Mixtures (DDPM).
Findings:
- Accuracy: EDDPM achieved the lowest Bias and Root Mean Square Error (RMSE) across all scenarios.
- Interval Estimation: EDDPM provided 95% credible intervals with coverage probabilities (CP) closest to the nominal 95% level and shorter average lengths (higher efficiency) compared to other methods.
- Robustness: Standard DP models suffered from poor survival prediction due to partitioning issues, while parametric LM models were sensitive to misspecification.

Application: HF-ACTION Trial

The method was applied to the HF-ACTION trial to evaluate the effect of structured aerobic exercise on rehospitalizations in chronic heart failure patients.

Model Fit: The EDDPM had the highest Log Pseudo Marginal Likelihood (LPML), indicating superior predictive fit compared to LM, DPM, and DDPM.
Causal Findings (Primary Analysis, $\rho=0.5$ ):
- Exercise training reduced the expected cumulative number of rehospitalizations among always-survivors.
- The benefit was time-dependent: The reduction in hospitalization burden became more pronounced over time (e.g., ~4% reduction at 1 year vs. ~10% at 4 years).
- Double-Index Insight: Fixing the survivorship stratum ( $r$ ) revealed that treatment effects varied across different risk strata. For instance, the evidence for benefit was stronger for patients who would survive to 4 years ( $r=1440$ ) compared to those surviving only to 3 years ( $r=1080$ ), a nuance obscured by single-index methods.
Sensitivity Analysis:
- Results were robust for moderate cross-world correlations ( $\rho \in [0.1, 0.7]$ ).
- Under an extreme assumption of near-perfect cross-world dependence ( $\rho=0.9$ , approximating a single shared frailty), the evidence for benefit was attenuated, and intervals included the null. This highlights the importance of the dual-frailty formulation in avoiding the masking of treatment effects.

5. Significance

Methodological Advancement: The paper bridges the gap between causal inference (principal stratification) and complex recurrent event analysis, providing a rigorous framework for handling truncation by death without relying on restrictive parametric assumptions.
Clinical Relevance: By applying the method to HF-ACTION, the authors demonstrate that structured exercise training offers a cumulative, long-term reduction in hospitalization burden for heart failure patients, a conclusion that might be missed or misinterpreted by conventional time-to-first-event or single-index analyses.
Practical Utility: The proposed framework offers a template for analyzing other chronic disease trials where morbidity (recurrent events) and mortality are intertwined, emphasizing the need to model cross-world dependence explicitly rather than assuming a single shared frailty.