1 Introduction↩︎

In many clinical time-to-event studies, a non-negligible subset of patients never experience the event of interest because they are no longer susceptible to failure. For example, in breast cancer studies, patients may be successfully treated and enter long-term remission, effectively being “cured” of their cancer. In these settings, standard survival models can blur two distinct features of the outcome process: whether a patient remains at risk at all, and how quickly events occur among those who are still susceptible. One approach for addressing this setting is the use of a cure rate model, which allows the event time to have a positive probability of being infinite (see [1] and [2] for detailed treatments). This additional modeling is useful because it isolates the ways in which a treatment may provide benefit; for example, a treatment might shorten survival among individuals who are not cured while simultaneously increasing the probability of being cured.

The interpretation of cure models as distinguishing different ways in which a treatment may affect patient outcomes naturally suggests a causal perspective, yet there has been surprisingly little work in this direction. Comparing survival times only among observed failures across treatment groups is insufficient, because the population of individuals who experience an event under treatment may differ from the population who experience an event under control. Recently, [3] used principal stratification [4] to define causal effects within the subpopulation of patients who are “always uncured”, i.e., patients who would experience an event under either treatment arm. However, identification of these effects requires additional monotonicity assumptions on the treatment effect, and in many clinical trials, including the noninferiority CALGB 40101 [5] trial analyzed in this paper, such assumptions may be difficult to justify.

In this work, we define causal estimands tailored to the cure rate setting and develop flexible methods for estimating them. Our main contributions are summarized as follows:

1. We define useful parameters in the cure-rate setting to target using causal machine learning, including effects that decompose restricted mean survival time effects into effects attributable and not attributable to cure process.

2. We introduce BartCureas an effective method for obtaining estimates of population and conditional-average variants of these estimands. We then evaluate BartCureon both synthetic and real data, finding that it compares very favorably with other causal machine learning methods for estimating effects and detecting treatment effect heterogeneity.

3. We provide tools to use with BartCurefor interpreting the fitted model to understand the different sources of treatment effect heterogeneity.

Natural causal estimands in settings with a cured population are (i) the effect on the probability of an individual being cured, and (ii) the causal effect on the restricted mean survival time (RMST, [6]). To isolate the effect of the treatment on outcomes removing the effect of the treatment on cure, we also define a latency effect that decomposes the RMST effect into a part depending on the cure effect and a part depending on the latency effect; we also define a stochastic latency effect that we argue is more interpretable. We show that the latency effect is a well-defined contrast of distributions of potential outcomes on the same population, and therefore defines a causal effect; an advantage of the latency effect over the proposals of [3] is that it can be estimated with minimal assumptions beyond those typically made in causal survival analysis.

As a particular estimation strategy, we introduce BartCure, a promotion-time cure rate model [7], [8] based on Bayesian additive regression trees (BART, [9]), and provide evidence that it performs favorably relative to other causal machine learning techniques. BartCureis fully nonparametric in the sense that it can capture arbitrary covariate-dependent variation in the survival function over time. We provide a rigorous evaluation of its estimation accuracy and interval coverage for average and conditional average causal effects, as well as its ability to detect meaningful treatment effect heterogeneity. We believe that both the model itself and the novel data augmentation algorithm used to fit it are of independent methodological interest. An advantage of BartCure, and Bayesian machine learning more broadly, is the availability of direct uncertainty quantification through the posterior distribution. Empirically, these methods have demonstrated strong performance: BART has been a highly competitive prediction method since its introduction, and BART-based causal estimators have performed well in the ACIC data analysis competition [10], [11]. More recently, [12] showed that the strong empirical performance of BART-based methods extends to causal survival analysis.

We apply our methodology to a large clinical trial (CALGB 40101) that aimed to establish the noninferiority of the treatment paclitaxel (T, \(A_i = 1\)) for breast cancer relative to a baseline treatment cyclophosphamide+doxorubicin (CA, \(A_i = 0\)); because T has milder side effects than CA, this would allow T to be used as a new standard of care. To facilitate interpretation of causal effect estimates, we also show how to summarize the BartCureposterior. Our analysis shows that there is limited evidence of treatment effect heterogeneity in the relative performance of these treatments. We find weak evidence that the gap in survival is higher in older patients. We also see mild differences between CA and T over shorter time horizons, with the main treatment effect being that CA cures more individuals as opposed to lengthening the survival time of uncured individuals.

1.1 Related Work↩︎

This work contributes to a broader literature on causal machine learning in survival analysis. Beyond additional flexibility, we believe an important advantage of BartCureis that it reduces the “researcher degrees of freedom” [13] associated with analysts manually specifying nonlinearities and interaction effects. Several recent methods target treatment effect heterogeneity in time-to-event settings without explicitly modeling cure, including targeted learning approaches for conditional survival effects [14], Bayesian accelerated failure time models for individualized treatment effects [15], causal survival forests [16], and deep representation learning methods for treatment-specific hazards [17]. These methods provide flexible tools for handling censoring and treatment effect heterogeneity, but they generally treat survival as finite rather than explicitly modeling a cured subpopulation.

In the context of BART, several survival models for causal inference have been proposed. The best-known approach, introduced by [15], models the logarithm of the survival time as \(\log T_i(a) = r(a, X_i) + W_i\), where the error distribution \(W_i\) is modeled using a Dirichlet process mixture model [18]. This additive specification for \(T_i(a)\) induces an accelerated failure time (AFT) structure. The method is implemented in the AFTrees package in R. In extensive simulation studies, [12] showed that AFTrees compares very favorably with causal survival forests in both estimating heterogeneous causal effects on survival and quantifying the associated uncertainty. Outside the causal survival setting, BART has also been used for survival analysis with discrete-time hazards [19], semiparametric survival models [20], and competing-risk and recurrent-event models [21], [22]. Subsequent work has extended these ideas to clustered, interval-censored, and spatial survival settings [23], [24]; fully nonparametric survival models [25]–[27]; dynamic treatment regimes [28]; and relative survival models [29].

Several existing causal machine learning methods can accommodate a cured subpopulation. For example, [16] proposed causal survival forests (CSFs), which extend causal forests to survival settings and target finite-horizon survival functionals such as survival probabilities and restricted mean survival time. Under the assumption that there exists a time \(\tau\) beyond which individuals are considered cured, these methods remain applicable in the presence of a cured subpopulation. Another common approach for modeling survival data with a cured subpopulation is the mixture cure rate modeling; [30] adopt this framework, using a BART probit model for the probability of being cured and a Bayesian causal forest (BCF, [31]) model on log-survival times among uncured individuals. Unlike the promotion-time models considered here, mixture cure models parameterize the cure probability and the survival distribution among uncured individuals separately. Consequently, these approaches do not explicitly share information between the cure and latency components of the model.

2 Notation and Definition of Causal Parameters↩︎

We adopt standard notation from causal survival analysis. Let \(T_i\) denote the survival time, \(C_i\) a right-censoring time, \(X_i\) a vector of confounders and effect modifiers, and \(A_i\) a binary treatment indicator. Rather than observing \((T_i, C_i)\) directly, we observe only the event time \(Y_i = \min(T_i, C_i)\) and the event indicator \(\delta_i = 1(T_i \leq C_i)\). Unlike in standard survival settings, we allow individual \(i\) to be cured, in which case \(T_i = \infty\); such individuals are therefore necessarily censored. We use the potential outcomes framework to define causal effects [32]. For each treatment level \(a \in \{0,1\}\), let \(T_i(a)\) denote the potential survival time under treatment \(a\). Under the consistency assumption [33], the realized outcome is \(T_i = T_i(A_i)\). We also define the propensity score \(e(x) = \Pr(A_i = 1 \mid X_i = x)\).

2.0.0.1 Causal Effects

