Causal Inference with BART

Introduction

BART is a popular choice for causal inference because it allows one to fit the nuisance functionals required for effect estimation (the relationships among the outcome, the treatment, and the potential confounders) in a flexible manner (Hill 2011). BART has repeatedly been shown to outperform other effect estimation methods in competitions (Dorie et al. 2019). In this vignette, we demonstrate how to use BART implemented in bartisan to estimate treatment effects under the assumption of no unmeasured confounding. This includes both standard BART as well as Bayesian causal forests (BCF), which is an improvement on traditional BART that works by fitting separate BART models for the outcome under control and for the treatment effect itself as a varying coefficient model.

However, it’s important to remember that BART is not a “causal inference method”; it is just a method of estimating certain quantities, which often are interpretable as associations. The assumptions that turn an association into a causal effect are assumptions about the design, and no model supplies them.

In addition to its use here, BART can be used with instrumental variables analysis (McCulloch et al., n.d.) and regression discontinuity (Alcantara et al. 2024).

The example we use here is a dataset used to answer whether right heart catheterization helps or harms critically ill patients (Connors et al. 1996).

library(bartisan)

data(rhc)

set.seed(2026)

Confounding

Catheterization was not randomized. Doctors chose it, and they chose it more often for patients who were already doing badly. Those patients were also more likely to die. We can use cobalt::bal.tab() to see how the patients differ between groups:

bal.tab(rhc ~ age + sex + race + edu + aps + meanbp + resp +
          hema + pafi + paco2 + crea + surv2m + card,
        data = rhc, stats = c("m", "ovl"), disp = "m")
#> Balance Measures
#>               Type  M.0.Un  M.1.Un Diff.Un OVL.Un
#> age        Contin.  61.816  60.756  -0.065  0.101
#> sex_male    Binary   0.524   0.591   0.067  0.067
#> race_white  Binary   0.793   0.768  -0.024  0.024
#> race_black  Binary   0.157   0.168   0.011  0.011
#> race_other  Binary   0.050   0.064   0.013  0.013
#> edu        Contin.  11.568  11.767   0.061  0.048
#> aps        Contin.  51.369  62.099   0.530  0.223
#> meanbp     Contin.  85.302  66.871  -0.505  0.231
#> resp       Contin.  28.841  27.465  -0.097  0.055
#> hema       Contin.  32.545  30.237  -0.290  0.159
#> pafi       Contin. 238.232 182.995  -0.504  0.208
#> paco2      Contin.  40.045  36.871  -0.262  0.096
#> crea       Contin.   1.905   2.506   0.292  0.179
#> surv2m     Contin.   0.604   0.559  -0.229  0.105
#> card_yes    Binary   0.299   0.405   0.106  0.106
#> 
#> Sample sizes
#>     Control Treated
#> All     935     565

The M.0.Un column indicates the mean for each variable in the control group and the M.1.Un column indicates the mean for each variable in the treated group. Diff.Un and OVL.Un are measures of the distributional difference between the groups for each covariate; values far from 0 indicate imbalance due to differential selection into treatment. In particular, we can see that patients with higher values of aps and crea and lower values of meanbp, hema, pafi, paco2, and surv2m are overrepresented among treated units.

Given that sicker patients are both more likely to die and more likely to receive RHC, it would not be unexpected to see that patients who receive RHC are more likely to die, and we do see that:

with(rhc, tapply(death, rhc, mean))
#>      0      1 
#> 0.6193 0.7115

However, that doesn’t mean RHC causes death; to disentangle the effects of RHC from the confounding effects of patients’ characteristics, we need to adjust for these characteristics.

Assumptions for Causal Inference

Several assumptions are required to interpret an adjusted effect estimate as causal, none of which the fit can check:

1. No unmeasured confounding. Every common cause of catheterization and death is in the model. Here that is the crux: the covariates include a physiological profile and the study’s own prognostic score, which is a serious attempt, but a doctor’s judgment at the bedside may not be fully captured by these thirteen variables.

2. Positivity. Every kind of patient could have received the procedure or not. This one is partly checkable and is checked below.

3. Consistency. “Catheterization” names a single well defined intervention.

4. No interference. One patient’s treatment does not affect another’s outcome.

The first is the one that usually fails and the one to be explicit about. State it as an assumption in what you write, rather than letting the interval imply it has been handled.

ps_fit <- bartisan(
  rhc ~ age + sex + race + edu + aps + meanbp + resp + hema + pafi +
    paco2 + crea + surv2m + card,
  data = rhc, family = binomial()
)

prop_score <- fitted(ps_fit)

bal.plot(rhc ~ prop_score, data = rhc, type = "hist", mirror = TRUE)