To define the causal effects of interest, it is useful to define some additional potential outcomes. First, we define the potential outcomes for the survival status indicator $ D_i(a, t) = 1{T_i(a) > t},$ for \(t \in [0, \infty]\). We then define the average causal effect on survival status at time \(t\) [34] and the causal effect on the cure event as \[\begin{align} \Delta_S(t) = \mathbb{E}\{D_i(1, t) - D_i(0, t)\} \quad \text{and} \quad \Delta_C = \mathbb{E}\{D_i(1, \infty) - D_i(0,\infty)\}, \end{align}\] respectively. Note that \(\mathbb{E}\{D_i(a, t)\} = \Pr(T_i(a) > t) = S_a(t)\), so \(\Delta_S(t)= S_1(t) - S_0(t)\) is the average causal effect on the probability of surviving beyond time \(t\). Note also that \(D_i(a, \infty)\) is never observed, and \(\Delta_C\) therefore is only recoverable under assumptions about the tail behavior of \(S_a(t)\). As a measure of lifetime gained up to time \(t\) due to treatment, we define the causal effect on restricted mean survival time (RMST) [34] as \[\begin{align} \Delta_R(t) = \mathbb{E}\{R_i(1, t) - R_i(0, t)\} \quad \text{where} \quad R_i(a,t) = \min\{T_i(a), t\}. \end{align}\] RMST can be interpreted as the life-expectancy-up-to-time-\(t\) [6] and so \(\Delta_R(t)\) represents the average amount of lifetime gained due to treatment up to time \(t\). Alternatively, one can write the causal RMST effect as the integrated difference in survival functions \(\Delta_R(t) = \int_0^t [\Pr\{T_i(1) > u\} - \Pr\{T_i(0) > u\}] \;du = \int_0^t \{S_1(u) - S_0(u)\} \;du\).

Conditional average versions of each of these effects can also be defined straightforwardly. For example, the conditional average treatment effect (CATE) version of \(\Delta_R(t)\) is \(\Delta_R(t,x) = \mathbb{E}\{R_i(1, t) - R_i(0, t) \mid X_i = x\}\). Similarly, we can define \(\Delta_C(x)\), \(\Delta_S(t, x)\), and so forth.

2.1 Identification Assumptions↩︎

Throughout, we make the standard positivity, ignorability, stable unit treatment value, consistency, and ignorable censoring assumptions that are standard for identifying causal effects in causal survival analysis. For completeness, we list these in the Supplementary Material as Assumptions 1—4. We also make the following assumption to identify the cured subpopulation.

2.1.0.1 Assumption 5: Cure Threshold.

There exists a known time \(\tau < \infty\) such that the hazard function satisfies \(h(t \mid a, x) = 0\) for all \(t > \tau\), \(a \in \{0,1\}\), and \(x\) in the support of \(X_i\). Equivalently, survival beyond \(\tau\) implies cure. Additionally, we require sufficient follow-up in the sense that \(\Pr(C_i > \tau \mid X_i = x) > 0\) for all \(x\).

Some tail behavior assumption for the survival function is required in order to conclude that an apparent asymptote in the survival curve truly reflects the presence of a cured subpopulation, since the event \(\{T_i = \infty\}\) is never directly observed. Assumption 5, which links \(\{T_i = \infty\}\) to the observable event \(\{T_i > \tau\}\), is a pragmatic option that is frequently adopted in the semiparametric literature on cure rate models [2]. While not strictly required, it has several appealing consequences for finite-sample identifiability and estimation stability, both of which can otherwise be problematic even when the model is formally identified [35].

Our methodology can also accommodate the weaker assumption that there exists some \(\tau > 0\) such that the hazard is constant beyond \(\tau\), i.e., \(\widetilde{h}(t \mid a, x) = c(a, x)\) for some \(c(\cdot, \cdot)\) and \(t > \tau\), which also identifies the causal effects. This replaces the sharp cutoff beyond \(\tau\) with an exponential tail. See [2] for a discussion of this point. The key requirement for identifiability is that the tail behavior of \([T_i \mid T_i < \infty, A_i, X_i]\) is sufficiently constrained for the model to determine where the survival curves level off.

2.2 The Latency and Stochastic Latency Effects↩︎

The effects \(\Delta_S(t)\) and \(\Delta_R(t)\) are standard causal estimands in survival analysis, whereas \(\Delta_C\) is specific to cure rate models. Defining a companion causal effect to \(\Delta_C\) that removes the contribution of cure is desirable, but doing so is not straightforward and involves potential pitfalls. For example, the uncured average treatment effect (UATE) proposed by [30], \(\Delta_{\text{uncured}}(x) = \mathbb{E}\{\log T_i(1) \mid X_i = x,\, T_i(1) < \infty\} - \mathbb{E}\{\log T_i(0) \mid X_i = x,\, T_i(0) < \infty\}\), does not correspond to a standard causal estimand because it compares potential outcomes across two different latent populations: individuals who would remain uncured under treatment and individuals who would remain uncured under control. [30] discuss further assumptions required for \(\Delta_{\text{uncured}}\) to admit a causal interpretation.

As a first step in the direction of defining a causal effect that “removes” the effect of the cure process from RMST we define the latency effect \[\begin{align} \Delta_L(t) = \mathbb{E}\{L_i(1, t) - L_i(0, t)\} \quad \text{where} \quad L_i(a, t) = R_i(a, t) - t \, D_i(a, \infty). \end{align}\] Equivalently, one can show that $ {L_i(a, t)} = _0^t (u < T_i(a) < ) , du,$ so that \(\Delta_L(t) = \int_0^t \big[\Pr(u < T_i(1) < \infty) - \Pr(u < T_i(0) < \infty)\big] \, du\) represents the integrated difference in the probability of being alive and uncured at time \(u\), for \(u \in [0, t]\). This implies that \(\Delta_L(t)\) measures how much more time, on average up to horizon \(t\), a treated individual spends alive and uncured relative to an untreated individual.

The value of defining the latency effect is that it allows for decomposition of \(\Delta_R(t)\) into a part that is directly attributable to the probability of being cured and a part that is not; in particular, \(\Delta_R(t) = \Delta_L(t) + t \, \Delta_C\). Moreover, \(\Delta_L(t)\) is a valid causal effect in the sense that it contrasts two well-defined potential outcomes.

Unfortunately, \(\Delta_L(t)\) does not have a clean clinical interpretation. We investigate this by introducing the principal strata [4] \[\begin{align} UU &= \{i : T_i(0) < \infty, T_i(1) < \infty\}, & UC &= \{i : T_i(0) < \infty, T_i(1) = \infty\}, \\ CU &= \{i : T_i(0) = \infty, T_i(1) < \infty\}, & CC &= \{i : T_i(0) = T_i(1) = \infty\}. \end{align}\] A common assumption in this setting is that the treatment is beneficial in the sense that the stratum \(CU\) is empty, i.e., nobody who died under treatment would have survived under control. Under this assumption, we can express \(\Delta_L(t)\) in terms of its contribution to the causal effect on RMST in \(UU\). The following proposition relates these two quantities; its main implication is that individuals who are cured by treatment contribute \(-\Delta_C \, \mathbb{E}\{R_i(0,t) \mid i \in UC\}\) to \(\Delta_L(t)\), which can make the treatment appear harmful (as measured by \(\Delta_L(t)\)) simply because the treatment is effective at curing individuals.

All proofs are deferred to the Supplementary Material. This flaw in \(\Delta_L(t)\) motivates us to find an alternate decomposition of \(\Delta_R(t)\) that does not have this flaw. Let \(p_a(x) = \Pr\{T_i(a) < \infty \mid X_i = x\}\) denote the probability of not being cured under treatment \(a\), and define the conditional RMST among uncured individuals as \(m_a(t,x) = \mathbb{E}\{R_i(a,t) \mid T_i(a) < \infty, X_i = x\}\). Equivalently, if we let \(G_a(\cdot \mid x)\) denote the survival function of \(T_i(a)\) conditional on not being cured, then \(m_a(t,x) = \int_0^t G_a(u \mid x) \;du\). With these definitions, we have \(\mathbb{E}\{R_i(a,t) \mid X_i = x\} = \{1 - p_a(x)\} \times t + p_a(x) \times m_a(t,x)\).