The Estimator: G-Computation

To estimate the average treatment effect \(\tau_{\text{ATE}} = E[Y(1)] - E[Y(0)]\), which is a function of the unobserved potential outcomes \(Y(1)\) and \(Y(0)\), we can use the assumptions above to express it as a function of observed quantities:

\[E[Y(1)] - E[Y(0)] = E \left[ E[Y | X, A = 1] \right] - E \left[ E[Y | X, A = 0] \right]\] To estimate \(E \left[ E \left[ Y | X, A = a \right] \right] = \theta_a\), we use g-computation (Snowden et al. 2011):

\[\hat{\theta}_a = \frac{1}{n}\sum_{i=1}^n {\mu(a, x_i)}\]

where \(\mu(a, x_i)\) is the predicted value from a regression of \(Y\) on \(A\) and \(X\) for a unit with covariate profile \(X=x_i\) and treatment \(A\) set to \(a\). This estimator is sometimes also known as the “regression estimator” or the “plug-in estimator”.

Here, we use traditional BART and BCF to model \(\mu(A, X)\). The rest of the analysis comes from the definitions above.

The Outcome Model

We can fit the outcome model as a binary BART regression of the outcome on the treatment and covariates, as in Hill (2011), with the propensity score estimated above added as one more covariate:

rhc$prop_score <- prop_score

fit <- bartisan(death ~ rhc + age + sex + race + edu + aps + meanbp + resp +
                  hema + pafi + paco2 + crea + surv2m + card + prop_score,
                data = rhc, family = binomial())

We include the treatment, the covariates, and the propensity score in the formula’s right-hand side, and set family = binomial() to model the binary outcome with logistic regression. Including the propensity score is recommended by Carnegie (2019), though Souto and Louzada (2024) report that it often makes little difference to the estimate. Linero (2024) gives the Bayesian reason for including it: when the outcome model is flexible and its prior is independent of the model for the treatment, the prior it implies for the amount of confounding bias is concentrated near zero, strongly enough that the data may not overcome it, and adjusting for the propensity score is one way to relax that. Because the score is a function of the covariates, adding it does not change what the model conditions on, and the g-computation below holds it at each patient’s own value.

The splitting prior is left at its default, which weights the predictors equally and cannot drop any of them. Setting sparsity = TRUE would place a variable-selection prior on them instead, which can drop the treatment from the forest in some draws and put a point mass at exactly zero in the posterior of its effect, so it is not the setting for a model whose purpose is a contrast on one predictor; ?bartisan_control explains the trade-off and vignette("effects") covers what the point mass does to an estimate. The propensity score model is a prediction problem with no contrast read off it, so the sparsity prior can be turned on there when there are many covariates.

The fits in this vignette use the default single chain for brevity. Normally, we would run several chains, possibly change the number of burn-in and retained draws, and examine convergence diagnostics to make sure the model was fit correctly; see vignette("diagnostics") for how to set these fit controls and check them.

Potential Outcomes

The quantities underlying the effect estimate are the two average potential outcomes: the proportion who would die if every patient were catheterized, and the proportion who would die if none were. estimate_effect() computes them on the way to the effect and prints them beneath it. On a fit from bartisan(), the treatment has to be named, since nothing in the formula marks one predictor as the treatment:

ate <- estimate_effect(fit, treat = "rhc")

ate
#> Average treatment effect (difference)
#> 
#> Treatment: `rhc`
#> Averaged over 1500 units
#> 
#>     contrast estimate  lower upper    n
#>  Y[1] - Y[0]    0.063 0.0145 0.113 1500
#> 
#> Average potential outcomes
#> 
#>  quantity estimate lower upper
#>      Y[0]    0.631 0.600 0.659
#>      Y[1]    0.694 0.655 0.730
#> 
#> ℹ estimate is the posterior mean; lower and upper bound the 95% equal-tailed
#>   credible interval.
#> ℹ Y[a] is the average response with `rhc` set to "a".

Under the assumptions above this is the average treatment effect: catheterization raises the probability of death by about 6 percentage points, with an interval running from roughly 1.4 to 11.3.

The two rows below the contrast are the estimates of \(E[Y(0)]\) and \(E[Y(1)]\), each averaged over the observed covariate distribution. Reporting both is often more informative than reporting their difference alone, because a difference of a few percentage points means something different against a baseline of 63% than it would against 5%. Setting potential_outcomes = FALSE in the print() call drops them where the difference is all that is wanted.

estimate_effect(fit, treat = "rhc", comparison = "lnor")
#> Average treatment effect (log odds ratio)
#> 
#> Treatment: `rhc`
#> Averaged over 1500 units
#> 
#>                contrast estimate  lower upper    n
#>  log(O(Y[1]) / O(Y[0]))    0.284 0.0633 0.513 1500
#> 
#> Average potential outcomes
#> 
#>  quantity estimate lower upper
#>      Y[0]    0.631 0.600 0.659
#>      Y[1]    0.694 0.655 0.730
#> 
#> ℹ estimate is the posterior mean; lower and upper bound the 95% equal-tailed
#>   credible interval.
#> ℹ Y[a] is the average response with `rhc` set to "a", and O(y) is the odds
#>   `y/(1-y)`.

options(marginaleffects_posterior_center = mean)

avg_comparisons(fit, variables = "rhc")
#> 
#>  Estimate  2.5 % 97.5 %
#>     0.063 0.0145  0.113
#> 
#> Term: rhc
#> Type: response
#> Comparison: 1 - 0

Using Bayesian Causal Forests

Traditional BART shrinks the fitted function toward a constant; that shrinkage applies to everything the forest fits, including the part of the outcome that depends on the exposure. When the exposure is strongly predicted by the covariates, the forest can explain the outcome using the covariates alone, leaving little for the exposure to explain, and shrink the estimated effect toward zero. The interval shrinks with it, so the result is a confident estimate biased toward no effect.

Hahn et al. (2020) identified this mechanism; their remedy is to give the treatment effect its own forest with its own prior, so that shrinking the confounding part does not shrink the effect. This is the Bayesian causal forest (BCF) model, a special case of the varying coefficients BART model; bcf() fits it.

One specifies the control function in the model formula and identifies the treatment in the treat argument. bcf() then fits a varying coefficient BART model, the BCF. By default, bcf() estimates a propensity score using a logistic BART model and includes that as a covariate in the control function, as recommended by Hahn et al. (2020) and for the reason given above for the plain BART model.

fit_bcf <- bcf(death ~ age + sex + race + edu + aps + meanbp + resp + hema +
                 pafi + paco2 + crea + surv2m + card,
               treat = ~ rhc, data = rhc,
               family = binomial())

fit_bcf
#> Generalized BART
#> 
#> Call:
#> bcf(formula = death ~ age + sex + race + edu + aps + meanbp + 
#>     resp + hema + pafi + paco2 + crea + surv2m + card, treat = ~rhc, 
#>     data = rhc, family = binomial())
#> 
#> Family: "binomial" with the "logit" link
#> Observations: 1500
#> Structure: 2 forests of 50 and 25 trees, soft decision rules
#> Draws: 800 kept after 200 warmup
#> 
#> Posterior means: b.rhc.0 = -0.13, b.rhc.1 = -0.209
#> 
#> Treatment: "rhc"
#> Effect moderators: "age", "sex", "race", "edu", "aps", "meanbp", "resp", "hema", "pafi", "paco2", "crea", "surv2m", and "card"
#> ℹ `estimate_effect()` reports the treatment effect, with the average potential
#>   outcomes beside it; `plot()` draws the conditional ones.

Using bcf() rather than writing the varying coefficient BART model by hand sets four settings to improve effect estimation:

1. the effect gets a forest of its own
2. the estimated propensity score goes into the control function and not into the effect forest
3. the effect forest gets fewer trees because effect heterogeneity is usually simpler than a prognostic surface
4. a binary treatment’s coding is drawn rather than fixed, so the answer does not depend on which arm was written as 1.

The splitting prior stays at its default of off for both forests. Because the treatment is the coefficient rather than a predictor the forest splits on, nothing can drop it, so the sparsity prior can be turned on for the effect forest alone with sparsity = c(FALSE, TRUE), which asks that forest to select among the moderators; ?bcf discusses when that is worth doing.

Because the treatment is named in the call, estimate_effect() needs nothing else:

estimate_effect(fit_bcf)
#> Average treatment effect (difference)
#> 
#> Treatment: `rhc`
#> Averaged over 1500 units
#> 
#>     contrast estimate   lower upper    n
#>  Y[1] - Y[0]   0.0563 0.00241 0.108 1500
#> 
#> Average potential outcomes
#> 
#>  quantity estimate lower upper
#>      Y[0]    0.633 0.605 0.662
#>      Y[1]    0.689 0.650 0.728
#> 
#> ℹ estimate is the posterior mean; lower and upper bound the 95% equal-tailed
#>   credible interval.
#> ℹ Y[a] is the average response with `rhc` set to "a".

The two potential outcomes are printed beneath the contrast, since a difference of a few percentage points means one thing against a baseline of 63% and another against 5%; setting potential_outcomes = FALSE in print() suppresses them. The contrast’s label names the quantity rather than leaving it to the heading, which matters once a ratio is asked for: Y[1] - Y[0] is a difference of average responses where log(O(Y[1]) / O(Y[0])) is a log odds ratio, and the note beneath the table says what Y[a] is, adding what O(y) is once an odds ratio is asked for.