Next we define a synthetic conditional distribution \(Q_{aa'}(\cdot \mid x)\) by combining the cure probability under treatment \(a\) with the finite-event-time distribution under \(a'\) as \(Q_{aa'}(\cdot \mid x) = \{1 - p_a(x)\} \, \delta_\infty(\cdot) + p_{a}(x) \, G_{a'}(\cdot \mid x)\). Let \(T_i(a, a') \sim Q_{aa'}(\cdot \mid X_i)\) denote a draw from this distribution conditional on \(X_i\) and define \(R_i(a,a',t) = \min\{T_i(a,a'), t\}\). The corresponding stochastic-intervention RMST is $ _t(a,a’) = {R_i(a,a’,t)} = [{1 - p_a(X_i)} t + p_a(X_i) m_{a’}(t, X_i)]$ and we define the stochastic cure component of the RMST effect as \[\begin{align} \Delta_{SC}(t) &= \frac{\vartheta_t(1,0) - \vartheta_t(0,0)}{2} + \frac{\vartheta_t(1,1) - \vartheta_t(0,1)}{2} \\ &= \mathbb{E}\left[ \{p_0(X_i) - p_1(X_i)\} \left\{ t - \frac{m_0(t,X_i) + m_1(t,X_i)}{2} \right\} \right]. \end{align}\] The intuition behind this definition is that we imagine that our intervention can separately be used to replace \(p_1(x)\) with \(p_0(x)\) (i.e., modifying the effect of the cure) and to replace \(G_1(\cdot \mid x)\) with \(G_0(\cdot \mid x)\) (i.e., modifying the effect on survival among those who are uncured). It is stochastic in the sense that the effect can be expressed as \(\mathbb{E}\{\vartheta_t(1, A') - \vartheta_t(0, A')\}\) where \(A' \sim \operatorname{Bernoulli}(1/2)\), so that we randomly choose which survival distribution we assign.

The stochastic cure effect has a complementary latency component, defined as \[\begin{align} \Delta_{SL}(t) = \mathbb{E}\left[ \{m_1(t, X_i) - m_0(t, X_i)\} \left\{\frac{p_0(X_i) + p_1(X_i)}{2} \right\} \right]. \end{align}\] Similarly, it can be expressed as \(\Delta_{SL}(t) = \mathbb{E}\{\vartheta_t(A',1) - \vartheta_t(A',0)\}\), where \(A' \sim \operatorname{Bernoulli}(1/2)\). We write \(\Delta_{SC}(t,x)\) and \(\Delta_{SL}(t,x)\) for the corresponding conditional quantities. As these effects are functionals of the observed data distribution, both are identified under Assumptions 1—5. Moreover, they decompose the RMST effect.

The proposed stochastic intervention effects attempt to address the main interpretational difficulty of \(\Delta_L(t)\), which implicitly attributes the entire \(t\) units of RMST in a cured individual to the effect of cure. By contrast, \(\Delta_{SC}(t)\) treats the gain from cure as the difference between \(t\) and a reference finite-event-time RMST. Thus, if treatment increases the cure probability but the finite-event-time distribution is otherwise unchanged, \(\Delta_{SC}(t)\) captures the full RMST effect, while conversely if the treatment does not change the cure probability then \(\Delta_{SC}(t) = 0\) and the entire RMST effect is attributed to changes in the finite-event-time distribution.

Without further assumptions, we view \(\Delta_{SC}(t)\) and \(\Delta_{SL}(t)\) as model-based standardizations of the counterfactual survival distributions rather than as literal mechanistic effects. A mechanistic interpretation would require one to imagine treatment components that can separately alter cure incidence and post-susceptibility latency, bringing the estimands close in spirit to separable effects or interventional mediation effects [36], [37]. To link the stochastic effects to more mechanistic effects, \(\Delta_{SL}(t)\) can also be related to effects in the principal strata.

Relative to Proposition 1, \(\Delta_{SL}(t,x)\) is, informally, a “doubly-robust” proxy for \(\Delta_{UU}(t,x)\): it can be proportional to \(\Delta_{UU}(t,x)\) by either (i) having \(\Delta_C(x) = 0\) (i.e., there is no cure effect) or (ii) by satisfying \(B_{UU}(t,x) = B_{UC}(t,x)\). The latter is a principal ignorability condition stating that, for individuals who are uncured under control, knowing whether they would be cured under treatment does not inform their RMST under control. This type of principal ignorability is discussed by [3]. Consequently, the relationship in Proposition 3 is much more desirable than the relationship in Proposition 1. Proposition 3 might also form the basis of a sensitivity analysis by varying the sensitivity parameter \(\rho(x) = B_{UU}(t,x) - B_{UC}(t,x)\) over some plausible range of values, although we do not pursue this.

2.2.0.1 Share of RMST Attributable to Cure

A useful summary of the overall effect of the cure probability is the signed cure share \(\Delta_{SC}(t) / \Delta_R(t)\), provided that \(\Delta_R(t)\) is not close to zero. Because the two components may have opposite signs, this ratio need not lie between zero and one. If a bounded descriptive measure is desired, one can instead report \(|\Delta_{SC}(t)| / \{|\Delta_{SC}(t)| + |\Delta_{SL}(t)|\}\), which can be interpreted as the relative magnitude of the cure component rather than as a percentage of the net effect.

3 Bayesian Machine Learning for Cure Rate Models↩︎

To estimate the causal effects described in Section 2, we propose machine learning methods based on the promotion-time cure model [38] \[\begin{align} \label{eq:promotion-time} S(t \mid a, x) = \exp\left\{ -\theta(a,x) \, \widetilde{F}(t \mid a,x) \right\}, \end{align}\tag{1}\] where \(\widetilde{F}(t \mid a, x)\) is a cumulative distribution function. Because \(\lim_{t \to \infty} \widetilde{F}(t \mid a, x) = 1\), this implies that the cure probability is given by $ _{t } S(t a, x) = e^{-(a,x)}.$ The form of the promotion-time model is motivated by assuming a latent number of competing risks \(N_i(a) \sim \text{Poisson}(\theta(a,x))\) for individual \(i\) under treatment \(a\). If \(N_i(a) = 0\), the individual is cured and \(T_i(a) = \infty\); otherwise, each of the \(N_i(a)\) latent failure times are drawn i.i.d. from the distribution with distribution function \(\widetilde{F}(t \mid a, x)\), and the survival time \(T_i(a)\) is the minimum of these latent failure times (see [7] for a derivation). An advantage of the promotion-time model is that it more effectively borrows information across the cure probability \(e^{-\theta(a,x)}\) and the survival function conditional on not being cured $ G_a(t x) = {e^{-(a,x) , (t a, x)} - e^{-(a,x)}} / {1 - e^{-(a,x)}}\(, because both quantities are determined by the same underlying functions\)(a,x)$ and \(\widetilde{F}(t \mid a,x)\). By contrast, nonparametric mixture cure models such as the one used by [30] parameterize the mixture probability separately from the survival distribution, and consequently do not share information across these components.

3.1 A Brief Review of BART↩︎

To model the unknown functions \(\theta(a,x)\) and \(\widetilde{F}(t \mid a, x)\) we use Bayesian additive regression trees (BART) models. A general BART model expresses a function of interest as a sum of regression trees $ r(a, x) = _{m=1}^M g(a, x; _m, _m),$ where \(\mathcal{T}_m\) denotes the tree shape and splitting rules of tree \(m\) and \(\mathcal{M}_m = \{\mu_{\ell m} : \ell \in \operatorname{Leaves}(\mathcal{T}_m)\}\) consists of the leaf node predictions such that \(g(a, x ; \mathcal{T}_m, \mathcal{M}_m) = \mu_{\ell m}\) if \((a, x)\) is associated to leaf node \(\ell\) of \(\mathcal{T}_m\). Each \(\mathcal{T}_m\) is a binary tree whose internal (non-terminal) nodes are equipped with splitting rules of the form \([x_j \le c]\), where \(x_j\) is a component of the input vector \((a, x)\) and \(c\) is a cutpoint. An input arriving at an internal node is associated to the left child if the splitting condition is satisfied, and to the right child otherwise, until a terminal (leaf) node is reached.

BART was introduced by [9] for nonparametric regression with normal errors, as well as probit regression. In causal inference, [39] used BART to estimate heterogeneous treatment effects in observational studies, with the main goal of flexibly adjusting for confounding and recovering individualized counterfactual outcome surfaces without having to prespecify interactions or nonlinear main effects. Extending this, [31] proposed the Bayesian causal forest (BCF), which reparameterizes the outcome model into separate prognostic and treatment-effect components so that treatment effect heterogeneity can be regularized more aggressively.

3.2 The One-Forest and Two-Forest Models↩︎

We now describe the BartCuremodel for estimation of the causal estimands defined in Section 2. We introduce two variants of BartCure, which we refer to as the one-forest and two-forest models, respectively. These models are fully nonparametric, in the sense that they can approximate any conditional survival function \(S(t \mid a, x)\). The two-forest approach directly applies the BART model of [9] to 1 by setting \[\begin{align} \log \theta(a,x) = \sum_{m = 1}^M g(a, x; \mathcal{T}^\theta_m, \mathcal{M}^\theta_m) \;\;\text{and} \;\; \widetilde{h}(t \mid a,x) = \lambda_0(t) \exp\left\{ \sum_{m = 1}^M g\big(B(t), a, x; \mathcal{T}^h_m, \mathcal{M}^h_m\big) \right\}, \end{align}\] where \(B(\cdot): \mathbb{R}^+ \to \{1, \ldots, K\}\) maps a given time \(t\) to a bin in which the hazard function is constant, \(\lambda_0(t)\) is a baseline hazard function, and \(\widetilde{h}(t \mid a, x) = -\frac{d}{dt} \log \{1 - \widetilde{F}(t \mid a, x)\}\) is the hazard function associated to \(\widetilde{F}(t \mid a, x)\). We adopt the shorthand \(r^\theta(a,x) = \log \theta(a,x)\) and \(r^h(t, a, x) = \sum_m g\big(B(t), a, x; \mathcal{T}_m^h, \mathcal{M}_m^h \big)\).

The one-forest model instead reparameterizes the model as \[\begin{align} \label{eq:one-forest} \theta(a,x) \, \widetilde{f}(t \mid a, x) = \lambda_0(t) e^{r(t, a, x)} \quad \text{where} \quad r(t, a, x) = \sum_{m = 1}^M g\big( B(t), a, x ; \mathcal{T}_m, \mathcal{M}_m \big). \end{align}\tag{2}\] and where \(\widetilde{f}(t \mid a, x) = \frac{\partial}{\partial t} \widetilde{F}(t \mid a, x)\) is the conditional density of \(\widetilde{F}(t \mid a, x)\). From this, we can recover \(\theta(a,x) = H(\infty \mid a, x)\) and \(\widetilde{F}(t \mid a, x) = H(t \mid a, x) / H(\infty \mid a, x)\) where \(H(t \mid a, x) = \int_0^t \lambda_0(u) e^{r(u,a,x)} \;du\). The one-forest model is less faithful to the promotion-time parameterization, but retains its important benefits: it is computationally simple, and allows the cure probability to be a priori correlated with better survival prognosis. A potential drawback of the one-forest model is that sharing information between \(\theta(a,x)\) and \(\widetilde{f}(t \mid a, x)\) may not always be desirable.

3.2.0.1 Prior on \(\mathcal{T}_m\)

All of the tree structures are assigned a branching process prior as described by [9], [40]. Each node at depth \(d\) is non-terminal with probability $ ,$ where \(\gamma \in (0, 1)\) and \(\beta > 0\) control the tree size. The default values \(\gamma = 0.95\) and \(\beta = 2\) encourage trees with few terminal nodes, so that each tree captures only a small piece of the overall signal. Because treatment effect heterogenity only arises through splitting on \(a\) and one other variable, this prior naturally penalizes treatment effect heterogeneity.

3.2.0.2 Prior on \(\mathcal{M}_m\)

To facilitate computations, log-gamma priors are specified for the leaf node predictions \(\mu_{\ell m } \sim \log \operatorname{Gam}(a_\mu, b_\mu)\). In this setting, as well as other survival settings, the log-gamma prior is conditionally conjugate [25], [26], [29]. Following [29], we choose \(a_\mu\) and \(b_\mu\) so that \(\mathbb{E}(\mu_{\ell m}) = 0\) and \(\operatorname{Var}(\mu_{\ell m}) = \sigma^2_\mu / M\) with respect to the prior, with \(\sigma_\mu = 1.5\) in all of our examples.

3.2.0.3 Prior for \(\lambda_0(t)\)

For both the one-forest and two-forest models, we model the baseline hazard \(\lambda_0(t)\) as piecewise-constant with cut points \(0 = \alpha_0 < \alpha_1 < \cdots < \alpha_{K - 1} < \alpha_K = \infty\) such that \(\lambda_0(t) = \sum_{k = 1}^K \lambda_k \times 1(\alpha_{k-1} \le t < \alpha_k)\). We impose Assumption 5 by setting \(\alpha_{K -1} = \tau\) and \(\lambda_K = 0\), although this is not strictly required to fit the two-forest model. We then set \(\lambda_k \stackrel{\textrm{iid}}{\sim}\operatorname{Gam}(a_\lambda, b_\lambda)\) with \(a_\lambda = 1\) and \(b_\lambda \sim \operatorname{Gam}(1, 1)\) by default.

3.2.0.4 Targeted Selection

In observational studies, treatment assignment may depend strongly on the baseline prognosis, a phenomenon referred to as targeted selection [31]. To mitigate the impact of targeted selection, it is recommended to include an estimate of the propensity score as a covariate in the model. We follow this approach by incorporating an estimated propensity score into both the one-forest and two-forest BartCuremodels.

3.2.0.5 Computation

Inference for both the one-forest and two-forest models is performed via Gibbs sampling using the generalized Bayesian backfitting algorithm of [41]. In the case of the one-forest model, the algorithm is essentially identical to the algorithm of [26] except that the cure threshold is taken into account. A Gibbs sampler for the two-forest model is not immediate from prior work, but is feasible after augmenting the latent number of unobserved events \(N_i \sim \operatorname{Poisson}\{\theta(A_i, X_i) \, [1 - \widetilde{F}(Y_i \mid A_i, X_i)]\}\). Importantly for computational feasibility, the latent events themselves do not need to be sampled. Details are deferred to the Supplementary Material.

4 Experiments↩︎

We now conduct simulation experiments to evaluate the performance of the one-forest and two-forest BartCuremodels relative to baseline methods. Our goal is to determine whether BartCureprovides reliable estimation of average causal effects, conditional average causal effects, and treatment effect heterogeneity, in the presence of a cured subpopulation.

4.1 Estimating Causal Effects↩︎

We compare BartCureto competing methods on a suite of data generating processes (DGPs) that were introduced in prior work on causal survival analysis, augmented (or not) to have a cured subpopulation.

4.1.0.1 Methods for Comparison

We compare the one-forest model and two-forest models, which we label BartCure(1) and BartCure(2). We also consider two existing causal machine learning survival analysis methods: causal survival forests (CSF, [16]) and an accelerated failure time BART model (IndivAFT, [15]). CSF is implemented in the grf package [42], while IndivAFT is implemented in the AFTrees package. IndivAFT sets \(\log T_i(a) = \mu_a(X_i) + \epsilon_i\) with a flexible Dirichlet process mixture model specified for the error distribution. To address targeted selection [31] all BART methods incorporated the estimate propensity score as a covariate, while the CSF uses the propensity score via inverse propensity score weighting. All BART methods use 200 trees and inferences are based on 2,000 draws from a Gibbs sampler targeting the posterior.

4.1.0.2 Data Generating Processes

We consider three families of DGPs drawn from the causal survival analysis literature [15], [16], [43]. The same DGPs were analyzed by [12]; full details on these DGPs are in the Supplementary Material. All experiments used \(N = 1000\) samples, which is somewhat smaller than the sample size of the CALGB 40101 breast cancer trial. Each of these DGPs was extended to include a cured fraction via the mixture cure model $ S(t a, x) = {1 - p_a(x)} + p_a(x) G_a(t x)$ where \(G_a(t \mid x)\) is the survival function of the uncured sub-population, obtained by truncating the base survival at the time horizon \(t\). We set the cure probability to \(\operatorname{expit}\{\alpha + z_a(x)\}\) where \(z_a(x)\) is \(\operatorname{logit}G_a(t \mid x)\) standardized to have mean \(0\) and variance \(1\), and \(\alpha\) is chosen so that the marginal cure rate is \(1/2\) under all settings. This implies that individuals with a high probability of surviving until the follow-up horizon \(t\) in the uncured population also have a high probability of being cured.

4.1.0.3 Comparison Measures

For the simulation settings without a cured subpopulation, we estimate \(\Delta_R(t), \Delta_R(t, x)\), \(\Delta_S(t)\), and \(\Delta_S(t, x)\) at the fixed time horizon \(t\). For settings with a cured subpopulation, we monitor \(\Delta_C\) and \(\Delta_C(x)\) rather than \(\Delta_S(t)\) and \(\Delta_S(t,x)\). We then compute the bias, root mean-squared error (RMSE), and the length and coverage of nominal 95% confidence/credible intervals, averaged over all simulated datasets and individuals.

4.2 Results for Estimating Causal Effects↩︎

Table 1 and Table 2 report the results for the causal effects in settings with and without a cured subpopulation, respectively. We defer simulation results for conditional average treatment effects to the Supplementary Material, with results in Figures 7—12.

4.2.0.1 Average Treatment Effects: With a Cured Subpopulation

Across all settings, methods perform similarly for estimating average treatment effects across all metrics, with the exception of BartCure(2), which appears less stable specifically for estimating \(\Delta_C\) and so has wider intervals. Aside from this, both IndivAFT and BartCure(2)perform quite poorly on the Hu 3 and Hu 4 DGPs for estimating \(\Delta_C\), and IndivAFT performs poorly for \(\Delta_R(t)\) as well; in the case of IndivAFT, this occurs because the model does not account for the cured subpopulation, and it also performs poorly for estimating \(\Delta_C\) on the Henderson 2 setting. In view of this, the best performing methods for estimating ATEs are BartCure(1) and the CSF; BartCure(1)performs marginally better by RMSE across most settings and better by coverage on the Hu DGPs, and so performs best overall.

4.2.0.2 Average Treatment Effects: Without a Cured Subpopulation

In the absence of a cured subpopulation, the basic findings are the same, with the exception that IndivAFT no longer performs poorly on any of the settings. We see that it is extremely difficult for any method to estimate \(\Delta_S(t)\) under the Hu 3 and Hu 4 settings. This occurs because the survival functions are very close to zero at these follow-up times, and so small amounts of bias in estimating \(\Delta_S(t)\) dominate inference; coverage for \(\Delta_R(t)\), which is determined by the full survival curve, is much better. Interestingly, CSF suffered from substantial instability across all Henderson settings due to its use of inverse weighting highlighting that the BART-based methods can be much more stable.

4.3 Directionality of Heterogeneity↩︎

To complement our results on estimation of individual effects, we also give an assessment of how well the one-forest BartCuremodel and CSFs perform on the task of estimating the directionality of an effect. To test for the existence of individual differences in effect, we constructed for each individual an 80% confidence/credible interval and checked whether this interval contains \(0\). For each individual, we then computed the centered treatment effect \(\Delta_R(\tau, x) - \Delta_R(\tau)\) so that positive values indicate above-average benefit. We recorded the discovery rate, the correct sign rate, the Type-S error rate, and the net directional score, \[\begin{align} \text{DR}= \frac{N_D}{N}, \quad \text{CS}= \frac{N_{CD}}{N}, \quad \text{Type-S} = \frac{N_{ID}}{N_D} 1(N_D \ne 0), \quad \text{NDS}= \frac{N_{CD} - N_{ID}}{N} \end{align}\] where \(N\) is the sample size, \(N_D\) is the number of discoveries, \(N_{CD}\) is the number of correct discoveries, and \(N_{ID}\) is the number of incorrect discoveries. A method that discovers few effects will have a low correct sign rate but also a low Type-S error rate, while a method that declares many effects aggressively may achieve a high correct sign rate at the cost of elevated Type-S errors; the net directional score summarizes the tradeoff. Higher values of the net directional score indicate that a method is both willing and able to recover the direction of treatment-effect heterogeneity.

Table 3 shows that BartCureis generally more effective at recovering the direction of RMST heterogeneity. In the Henderson settings, BartCuremakes substantially more discoveries than CSF, with discovery rates around \(35\%\) compared to around 6%, while keeping Type-S error rates small. This leads to much larger net directional scores for BartCureacross all four Henderson DGPs. The same pattern holds for Cui 2 and Hu 3, where BartCurehas both a higher correct sign rate and lower Type-S error rates than CSF.

The main exceptions to this trend occur at Cui 1, Hu 1, and Hu 4. For Cui 1, both methods make very few discoveries, and the net directional scores are close to zero, indicating that neither method reliably detects directional heterogeneity; given the small amount of treatment effect heterogeneity in this setting, both methods appear to be behaving reasonably in failing to detect heterogeneity across the board. For Hu 1 and Hu 4, CSF is slightly more aggressive and obtains marginally larger net directional scores, although this comes with higher Type-S error rates. Overall, these results suggest that BartCureis conservative in weak-signal settings but, when heterogeneity is identifiable, it tends to recover directional effects with a favorable balance between discovery and sign accuracy.

5 Application to CALGB 40101 Trial↩︎

We now apply BartCureto data from 3,871 participants in the CALGB 40101 trial, a Phase III randomized trial designed to evaluate the noninferiority of paclitaxel (T) relative to cyclophosphamide+doxorubicin (CA) for the treatment of breast cancer. Our primary endpoint is disease-free survival (DFS). Covariates of interest include age, race, ethnicity, menopause status (stra1), hormone receptor status (stra2), tumor size, assigned treatment duration, and number of positive lymph nodes. We restrict attention to 3,864 subjects with complete covariate information. As suggested by [2], we set \(\tau\) in Assumption 5 to just-after the last observed failure time (\(\tau = 115\) months).

5.0.0.1 Average Effects and Model Comparison

Table 4 compares average causal effect estimates from different models for RMST and cure probability. Because treatment assignment is randomized, estimates based on the Kaplan–Meier estimator (displayed in Figure 1) provide a useful benchmark. The Kaplan–Meier-based approach estimates a \(-2.39\)-month RMST difference (95% CI: \(-4.31\), \(-0.47\)) and a \(-4.7\%\) cure-probability difference (95% CI: \(-10.7\%\), \(1.3\%\)). BartCureis closest to this benchmark for \(\Delta_R(115)\) (\(-2.24\) months), with intervals excluding zero for both \(\Delta_R(115)\) and \(\Delta_C\); this is sensible because, while BartCureis nonparametric, it applies some amount of shrinkage to the treatment effect. CSF estimates an RMST contrast in the same direction and of slightly larger magnitude, but with much wider intervals that include zero; its cure-probability estimate is also directionally similar but imprecise. IndivAFT gives the most attenuated estimates for both RMST (\(-1.94\) months) and cure probability (\(-2.8\%\)), with comparatively narrow intervals; this may be due either to the imposition of the accelerated failure time assumption or to the fact that IndivAFT does not account for the existence of a cured population. Overall, the model-based methods agree on the direction of the effect, but BartCuremost closely reproduces the nonparametric RMST estimate.

The posterior distributions of \(\Delta_R(t)\), \(\Delta_{SL}(t)\), and \(\Delta_C\) for BartCureare given in Figure 2. Because \(\Delta_{SL}(t)\) is indistinguishable from \(0\), the RMST effect appears primarily attributable to T curing fewer individuals rather than to latency differences among uncured individuals

5.0.0.2 Individual Effects

The posterior distributions of \(\Delta_R(115, X_i)\) and \(\Delta_C(X_i)\) for the individuals in the sample are given in Figure 3; individuals are ordered by the posterior mean effect, and bands show pointwise 95% credible intervals for the effects. We see that RMST and cure probability are very closely related. All individual effect estimates favor CA over T, with differing amounts of supporting evidence. Individual posterior probabilities of negative effects of T are given in Table 5; we see that roughly 75% of individuals are estimated to have a posterior probability of a negative treatment effect of at least 80%, and around 18% have at least a 95% posterior probability of a negative effect.

For comparison, boxplots of individual RMST effect estimates for BartCureand CSF are given in Figure 4. We see that CSFs produce much larger estimates of treatment effect heterogeneity, with many individuals estimated to have a positive T-versus-CA effect. By contrast, BartCureregularizes towards a nearly-homogeneous treatment effect.

5.0.0.3 Subgroup Effects

Following [31], [44], we constructed a decision tree summary to find optimal decision-theoretic subgroups that maximize treatment effect heterogeneity. Figure 5 provides a decision tree summary of \(\Delta_R(t,x)\). This summary identifies two variables that the treatment may differ across: the age of an individual (which was stratified into 6 groups, with groups 5 and 6 corresponding to individuals above the age of 60), and stra2 (which stratifies on hormone receptor status, with 1 for positive and 2 for negative). Generally, older individuals with negative receptor status have more negative RMST contrast estimates across arms.

Figure 6 shows posterior distributions for \(\Delta_R(t)\) and \(\Delta_C\) within these subgroups, centered by the population average treatment effect. For all subgroups, there is at-most weak evidence supporting treatment effect heterogeneity, with sizeable uncertainty. We see that older individuals have an RMST difference estimated to be roughly 1.5 months lower and a cure probability effect about 2.5% lower than younger individuals. We also see a negative effect of stra2 = 2, and both of these effects are roughly additive when looked at jointly.

5.0.0.4 Additional Results

Further analysis of this dataset is available in the Supplementary Material. This includes an assessment of variable importance and analysis of the proportion of the RMST effect attributable to cure.

6 Discussion↩︎

In this work, we proposed a causal machine learning framework for survival settings with a cured subpopulation. We defined causal effects on finite-time survival, restricted mean survival time (RMST), and cure probability, and considered decompositions of the RMST effect into contributions arising from cure and from delayed failure among uncured individuals. This cure-rate formulation provides clinically relevant information that standard causal analyses do not directly capture. We first considered a naive latency contrast, which is a well-defined causal estimand, but argued that its interpretation is unsatisfactory. To address this issue, we introduced stochastic cure and latency effects that separate the contribution of the cure probability from that of the finite-event-time distribution. We also showed how the stochastic latency effect relates to a principal-strata effect among individuals who would remain uncured under either treatment arm. For estimation, we proposed BartCure, a BART-based promotion-time cure model targeting both average and conditional versions of these effects.

The main strength of BartCureis its ability to estimate clinically interpretable cure-rate causal effects while retaining the flexibility of Bayesian nonparametric regression. In simulations, BartCureperforms well relative to causal survival forests and an accelerated failure time BART model, both in settings with a cured subpopulation and in settings where all survival times are finite. A second strength is its behavior when estimating heterogeneous treatment effects. In many applications, overstating heterogeneity can be more harmful than failing to detect weak heterogeneity, since subgroup claims often influence clinical interpretation and the design of future studies. The prior used by BartCureshrinks treatment effects toward homogeneity, and the simulation results suggest that this induces conservative heterogeneity detection with favorable Type-S error properties. Consequently, when the method identifies strong subgroup patterns, those patterns have persisted despite the model’s regularization toward homogeneous effects.

The proposed estimands should be interpreted with care. The stochastic latency effect is identified as a functional of the observed data distribution under the stated assumptions, but connecting it to the principal-strata effect \(\Delta_{UU}\) requires additional conditions; the connection is most direct under monotonicity and a principal ignorability condition relating the control-arm RMST distribution in the \(UU\) and \(UC\) strata. These assumptions may be plausible in some clinical settings, but they are not fully testable. A second limitation is the reliance on a cure-threshold assumption. This assumption is common for semiparametric cure models, and is useful for stabilizing estimation of the cured fraction, but it remains a substantive modeling choice. If the threshold is chosen too small, we may incorrectly classify uncured individuals as cured; if it is chosen too large, the cure probability may be weakly identified in small samples. In practice, the threshold should therefore be guided by clinical knowledge and sensitivity analyses should be conducted whenever possible.

6.0.0.1 Acknowledgements

This work was supported under NSF grant DMS-2144933.

6.0.0.2 Supplementary Material

Supplementary files include formal identification assumptions, proofs of all propositions, full details of the Gibbs sampling algorithm used to fit the models, additional simulation details and experiments, additional analyses CALGB 40101, and code reproducing figures and tables.

7 Causal Identification Assumptions↩︎

We require the following assumptions to identify the causal parameters. Assumptions 1–3 are standard assumptions in the causal inference literature [32], which state that the treatment is well-defined with no interference between units, all individuals have some non-negligible probability of being assigned to either treatment or control, and there are no unmeasured common causes of the treatment assignment mechanism and the potential outcomes; in CALGB 40101, the treatment assignment was randomized, and so these assumptions are known to hold. Assumption 4 is also standard in survival analysis, and is plausible when censoring is primarily administrative.

Assumption 1:

Stable Unit Treatment Value (SUTVA). The observed failure time satisfies \(T_i = T_i(A_i)\), and the potential outcomes of individual \(i\) do not depend on the treatment assignments of other individuals [33].

Assumption 2:

Positivity. There exists a \(\delta > 0\) such that \(\delta < e(x) < 1 - \delta\) for all \(x\) in the support of \(X_i\), where we recall that \(e(x) = \Pr(A_i = 1 \mid X_i = x)\) is the propensity score.

Assumption 3:

Ignorability. Treatment assignment is independent of the potential outcomes given the observed covariates: $ {T_i(0), T_i(1)} A_i X_i.None$

Assumption 4:

Ignorable Censoring. The censoring mechanism is independent of the observed failure times given the treatment and covariates: $ C_i T_i A_i, X_i.None$

Assumption 5:

Cure Threshold. There exists a known time \(\tau < \infty\) such that the hazard function satisfies \(h(t \mid a, x) = 0\) for all \(t > \tau\), \(a \in \{0,1\}\), and \(x\) in the support of \(X_i\). Equivalently, survival beyond \(\tau\) implies cure. Additionally, we require sufficient follow-up in the sense that \(\Pr(C_i > \tau \mid X_i = x) > 0\) for all \(x\).

8 Proofs of Propositions↩︎

8.0.0.1 Proof of Proposition 1.

Because \(L_i(a,t) = R_i(a,t)1\{T_i(a) < \infty\}\), \[\begin{align} \Delta_L(t) &= \mathbb{E}\left[R_i(1,t)1\{T_i(1) < \infty\} - R_i(0,t)1\{T_i(0) < \infty\}\right]. \end{align}\] The condition \(\Pr(i \in CU \mid X_i) = 0\) implies \(\Pr(i \in CU)=0\), so \(1\{T_i(1)<\infty\}=1\{i \in UU\}\) and \(1\{T_i(0)<\infty\}= 1\{i \in UU \cup UC\}\). Hence \[\begin{align} \Delta_L(t) &= \Pr(i \in UU)\mathbb{E}\{R_i(1,t)-R_i(0,t)\mid i \in UU\} - \Pr(i \in UC)\mathbb{E}\{R_i(0,t)\mid i \in UC\}. \end{align}\] Under the same condition, \(p_1=\Pr(i\in UU)\) and \(\Delta_C=\Pr\{T_i(1)=\infty\}-\Pr\{T_i(0)=\infty\} =\Pr(i\in UC)=p_0-p_1\). Substitution gives the result.

8.0.0.2 Proof of Proposition 2.

Fix \(x\). By definition, \[\begin{align} \Delta_R(t,x) &= \{1-p_1(x)\}t + p_1(x)m_1(t,x) - \{1-p_0(x)\}t - p_0(x)m_0(t,x) \\ &= \{p_0(x)-p_1(x)\}t + p_1(x)m_1(t,x)-p_0(x)m_0(t,x). \end{align}\] Also, \[\begin{align} \Delta_{SC}(t,x)+\Delta_{SL}(t,x) &= \{p_0(x)-p_1(x)\}\left\{t-\frac{m_0(t,x)+m_1(t,x)}{2}\right\} \\ &\quad+ \{m_1(t,x)-m_0(t,x)\}\left\{\frac{p_0(x)+p_1(x)}{2}\right\} \\ &= \{p_0(x)-p_1(x)\}t + p_1(x)m_1(t,x)-p_0(x)m_0(t,x), \end{align}\] which equals \(\Delta_R(t,x)\).

8.0.0.3 Proof of Proposition 3.

First, note that because \(\Pr(i \in CU \mid X_i = x) = 0\), an individual who is uncured under treatment must belong to the stratum \(UU\). Hence $ {T_i(1)<X_i=x} = (iUUX_i=x),$ so that \(\Pr(i\in UU\mid X_i=x)=p_1(x)\). Similarly, an individual who is uncured under control belongs to either \(UU\) or \(UC\), and therefore $ p_0(x) = (iUUX_i=x) + (iUCX_i=x).$ It follows that \(\Pr(i\in UC\mid X_i=x) = p_0(x)-p_1(x) = \Delta_C(x)\).

Now consider the conidtional RMST among uncured individuals. Because \(T_i(1) < \infty\) is equivalent to \(i \in UU\), we have \[\begin{align} m_1(t,x) = \mathbb{E}\{R_i(1,t)\mid i\in UU, X_i=x\}. \end{align}\] On the other hand, rearranging the expression \(\mathbb{E}\{R_i(a,t) \mid X_i = x\} = \{ 1 - p_a(x)\} \times t + p_a(x) \times m_a(t,x)\) with \(a = 0\) and applying \(\mathbb{E}\{R_i(0,t) \mid X_i = x\} = \Pr(i \in UU \mid X_i = x) \times \mathbb{E}\{R_i(0,t) \mid i \in UU, X_i = x\} + \Pr(i \in UC \mid X_i = x) \times \mathbb{E}\{R_i(0,t) \mid i \in UC, X_i = x\}\) gives \[\begin{align} m_0(t,x) &= \frac{ p_1(x)\mathbb{E}\{R_i(0,t)\mid i\in UU, X_i=x\} + \Delta_C(x)\mathbb{E}\{R_i(0,t)\mid i\in UC, X_i=x\} }{p_0(x)}. \end{align}\] Using the notation $ B_{UU}(t,x)={R_i(0,t)iUU, X_i=x}, B_{UC}(t,x)={R_i(0,t)iUC, X_i=x},$ and $ _{UU}(t,x) = {R_i(1,t)-R_i(0,t)iUU, X_i=x},$ we can write \[\begin{align} m_1(t,x)=B_{UU}(t,x)+\Delta_{UU}(t,x) \quad \text{and} \quad m_0(t,x) = \frac{p_1(x)B_{UU}(t,x)+\Delta_C(x)B_{UC}(t,x)}{p_0(x)}. \end{align}\] Therefore, \[\begin{align} m_1(t,x)-m_0(t,x) &= B_{UU}(t,x)+\Delta_{UU}(t,x) - \frac{p_1(x)B_{UU}(t,x)+\Delta_C(x)B_{UC}(t,x)}{p_0(x)} \\ &= \Delta_{UU}(t,x) + \frac{\{p_0(x)-p_1(x)\}B_{UU}(t,x)-\Delta_C(x)B_{UC}(t,x)}{p_0(x)} \\ &= \Delta_{UU}(t,x) + \frac{\Delta_C(x)}{p_0(x)} \{B_{UU}(t,x)-B_{UC}(t,x)\}. \end{align}\] Finally, by definition, \[\Delta_{SL}(t,x) = \{m_1(t,x)-m_0(t,x)\} \frac{p_0(x)+p_1(x)}{2}.\] Substituting the preceding display gives the stated identity. If either \(\Delta_C(x)=0\) or \(B_{UU}(t,x)=B_{UC}(t,x)\), the second term in braces vanishes, so \(\Delta_{SL}(t,x)\) is proportional to \(\Delta_{UU}(t,x)\). This proves the final claim.

9 Full Details of Gibbs Sampler↩︎

For the one-forest model, the Gibbs sampler used is exactly the one described by [26], with the additional constraint that \(\lambda_K \equiv 0\) is enforced at the end of each update of the \(\lambda\)’s.

For the two-forest model, consider the promotion-time survival model with survival function \[\begin{align} S(t \mid x) = \exp\left\{ -\theta(x) \widetilde{F}(t \mid x) \right\}. \end{align}\] This model has hazard function \[\begin{align} h(t \mid x) = \theta(x) \, \widetilde{f}(t \mid x). \end{align}\] Adopt the shorthand \(S = S(t \mid x)\), \(\widetilde{F}= \widetilde{F}(t \mid x)\), \(\theta = \theta(x)\), and \(\widetilde{f}= \widetilde{f}(t \mid x)\); adding the subscript \(i\) is shorthand for plugging in \(t = Y_i\) and \(x = X_i\) so \(S_i = S(Y_i \mid X_i)\). We also write \(\widetilde{S}= 1 - \widetilde{F}\) and \(\widetilde{h}= -\frac{\partial}{\partial t} \log \widetilde{S}\).

The likelihood in the promotion-time model is \[\begin{align} \prod_i (\theta_i \widetilde{f}_i)^{\delta_i} \exp\left\{ -\int_0^{Y_i} \theta_i \, \widetilde{f}\;dt \right\} = \prod_i \left( \theta_i \, \widetilde{h}_i \, \widetilde{S}_i \right)^{\delta_i} \exp\left( -\theta_i \, \widetilde{F}_i \right). \end{align}\] Now, augment \(K_i \sim \operatorname{Poisson}(\theta_i \, \widetilde{S}_i)\). The augmented likelihood becomes \[\begin{align} \prod_i \frac{\theta_i^{K_i + \delta_i} e^{-\theta_i}}{K_i!} \times \widetilde{h}_i^{\delta_i} \widetilde{S}_i^{\delta_i + K_i}. \end{align}\] Next, we consider the model $ = r^{}(x),$ with shorthand \(r^\theta\) and \(r^\theta_i\). We also consider the model \(\widetilde{h}= \lambda_{B(t)} \exp\{r^h(B(t), x)\}\), also with shorthand \(\widetilde{h}= \lambda \, \exp(r^h)\). Both of these terms can be handled in a fully conjugate fashion using log-linear BART models with log-gamma priors. For example, we could set \[\begin{align} r^\theta = \sum_{m = 1}^M g(x; \mathcal{T}^\theta_m, \mathcal{M}^\theta_m) \end{align}\] with leaf node parameters \(\mu^\theta_{\ell m} \sim \log \operatorname{Gam}(a^\theta, b^\theta)\). Similarly, we can use the model $ r^h(B(t), x) = r^h = _{m = 1}^M g(B(t), x ; ^r_m, ^r_m)$ with leaf node parameters \(\mu^r_{\ell m}\); the interpretation here is that we have \(K_i\) censored values above \(Y_i\) in the usual piecewise exponential model, with \(Y_i\) itself either denoting an event time \((\delta_i = 1)\) so that we have observed an event exactly at \(Y_i\), or itself might correspond to censoring (so that we have \(K_i + 1\) events to consider above \(Y_i\); this shouldn’t be interpreted literally in the context of the cure model, but is just how one can interpret this component of the likelihood).

The Poisson log-linear model updates are described by [45] and, while the survival updates are also relatively straightforward given the above development, we will also outline them here.

Fix leaf \(\ell\) of tree \(m\) and define \(\eta_{ib} = r^h(b, X_i) - \mu_{\ell m}\). Define \(Z_{ib} = 1(Y_i \ge c_b) (c_b - c_{b-1}) + 1(c_{b-1} \le Y_i < c_b)(Y_i - c_{b-1})\) and \(\delta_{ib} = 1(\delta_i = 1 \wedge Y_i \in (c_{b-1}, c_b))\). Then the likelihood associated to this leaf node is \[\begin{align} \prod_{i, b : (b, X_i) \leadsto (\ell, m)} \exp\left\{ \delta_{ib} (\eta_{ib} + \mu_{\ell m}) - (\delta_i + K_i)\lambda_b \, Z_{ib} \, e^{\eta_{ib}} e^{\mu_{\ell m}} \right\}, \end{align}\] where \((b,X_i) \leadsto (\ell, m)\) denotes that \((b,X_i)\) is associated to leaf node \(\ell\) of tree \(m\). Integrating out \(\mu_{\ell m}\) against its prior gives the integrated likelihood \[\begin{align} m(\mathcal{T}_m) = \prod_{\ell} \frac{\Gamma(a_\mu + A_\ell)}{(b_\lambda + B_\ell)^{a_\mu + A_\ell}} \times \frac{b_\mu ^{a_\mu }}{\Gamma(a_\mu)} \end{align}\] where \(A_\ell = \sum_{(i,b) : (b, X_i) \leadsto (\ell, m)} \delta_{ib}\) and \(B_\ell = \sum_{(i,b) : (b, X_i) \leadsto (\ell, m)} (\delta_i + K_i) \lambda_b Z_{ib} e^{\eta_{ib}}\). This integrated likelihood can be used, as described by [41], as the basis for a Metropolis-Hastings update for \(\mathcal{T}_\ell\). Given a proposal distribution \(\Lambda(\mathcal{T}\mid \mathcal{T}')\), we accept a proposed \(\mathcal{T}' \sim \Lambda(\mathcal{T}\mid \mathcal{T}_m)\) with probability \[\begin{align} 1 \wedge \frac{m(\mathcal{T}') \, \pi_\mathcal{T}(\mathcal{T}') \, \Lambda(\mathcal{T}_m \mid \mathcal{T}')}{m(\mathcal{T}_m) \, \pi_\mathcal{T}(\mathcal{T}_m) \, \Lambda(\mathcal{T}' \mid \mathcal{T}_m)}. \end{align}\] Reasonable proposals for new trees are the grow, prune, and change moves described by [40]. After updating \(\mathcal{T}_\ell\) by Metropolis-Hastings, we then sample \(\mu_{\ell m} \sim \log \operatorname{Gam}(a_\mu + A_\ell, b_\mu + B_\ell)\).

Finally, after all of the \((\mathcal{T}_m, \mathcal{M}_m)\)’s have been updated, the baseline hazard parameters are updated similarly as $ b (a+ A_b, b_+ B_b)$ where \(A_b = \sum_i \delta_{ib}\) and \(B_b = \sum_i (\delta_i + K_i) Z_{ib} e^{r^h_i}\).

10 Additional Simulation Details↩︎

10.1 More Details on the DGPs↩︎

We now provide the full details for the different DGPs used in the simulation experiments of Section 4.

10.1.0.1 [15]

This setting is based on real data taken from the SOLVD trial [15]. Covariates and treatment data are taken directly from the trial data, while the outcome \(T_i(a)\) follows an accelerated failure time model with \(\log T_i(a) = \mu_a(x) + \epsilon_i\), where \(\epsilon_i\) has mean zero and unit variance following either a normal (DGP 1), Gumbel (DGP 2), standardized gamma (DGP 3), or \(t\)-mixture (DGP 4) distribution. The propensity score is estimated by a logistic regression on \(X_i\), while the censoring time is sampled independently as \(C_i \sim \operatorname{Uniform}(6.5, 10)\). We use \(t = 6\) for the follow-up time.

10.1.0.2 [16]

Covariates for this setting are sampled by setting \(X_i = U_i V\) where \(U_{ij} \sim \operatorname{Uniform}(0,1)\) and \(V\) is given by the Cholesky decomposition of a matrix with \((j,k)^{\text{th}}\) entry \(0.5^{|j - k|}\). We consider two sub-settings:

1. Log-normal AFT: \(\log T_i(a) = \mu_a(X_i) + \epsilon_i\) where \(\epsilon_i \sim \operatorname{Normal}(0,1)\) and with \(\mu_0(x) = -1.85 - 0.8 \cdot 1(x_1 < 0.5) + 0.7\sqrt{x_2} + 0.2 x_3\) and \(\mu_1(x) = \mu_0(x) + 0.7 - 0.4 \cdot 1(x_1 < 0.5) - 0.4\sqrt{x_2}\). The propensity score is taken to be \(e(x) = (1 + f_\beta(x_1; 2, 4))/4\), where \(f_\beta\) is the Beta density. Censoring is covariate-dependent via a Weibull model with \(C_i(a) = \sqrt{-\log(U_i)/\exp\{f_{Ca}(X_i)\}}\), where \(U_i \sim \operatorname{Uniform}(0,1)\), \(f_{C0}(x) = -1.75 - 0.5\sqrt{\max(x_2,0)} + 0.2x_3\), and \(f_{C1}(x) = f_{C0}(x) + 1.15 + 0.5 \cdot 1(x_1 < 0.5) - 0.3\sqrt{\max(x_2,0)}\). The follow-up time is taken to be \(t = 1.5\).

2. Weibull proportional hazards: \(S(t \mid x, a) = \exp(-e^{f_a(x)} \, t^{1/2})\), with \(f_0(x) = x_1\) and \(f_1(x) = x_1 + x_2 - 0.5\). The propensity score is \(e(x) = (1 + f_\beta(x_2; 2, 4))/4\), while censoring is \(C_i \sim \operatorname{Uniform}(0, 3)\). The time horizon is \(t = 1.25\).

10.1.0.3 [43]

Covariates for this setting consist of five continuous variables \(X_1,\ldots,X_5 \stackrel{\textrm{iid}}{\sim}\operatorname{Normal}(0,0.35^2)\) and five binary variables \(X_6,\ldots,X_{10} \stackrel{\textrm{iid}}{\sim}\operatorname{Bernoulli}(0.5)\). The propensity score is \(e(x) = \operatorname{expit}(0.3 - 0.25x_1 - 2.25x_2 - 0.75x_3 - 0.25x_5 - 0.25x_6 - 0.50x_7 - x_9 + 1.25x_{10})\). Potential event times follow Weibull proportional hazards models with shape \(\eta = 2\) and scale parameters \(\kappa_a(x) = d_a \exp\{f_a(x)\}\), where \(d_0 = 1200\) and \(d_1 = 2000\), so that \(T_i(0) = \left[-\log(U_i) / \{1200\exp(f_0(X_i))\}\right]^{1/2}\) and \(T_i(1) = \left[-\log(U_i) / \{2000\exp(f_1(X_i))\}\right]^{1/2}\) with \(U_i \sim \operatorname{Uniform}(0,1)\). Censoring is independent with \(C_i \sim \operatorname{Gam}(1, 0.007)\). The functions \(f_0\) and \(f_1\) are given by

1. DGP 1: \(f_0(x) = 0.2 - 0.5x_1 - 0.8x_3 - 1.8x_5 - 0.9x_6 - 0.1x_7\) and \(f_1(x) = -0.2 + 0.1\operatorname{expit}(x_1) - 0.8\sin(x_3) - 0.1x_5^2 - 0.3x_6 - 0.2x_7\).

2. DGP 2: \(f_0(x) = -0.1 + 0.1x_1^2 - 0.2\sin(x_3) + 0.2\operatorname{expit}(x_5) + 0.2x_6 - 0.3x_7\) and \(f_1(x) = -0.2 + 0.1\operatorname{expit}(x_1) - 0.8\sin(x_3) - 0.1x_5^2 - 0.3x_6 - 0.2x_7\).

3. DGP 3: \(f_0(x) = -0.1 + 0.1x_1^2 - 0.2\sin(x_3) + 0.2\operatorname{expit}(x_5) + 0.2x_6 - 0.3x_7\) and \(f_1(x) = 0.5 - 0.1\operatorname{expit}(x_2) + 0.1\sin(x_3) - 0.1x_4^2 + 0.2x_4 - 0.1x_5^2 + 0.2\operatorname{expit}(x_5) + 0.2x_6 - 0.3x_7\).

4. DGP 4: \(f_0(x) = -0.2 + 0.5\sin(\pi x_1 x_3) + 0.2\operatorname{expit}(x_5) + 0.2x_6 - 0.3x_7\) and \(f_1(x) = 0.5 - 0.1\operatorname{expit}(x_2) + 0.1\sin(x_3) - 0.1x_4^2 + 0.2x_4 - 0.1x_5^2 - 0.3x_6\).

We use a time horizon of \(t = 0.05\) for all settings.

10.2 Results for Conditional Average Treatment Effects↩︎

10.2.0.1 Conditional Average Treatment Effects: With a Cured Subpopulation

Estimating conditional average treatment effects is a much more difficult problem than estimating average treatment effects, and methods differ much more strongly in their performance here. Results for when a cured subpopulation is present are given in Figures 7—9. We did not observe any consistent patterns in terms of overall bias, except that IndivAFT performs very poorly on the Hu settings. For coverage, interval length, and RMSE, BartCure(1) is the most consistent in the sense that it never performs particularly poorly and always outperforms CSF, but no method attains nominal coverage of confidence/credible intervals; an advantage of BartCure(1)even in this context is that, as shown in Section 4.3, it is generally conservative in identifying treatment effect heterogeneity, so in settings where it does detect treatment effect heterogeneity it will generally be underestimating it rather than overestimating it.

10.2.0.2 Conditional Average Treatment Effects: Without a Cured Subpopulation

Results are given in Figures 10—12. We again note no consistent patterns in terms of overall bias, and that BartCure(1) is overall the most robust method and outperforms CSF in terms of RMSE and performance of interval estimands. Interestingly, coverage for all of the effects is much better when there is no cured subpopulation. BartCure(2)remains a more unstable method than BartCure(1). Just as with the average effects, the CSF breaks down entirely on the Henderson settings, and IndivAFT is mostly competitive with BartCure(1) except for the Hu 3 and Hu 4 settings.

11 Additional Analyses of the CALGB 40101 Trial↩︎

11.0.0.1 Proportion of the RMST Effect Attributable to Cure

Plots of the signed and unsigned shares of the effect attributable to cure as defined in Section 2.2 are given in Figure 13. These plots provide further evidence that the majority of the effect is attributable to the cured subpopulation, and generally there is a substantial amount of uncertainty in the stochastic latency effect.

11.0.0.2 Variable Importance

Figure 14 displays the variable importances from the BartCuremodel, quantified using the average number of splitting rules involving each variable as suggested by [9]. The most influential variable in the model was time, which suggests that a proportional hazards variant of the model 2 that takes \(r(t,a,x) = r(a,x)\) would be inadequate, i.e., some variables have a time-varying effect on the hazard. Among baseline clinical covariates, receptor status, tumor size, menopausal status, and age category were the largest contributors to survival. Overall, the variable importance profile suggests that BartCuredistributes predictive importance across several clinically relevant covariates rather than concentrating heavily on a single factor. The prominence of receptor status, tumor size, and age is consistent with established prognostic factors in breast cancer survival [46].

References↩︎