The output of summary() on a bcf() fit is the same summary of the forests it is on any other fit, and says at the end where the effect is reported. That way the same call means the same thing whichever way the model was written.

cate <- estimate_effect(fit_bcf, estimand = "CATE", comparison = "or")

quantile(cate$estimate, probs = c(0, .25, .5, .75, 1))
#>    0%   25%   50%   75%  100% 
#> 1.141 1.324 1.420 1.523 1.839

plot(fit_bcf)

A Second Example: A Continuous Outcome, and the ATT

The catheterization question is about a whole population, so the ATE is the estimand. Many questions are not. When a program is offered to a particular group and the question is whether it helped them, the ATT is the estimand to report, and the outcome is often continuous rather than binary.

lalonde, from cobalt, is the standard example: a job training program, with earnings in 1978 as the outcome.

data("lalonde", package = "cobalt")

with(lalonde, tapply(re78, treat, mean))
#>    0    1 
#> 6984 6349

Participants earned less than non-participants. Taken at face value, the program looks harmful, but that doesn’t mean it actually is: participants were selected for being out of work, so they would have earned less anyway. This is confounding in the opposite direction to the catheterization example, and it is why the raw comparison is worth showing before the adjusted one.

Below, we fit the BCF model, this time setting family = tweedie(), which is useful for continuous outcomes with a point mass at 0, as we have here. See vignette("families") for other families that could be used here.

fit_earn_bcf <- bcf(
  re78 ~ age + educ + race + married + nodegree + re74 + re75,
  treat = ~ treat,
  data = lalonde, family = tweedie()
)

estimate_effect(fit_earn_bcf, estimand = "ATT")
#> Average treatment effect on the treated (difference)
#> 
#> Treatment: `treat`
#> Averaged over the 185 units in group "1"
#> 
#>     contrast estimate lower upper   n
#>  Y[1] - Y[0]     1120  -363  2540 185
#> 
#> Average potential outcomes
#> 
#>  quantity estimate lower upper
#>      Y[0]     5230  4240  6330
#>      Y[1]     6350  5410  7410
#> 
#> ℹ estimate is the posterior mean; lower and upper bound the 95% equal-tailed
#>   credible interval.
#> ℹ Y[a] is the average response with `treat` set to "a".

Setting estimand = "ATT" averages over the treated rather than over everyone.¹

The potential outcomes earn their place here: a difference of about a thousand dollars means something different against a baseline of about five thousand than it would against five hundred.

Because the link is the identity, the ATT is exactly the average of the conditional effects among the treated, which is worth checking once to see that the two agree:

att <- estimate_effect(fit_earn_bcf, estimand = "ATT")

cate_att <- estimate_effect(fit_earn_bcf, estimand = "CATE",
                            newdata = subset(lalonde, treat == 1))

c(ATT = att$estimate, mean_CATE = mean(cate_att$estimate))
#>       ATT mean_CATE 
#>      1125      1125

And the conditional effects, drawn:

plot(cate_att)

Most of the conditional effects are positive, as the ATT is, and none of their 185 intervals exclude zero. That is the usual picture: the effect at one unit’s covariates is estimated from far less information than an average over all of them, so its interval is much wider.

Interpreting the Credible Interval

The posterior interval is a credible interval for the estimand under the model and under the identification assumptions. It covers uncertainty about the outcome surface at the covariate values in the sample, which comes from having a finite sample and from not knowing the surface’s shape. As described in the section on conditional, sample, and population averages, it treats those covariates as fixed, so it does not include the additional uncertainty an interval for the population average would carry.

For a bcf() fit, the interval also takes the propensity score as known. The score is the posterior mean of a separate model for the treatment, estimated once before the outcome model is fit, so uncertainty in it does not reach the interval. Estimating the two models in stages like this is the usual practice, and it keeps the outcome from informing the propensity score, which a joint model of the two would allow; Li et al. (2023) discuss both approaches and the alternatives between them.

The interval does not cover uncertainty about whether the assumptions hold.

An unmeasured confounder does not widen the interval. It moves the estimate and leaves the interval where it was. This is why a sensitivity analysis, asking how strong a confounder would have to be to explain the result away, is a more informative addition than any refinement of the model.

Further Reading

?estimate_effect is the reference for the estimands used here, including the subgroup form and the contrast types. vignette("effects") covers what marginaleffects adds beyond them, which is a much wider class of quantities and the route to take when the question is not a treatment contrast. vignette("diagnostics") covers checking the fit, which should happen before any of this is interpreted.