Chapter 15
Bayesian nonparametrics for causal questions

Lectures 35–37 · priors on distributions and on functions, the Dirichlet process and the Bayesian bootstrap, Bayesian causal forests, and why nonparametric does not mean assumption-free

Every model so far has had a fixed, finite set of parameters: the mean and standard deviation of a normal, the coefficients of a regression, the means and weights of a mixture with $K$ components. When the family is right, that is efficient; when it is wrong, the posterior can be confident and mistaken. Chapter 13's job-training analysis needed a hurdle lognormal for earnings, and the answer for individuals depended on how well it fitted the tails. Ambulance response times in New York have a long right tail and several bumps. Treatment effects in a national education experiment may vary with a dozen characteristics of students and schools in ways nobody would think to write down.

Bayesian nonparametric models put priors on objects that are not described by a fixed number of parameters: on whole probability distributions, and on whole functions. The name is misleading, because such models have more parameters, not fewer, infinitely many, of which the data determine how many matter. This chapter gathers threads that have appeared through the course (the Gaussian processes of Chapter 8, the regression trees of Chapters 10 and 13, the question in Chapter 11 of how many clusters to fit) into one framework, and applies it to causal questions, where flexible models matter most and are most easily misunderstood.

• Lecture 35 is about priors on distributions. Why a parametric family can mislead; the Dirichlet process, its stick-breaking and restaurant descriptions; the Bayesian bootstrap; Dirichlet-process mixtures for the density of ambulance response times; and nonparametric Bayesian estimates of causal effects.
• Lecture 36 is about priors on functions. Gaussian processes and Bayesian additive regression trees as priors over response surfaces; regularisation-induced confounding; and Bayesian causal forests applied to the National Study of Learning Mindsets.
• Lecture 37 is about the limits. A causal inference competition with known answers; how often flexible models' intervals cover the truth; overlap and extrapolation; and how to report nonparametric analyses.

Lecture 35 — Priors on distributions

15.1 When the family is wrong

Chapter 6 analysed the response times of New York City's ambulances to the most serious medical emergencies, on the log scale, with a multilevel normal model. Before trusting any normal model for these data, look at their distribution. The data are the 2024 calls with the highest severity level from NYC Open Data, cached locally after the first download.

library(tidyverse)
library(dirichletprocess)
library(stochtree)
library(dbarts)
library(jsonlite)
options(mc.cores = parallel::detectCores())

harbour <- c(ink = "#10191d", muted = "#4a5f68", grid = "#d5e3e6",
             sky = "#1d5187", amber = "#805907", rose = "#ad3b28",
             moss = "#1d6b54", lilac = "#5b3a94")
theme_harbour <- function(base_size = 11) {
  theme_minimal(base_size = base_size) +
    theme(text = element_text(colour = harbour[["ink"]]),
          axis.text = element_text(colour = harbour[["muted"]]),
          panel.grid.major = element_line(colour = harbour[["grid"]]),
          panel.grid.minor = element_blank(),
          strip.text = element_text(colour = harbour[["ink"]], hjust = 0),
          legend.position = "top")
}
theme_set(theme_harbour())

f <- "data/ems/ems_2024_severity1.csv"
if (!file.exists(f)) {
  dir.create("data/ems", showWarnings = FALSE)
  url <- paste0("https://data.cityofnewyork.us/resource/76xm-jjuj.json?",
    "$select=incident_datetime,initial_call_type,borough,zipcode,incident_response_seconds_qy,valid_incident_rspns_time_indc",
    "&$where=initial_severity_level_code='1' AND incident_datetime between '2024-01-01T00:00:00' and '2024-12-31T23:59:59'",
    "&$limit=50000")
  write_csv(as_tibble(fromJSON(URLencode(url))), f) }
ems <- read_csv(f, show_col_types = FALSE) |>
  filter(valid_incident_rspns_time_indc == "Y", incident_response_seconds_qy > 0) |>
  mutate(minutes = incident_response_seconds_qy / 60, borough = str_to_title(str_remove(borough, "^RICHMOND / ")))
c(calls = nrow(ems), quantile(ems$minutes, c(0.01, 0.1, 0.5, 0.9, 0.99))) |> round(2)
set.seed(15); y <- sample(log(ems$minutes), 1500)            # a random 1,500 calls for the models below
c(skewness_of_log = mean((y - mean(y))^3) / sd(y)^3, excess_kurtosis_of_log = mean((y - mean(y))^4) / sd(y)^4 - 3) |> round(2)

   calls       1%      10%      50%      90%      99% 
24404.00     1.22     3.08     5.57     9.28    16.32 
       skewness_of_log excess_kurtosis_of_log 
                 -0.58                   2.89 

The 24,404 calls in 2024 with valid response times had a median response of 5.6 minutes; 10% took less than 3.1 minutes and 10% more than 9.3, and 1% took more than 16.3. On the log scale, where Chapter 6's model assumed normality, the distribution is skewed to the left (skewness −0.58) and has heavier tails than a normal (excess kurtosis 2.89): there are more very fast responses, and more very slow ones, than a lognormal allows. A normal model of log times would describe the middle well and the tails, where equity questions live, badly.

Which family should we use? We could try a few, lognormal, gamma, a mixture of two lognormals, and compare them, as in earlier chapters. The nonparametric answer is different: put a prior directly on the distribution, wide enough to contain all of these and much else, and let the data decide its shape.

15.2 The Dirichlet process

A prior on distributions must produce a random probability distribution $G$ each time it is sampled. The most important such prior, introduced by Thomas Ferguson (1973), is defined by what it says about every way of cutting the space into pieces.

Two consequences follow from properties of the Dirichlet distribution. The expected value of $G(A)$ is $G_0(A)$, so $G_0$ is the prior guess; and its variance is $G_0(A)(1-G_0(A))/(\alpha +1)$, so $\alpha$ says how tightly $G$ concentrates around the guess, like a prior sample size. The definition is abstract; three constructive descriptions make it concrete.

• Stick-breaking (Sethuraman, 1994). Take a stick of length 1. Break off a fraction $V_1\sim \text{beta}(1,\alpha )$; from what remains break off a fraction $V_2\sim \text{beta}(1,\alpha )$; and so on forever. The pieces have lengths $w_1=V_1$, $w_2=V_2(1-V_1)$, ..., and they sum to 1. Place each piece at a location ${\theta }_h$ drawn from $G_0$. Then $G={\sum }_h{w}_h{\delta }_{{\theta }_h}$, a discrete distribution with infinitely many atoms, most of them tiny. Small $\alpha$ makes the first pieces large: $G$ is dominated by a few atoms. Large $\alpha$ spreads the mass over many.
• The Pólya urn (Blackwell and MacQueen, 1973). Draw ${\theta }_1,{\theta }_2,\dots$ from $G$ with $G$ integrated out: ${\theta }_1\sim G_0$, and each later draw either repeats one of the earlier values, chosen in proportion to how often each has appeared, or, with probability $\alpha /(\alpha +n)$ after $n$ draws, is a new value from $G_0$.
• The Chinese restaurant process. The same urn, told as customers arriving at a restaurant with infinitely many tables: each new customer sits at an occupied table with probability proportional to the number already there, or at a new table with probability proportional to $\alpha$ (the widget of Section 11.8 simulates it). The rich get richer, and the number of occupied tables after $n$ customers grows like $\alpha \log (1+n/\alpha )$.

The posterior is as simple as the prior. Given observations $x_1,\dots ,x_n$ from $G$,

a Dirichlet process whose base measure mixes the prior guess with the empirical distribution of the data, weighted by $\alpha$ and $n$. It is the beta–binomial update of Chapter 1, carried out on every set at once.

15.3 The Bayesian bootstrap

Let the concentration go to zero. The posterior then puts all its weight on the observed values, and the weights it gives them are $\text{Dirichlet}(1,\dots ,1)$: the posterior of the unknown distribution is a random reweighting of the data. This is Donald Rubin's (1981) Bayesian bootstrap. To get the posterior of any summary of the distribution (a mean, a median, a 90th percentile, a share above a threshold), draw weights $w\sim \text{Dirichlet}(1,\dots ,1)$, compute the weighted summary, and repeat. It is Efron's bootstrap with smooth weights instead of whole-number resampling counts, and it shows that the bootstrap, usually justified by frequentist arguments, also has a Bayesian reading: it is the posterior under a model that assumes nothing about the shape of the distribution, only that the observed values are the possible ones.

bb <- function(x, stat, draws = 4000) {                       # Bayesian bootstrap posterior of stat(x, w)
  map_dbl(seq_len(draws), function(i) { w <- rexp(length(x)); stat(x, w / sum(w)) }) }
wq <- function(x, w, p) { o <- order(x); x[o][which(cumsum(w[o]) >= p)[1]] }
set.seed(15)
ems |> group_by(borough) |>
  group_modify(~ { m <- .x$minutes
    p90 <- bb(m, function(x, w) wq(x, w, 0.9), 1000); over10 <- bb(m, function(x, w) sum(w * (x > 10)), 1000)
    tibble(calls = length(m), p90_mean = mean(p90), p90_lo = quantile(p90, 0.05), p90_hi = quantile(p90, 0.95),
           over10_mean = mean(over10), over10_lo = quantile(over10, 0.05), over10_hi = quantile(over10, 0.95)) }) |>
  mutate(across(where(is.double), ~ round(.x, 3)))

# A tibble: 5 x 8
# Groups:   borough [5]
  borough       calls p90_mean p90_lo p90_hi over10_mean over10_lo over10_hi
  <chr>         <int>    <dbl>  <dbl>  <dbl>       <dbl>     <dbl>     <dbl>
1 Bronx          5483     9.41   9.25   9.58       0.078     0.072     0.083
2 Brooklyn       6539     8.90   8.78   9.02       0.058     0.053     0.063
3 Manhattan      5597     8.84   8.68   9          0.059     0.054     0.064
4 Queens         5528    10.1    9.98  10.3        0.106     0.099     0.112
5 Staten Island  1257     9.01   8.72   9.42       0.059     0.048     0.071

Queens stands out: its 90th percentile response time is 10.1 minutes (90% interval 10.0 to 10.3) and 10.6% of its calls took more than ten minutes (9.9% to 11.2%), against about 5.8% in Brooklyn and Manhattan and 7.8% in the Bronx. Staten Island, with only 1,257 calls, has the widest intervals. No distributional assumption was needed. One assumption remains: the Bayesian bootstrap treats calls as exchangeable draws, while calls on the same day share traffic and staffing, so these intervals are somewhat too narrow; a version that reweights days rather than calls would respect that dependence.

15.4 Mixtures with no fixed number of components

A distribution drawn from a Dirichlet process is discrete, which is fine for the Bayesian bootstrap but useless as a model for a continuous density. The remedy is to use $G$ as a mixing distribution: give each observation its own parameter drawn from $G$, and its value from a smooth kernel around that parameter,

Because $G$ is discrete, observations share parameter values, and the model is a mixture of normals with a random and unbounded number of components: a Dirichlet-process mixture. It answers the question of Chapter 11, how many components, by not fixing the number at all. We fit it to the log response times of 1,500 calls with the dirichletprocess package, which uses a Gibbs sampler built on the restaurant process and puts a gamma prior on $\alpha$, and compare the fitted density with a lognormal.

dp_path <- "fits/ch15_dpm.rds"
z <- (y - mean(y)) / sd(y)                                   # the package's default base measure suits standardised data
if (file.exists(dp_path)) dp <- readRDS(dp_path) else {
  set.seed(15); dp <- Fit(DirichletProcessGaussian(z), 1500, progressBar = FALSE); saveRDS(dp, dp_path) }
keep <- 501:1500
clusters <- map_int(dp$labelsChain[keep], ~ length(unique(.x)))
big <- map_int(dp$labelsChain[keep], ~ sum(table(.x) >= 0.05 * length(z)))  # components holding at least 5% of calls
c(alpha_mean = mean(dp$alphaChain[keep]), components_median = median(clusters), components_q05 = unname(quantile(clusters, 0.05)),
  components_q95 = unname(quantile(clusters, 0.95)), large_components_median = median(big)) |> round(2)

             alpha_mean       components_median          components_q05 
                   5.85                   54.00                   33.00 
         components_q95 large_components_median 
                  77.00                    5.00 

xg <- seq(min(z) - 0.5, max(z) + 0.5, length.out = 300)
pf <- PosteriorFrame(dp, xg, ndraws = 200, ci_size = 0.1)
dens <- tibble(log_min = mean(y) + sd(y) * xg, dp_mean = pf$Mean / sd(y), lo = pf$X5. / sd(y), hi = pf$X95. / sd(y),
               lognormal = dnorm(mean(y) + sd(y) * xg, mean(y), sd(y)))
ggplot() +
  geom_histogram(aes(x = y, y = after_stat(density)), bins = 60, fill = harbour[["grid"]], colour = NA) +
  geom_ribbon(data = dens, aes(log_min, ymin = lo, ymax = hi), fill = harbour[["rose"]], alpha = 0.25) +
  geom_line(data = dens, aes(log_min, dp_mean), colour = harbour[["rose"]], linewidth = 0.9) +
  geom_line(data = dens, aes(log_min, lognormal), colour = harbour[["sky"]], linewidth = 0.8, linetype = 2) +
  scale_x_continuous(breaks = log(c(1, 2, 5, 10, 20)), labels = c(1, 2, 5, 10, 20)) +
  labs(x = "response time (minutes, log scale)", y = "density")

The concentration parameter's posterior mean is 5.9, and the mixture uses a median of 54 components (90% interval 33 to 77), although only about five of them hold more than 5% of the calls each; the rest are small components that shape the tails. The fitted density follows the histogram closely, including a shoulder of fast responses that the lognormal misses and a heavier right tail. The number of components, though, should not be read as the number of kinds of call: as Miller and Harrison (2013) showed, a Dirichlet-process mixture keeps adding small components as the sample grows, even when the truth has only a few.

15.5 Nonparametric estimates of causal effects

The Bayesian bootstrap gives the simplest possible Bayesian analysis of a randomised experiment: a posterior for the difference in mean outcomes that assumes nothing about the shape of the outcome distributions. For the job-training experiment of Section 13.3, reweight the trainees and the controls separately and take the difference of weighted means.

gh <- "https://raw.githubusercontent.com/cran/"
tmp <- tempfile(); download.file(paste0(gh, "causaldata/master/data/nsw_mixtape.rda"), tmp, quiet = TRUE, mode = "wb"); load(tmp)
nsw <- nsw_mixtape |> transmute(treat, re78 = re78 / 1000)
set.seed(15)
wmean <- function(x) { w <- rexp(length(x)); sum(w * x) / sum(w) }
ate_bb <- map_dbl(1:4000, ~ wmean(nsw$re78[nsw$treat == 1]) - wmean(nsw$re78[nsw$treat == 0]))
c(mean = mean(ate_bb), quantile(ate_bb, c(0.05, 0.95)), prob_positive = mean(ate_bb > 0)) |> round(3)

         mean            5%           95% prob_positive 
        1.803         0.744         2.884         0.999 

The Bayesian bootstrap gives a posterior mean effect of $1,800, with a 90% interval from $740 to $2,880 and a probability of 0.999 that the effect is positive. That is centred near the hurdle-lognormal model's estimate of Section 13.3 ($1,690) but much narrower than its 90% interval ($90 to $3,380). The difference is instructive. The lognormal model allows earnings far larger than any observed, and the uncertainty about its tail feeds into the uncertainty about the mean; the Bayesian bootstrap assumes that only the observed values can occur, so it cannot imagine a trainee richer than the richest one seen. For heavy-tailed outcomes the nonparametric answer is not automatically the more cautious one (Exercise 3).

Observational studies need more: the effect must be computed by g-computation (Section 13.9), averaging an outcome model's predictions with and without treatment over the distribution of the covariates. A fully nonparametric Bayesian analysis models both parts flexibly, the outcome surface with BART or a Gaussian process (Lecture 36) and the covariate distribution with a Dirichlet process or its limit, the Bayesian bootstrap. Using Bayesian-bootstrap weights instead of the fixed sample average adds the uncertainty about the population the covariates come from, so the posterior describes the super-population effect, not just the effect in the sample.

The Bayesian bootstrap has a second derivation that leads somewhere new. A parameter such as a mean is a property of the whole population. If we could see the population, we would know the parameter; our uncertainty about it is uncertainty about the observations we have not seen. So, Edwin Fong, Chris Holmes and Stephen Walker (2023) argued, the posterior can be computed by imputing the missing population: draw $y_{n+1}$ from a predictive distribution given $y_{1:n}$, then $y_{n+2}$ given $y_{1:n+1}$, and so on, and compute the parameter from the completed sequence. Repeating the imputation gives a distribution for the parameter, the martingale posterior. It is a proper posterior if the predictive distributions form a martingale: the predictive for the next-but-one observation, averaged over the next, is the predictive for the next. Then the imputed sequence settles down (its empirical distribution converges), and the parameter has a well-defined limit.

The simplest predictive rule is the Pólya urn: the next value is a copy of one of the values seen so far, chosen at random, including the imputed ones. It is the Dirichlet-process predictive of this lecture with $\alpha \to 0$. Its martingale posterior is exactly the Bayesian bootstrap, as the simulation shows for the NSW earnings.

gh <- "https://raw.githubusercontent.com/cran/"
tmp <- tempfile(); download.file(paste0(gh, "causaldata/master/data/nsw_mixtape.rda"), tmp, quiet = TRUE, mode = "wb"); load(tmp)
nsw <- nsw_mixtape |> transmute(treat, re78 = re78 / 1000)
polya_impute <- function(x, N) { n <- length(x); v <- c(x, numeric(N))   # Polya urn: copy a random earlier value
  for (i in 1:N) v[n + i] <- v[sample.int(n + i - 1, 1)]; v }
t1 <- nsw$re78[nsw$treat == 1]; t0 <- nsw$re78[nsw$treat == 0]
set.seed(15)
mp <- map_dbl(1:2000, ~ mean(polya_impute(t1, 5000)) - mean(polya_impute(t0, 5000)))   # impute 5,000 more per arm
wmean <- function(x) { w <- rexp(length(x)); sum(w * x) / sum(w) }
bb <- map_dbl(1:2000, ~ wmean(t1) - wmean(t0))
rbind(martingale_posterior = c(mean = mean(mp), sd = sd(mp), quantile(mp, c(0.05, 0.95))),
      bayesian_bootstrap   = c(mean = mean(bb), sd = sd(bb), quantile(bb, c(0.05, 0.95)))) |> round(3)
set.seed(2); v <- polya_impute(t1, 20000)
round(c(cumsum(v) / seq_along(v))[length(t1) + c(0, 10, 100, 1000, 5000, 20000)], 3)   # running mean settles down

                      mean    sd    5%   95%
martingale_posterior 1.787 0.677 0.709 2.935
bayesian_bootstrap   1.788 0.685 0.714 2.981
[1] 6.349 6.304 6.132 6.275 6.231 6.072

The two posteriors for the effect of the training programme agree to within simulation error: mean $1,790, standard deviation about $680, 90% interval roughly $710 to $2,950. The running mean of one imputed population of treated earnings moves from $6,350 (the observed mean) to $6,070 after 20,000 imputations, changing less and less: the martingale has converged, and its limit is one draw from the posterior.

Why bother, if the answer is the Bayesian bootstrap? Because the predictive rule can be anything that forms a martingale, and it need not come from a likelihood and a prior. Fong, Holmes and Walker use a smooth rule, a sequence of bivariate copula updates, which gives a posterior for densities and regression functions with no prior specified at all, and which can be computed in parallel quickly. The approach puts the emphasis where Section 13.3 put it, on predicting the unseen, and it offers a principled answer to the question this lecture keeps returning to: how to be Bayesian when you do not believe any parametric model.

Lecture 36 — Priors on functions

15.6 Priors on functions

Causal effects are differences between values of a function. Write $\mu (x,t)=𝔼[Y\mid X=x,T=t]$ for the response surface; under ignorability (Section 13.2) the conditional average treatment effect is $\tau (x)=\mu (x,1)-\mu (x,0)$, and average effects are averages of $\tau (x)$ over the covariates. A linear model forces $\mu$ into a narrow family; a nonparametric prior lets the data shape it. Two such priors have already appeared in this course.

• A Gaussian process (Chapter 8) makes any finite set of function values jointly normal, with a covariance kernel saying how quickly the function can change. Its posterior is smooth, and far from the data it reverts to the prior mean with growing uncertainty.
• Bayesian additive regression trees (BART; Chipman, George and McCulloch, 2010; Chapter 10) write the function as a sum of many small regression trees, each kept shallow and its leaf values shrunk toward zero by the prior, so that each tree explains a little. BART adapts to interactions and thresholds without being told about them, and it needs almost no tuning; that is why Jennifer Hill (2011) proposed it as a default outcome model for causal inference, and why it has done well in benchmarks ever since.

Both are consistent: with enough data spread across the covariate space, their posteriors concentrate on the true function. But "spread across the covariate space" is the crux. The figure shows both priors fitted to simulated data in which treated units are missing from part of the range of the confounder: the situation that Section 13.9 called a failure of overlap.

set.seed(15)
x <- runif(150, -2, 1); f <- function(x) 1 + x + 0.5 * x^2; yv <- f(x) + rnorm(150, 0, 0.5)
xg <- seq(-2, 3, length.out = 120)
bfit <- bart(x.train = matrix(x), y.train = yv, x.test = matrix(xg), ndpost = 1000, nskip = 500, verbose = FALSE)
bart_band <- tibble(x = xg, mean = colMeans(bfit$yhat.test), lo = apply(bfit$yhat.test, 2, quantile, 0.05),
                    hi = apply(bfit$yhat.test, 2, quantile, 0.95), model = "BART")
gp_fit <- function(x, y, xs) {                                 # GP with squared-exponential kernel; hyperparameters by marginal likelihood
  nll <- function(p) { l <- exp(p[1]); s <- exp(p[2]); n <- exp(p[3])
    K <- s^2 * exp(-outer(x, x, "-")^2 / (2 * l^2)) + diag(n^2, length(x)); L <- chol(K)
    a <- backsolve(L, forwardsolve(t(L), y - mean(y))); 0.5 * sum((y - mean(y)) * a) + sum(log(diag(L))) }
  p <- exp(optim(c(0, 0, -1), nll)$par)
  K <- p[2]^2 * exp(-outer(x, x, "-")^2 / (2 * p[1]^2)) + diag(p[3]^2, length(x)); Ks <- p[2]^2 * exp(-outer(xs, x, "-")^2 / (2 * p[1]^2))
  Ki <- solve(K); m <- mean(y) + Ks %*% Ki %*% (y - mean(y)); v <- p[2]^2 - rowSums((Ks %*% Ki) * Ks)
  tibble(x = xs, mean = as.vector(m), lo = as.vector(m) - 1.645 * sqrt(pmax(v, 0)), hi = as.vector(m) + 1.645 * sqrt(pmax(v, 0)), model = "Gaussian process") }
bands <- bind_rows(bart_band, gp_fit(x, yv, xg)) |> mutate(truth = f(x))
ggplot(bands, aes(x)) + geom_ribbon(aes(ymin = lo, ymax = hi), fill = harbour[["sky"]], alpha = 0.25) +
  geom_line(aes(y = mean), colour = harbour[["sky"]], linewidth = 0.9) + geom_line(aes(y = truth), linetype = 2, colour = harbour[["ink"]]) +
  geom_point(data = tibble(x = x, y = yv), aes(x, y), size = 0.7, alpha = 0.5, colour = harbour[["muted"]], inherit.aes = FALSE) +
  geom_vline(xintercept = 1, colour = harbour[["grid"]]) + facet_wrap(~ model) +
  coord_cartesian(ylim = c(-2, 9)) + labs(x = "confounder x", y = "outcome under treatment")

Neither prior can know what happens where there are no data; each answers according to its own assumptions. BART's sum of trees is constant beyond the last split, so it extrapolates flat, and because its uncertainty comes mostly from the data it has seen, its bands stay too narrow. The Gaussian process returns to its prior mean, with uncertainty growing back to the prior's, which is more honest but still a prior, not knowledge. Wherever treated and control units do not overlap, a nonparametric model's estimate of $\tau (x)$ is an extrapolation dressed as an estimate.

15.7 Regularisation-induced confounding

Flexible models need regularisation, and in causal inference regularisation can do harm in a way it cannot in prediction. Richard Hahn, Jared Murray and Carlos Carvalho (2020) identified the mechanism. Suppose treatment is targeted: units with high expected outcomes are more likely to be treated, so the propensity score $\pi (x)=Pr(T=1\mid x)$ rises with $\mu (x,0)$. A BART model with $T$ as just another covariate must learn $\mu (x,0)$ from the covariates, but its prior shrinks toward simple functions, and when a simple function of $x$ fits less well than $T$, the model attributes to treatment what belongs to the confounding: regularisation-induced confounding. The cure is to give the model the propensity score as an extra covariate, so that the confounding is easy to represent, and, better still, to separate the prognostic surface from the treatment effect:

with separate BART priors on $\mu$ and $\tau$. This is the Bayesian causal forest (BCF). The prior on $\tau$ uses fewer, shallower trees, so it shrinks toward a homogeneous effect, and heterogeneity has to earn its place. We compare the three approaches on simulated data with targeted selection and a constant true effect of 3, following the design of Hahn, Murray and Carvalho (2020), over 30 data sets of 250 units. The propensity score is itself estimated by BART, so that it can capture the nonlinear targeting.

ric_path <- "fits/ch15_ric.rds"
sim_ric <- function(n = 250) {
  x <- matrix(rnorm(n * 3), n, 3); x4 <- rbinom(n, 1, 0.5); x5 <- sample(1:3, n, replace = TRUE)
  mu <- -6 + c(2, -1, -4)[x5] + 6 * abs(x[, 3] - 1)                 # prognostic surface
  pi <- 0.8 * pnorm(3 * mu / sd(mu) - 0.5 * x[, 1]) + 0.05 + runif(n) / 10   # targeted selection
  t <- rbinom(n, 1, pi); y <- mu + 3 * t + rnorm(n)                  # constant effect of 3
  list(X = data.frame(x1 = x[, 1], x2 = x[, 2], x3 = x[, 3], x4 = x4, x5 = x5), t = t, y = y, pi = pi) }
fit_ric <- function(d, seed) {
  set.seed(seed)
  ps <- colMeans(pnorm(dbarts::bart(as.matrix(d$X), d$t, ndpost = 500, nskip = 250, verbose = FALSE)$yhat.train))  # BART propensity score
  X1 <- cbind(d$X, t = d$t)
  ate_bart <- function(Xtr) { f <- stochtree::bart(X_train = Xtr, y_train = d$y, num_gfr = 10, num_burnin = 200, num_mcmc = 1000,
      general_params = list(random_seed = seed))
    p1 <- predict(f, X = transform(Xtr, t = 1))$y_hat; p0 <- predict(f, X = transform(Xtr, t = 0))$y_hat; colMeans(p1 - p0) }
  b <- stochtree::bcf(X_train = d$X, Z_train = d$t, y_train = d$y, propensity_train = ps, num_gfr = 10, num_burnin = 200, num_mcmc = 1000,
                      general_params = list(random_seed = seed))
  list(BART = ate_bart(X1), `BART with propensity` = ate_bart(cbind(X1, ps = ps)), BCF = colMeans(b$tau_hat_train)) }
if (file.exists(ric_path)) ric <- readRDS(ric_path) else {
  set.seed(15); ric <- map(1:30, function(r) fit_ric(sim_ric(), r)); saveRDS(ric, ric_path) }
map_dfr(ric, function(r) imap_dfr(r, ~ tibble(method = .y, estimate = mean(.x), lo = quantile(.x, 0.05), hi = quantile(.x, 0.95)))) |>
  summarise(mean_estimate = mean(estimate), bias = mean(estimate - 3), mcse_bias = sd(estimate) / sqrt(n()),
            rmse = sqrt(mean((estimate - 3)^2)), coverage_90 = mean(lo < 3 & hi > 3), interval_width = mean(hi - lo), .by = method) |>
  mutate(across(where(is.double), ~ round(.x, 3))) |> as.data.frame()

                method mean_estimate   bias mcse_bias  rmse coverage_90
1                 BART         3.096  0.096     0.030 0.186       0.933
2 BART with propensity         2.902 -0.098     0.033 0.201       0.933
3                  BCF         2.931 -0.069     0.035 0.200       0.867
  interval_width
1          0.649
2          0.752
3          0.721

All three methods recover the effect of 3 to within about 0.1 on average, with root mean squared errors of about 0.2. Plain BART overestimates it by 0.096 (Monte Carlo standard error 0.030): a small bias, but in the direction regularisation-induced confounding predicts, since part of the confounding has been credited to treatment. Adding the BART-estimated propensity score reverses the sign of the bias (−0.098), and BCF's bias is the smallest in size (−0.069), though with 30 data sets the three are not clearly distinguishable, and neither are their coverages (93%, 93% and 87%, each with a Monte Carlo standard error of about 5 points). With 250 units and a modern BART sampler, the prognostic surface in this design is learned well enough that the problem is mild. The mechanism bites hardest when the prognostic surface is hard to learn relative to the sample size and selection is strongly targeted on it; there, separating $\mu$ from $\tau$ and giving the model the propensity score matter more. The safe default is to do both, since they cost little when they are not needed.

15.8 Bayesian causal forests and the mindset study

A growth mindset is the belief that intellectual abilities can be developed. The National Study of Learning Mindsets randomised about 12,000 ninth-graders in a national probability sample of 65 US public high schools to a short online growth-mindset intervention or a control programme, and measured their grades at the end of the year. Its main questions were about heterogeneity: for whom, and in which schools, does such a cheap intervention work? To let methodologists study those questions without releasing the study's confidential data, its investigators created a synthetic data set with the same structure and similar effect sizes, about 10,000 students in 76 schools, used at the 2018 Atlantic Causal Inference Conference (Carvalho et al., 2019). Unlike the real study, the synthetic data mimic an observational study: treatment depends on the covariates. The data may be used only to study methods for treatment effects.

f <- "data/nlsm/synthetic_data.csv"
if (!file.exists(f)) { dir.create("data/nlsm", showWarnings = FALSE)
  download.file("https://raw.githubusercontent.com/grf-labs/grf/master/experiments/acic18/synthetic_data.csv", f) }
ms <- read_csv(f, show_col_types = FALSE)
# S3 expected success; C1 race; C2 gender; C3 first-generation; XC urbanicity; X1 school mindset norms;
# X2 school achievement; X3 minority share; X4 poverty; X5 school size
X <- ms |> transmute(S3, C1 = factor(C1), C2 = factor(C2), C3, XC = factor(XC), X1, X2, X3, X4, X5)
Xm <- as.data.frame(model.matrix(~ . - 1, X))
ps_ms <- glm(ms$Z ~ ., data = Xm, family = binomial)$fitted.values
c(students = nrow(ms), schools = n_distinct(ms$schoolid), treated_share = mean(ms$Z),
  naive_difference = mean(ms$Y[ms$Z == 1]) - mean(ms$Y[ms$Z == 0]), propensity_range = range(ps_ms)) |> round(3)
bcf_path <- "fits/ch15_bcf_mindset.rds"
if (file.exists(bcf_path)) tau_ms <- readRDS(bcf_path) else {
  fit_ms <- stochtree::bcf(X_train = Xm, Z_train = ms$Z, y_train = ms$Y, propensity_train = ps_ms,
                           num_gfr = 10, num_burnin = 500, num_mcmc = 1000, general_params = list(random_seed = 15))
  tau_ms <- fit_ms$tau_hat_train; saveRDS(tau_ms, bcf_path) }
ate_ms <- colMeans(tau_ms)
c(ATE = mean(ate_ms), q05 = unname(quantile(ate_ms, 0.05)), q95 = unname(quantile(ate_ms, 0.95))) |> round(3)

         students           schools     treated_share  naive_difference 
        10391.000            76.000             0.326             0.304 
propensity_range1 propensity_range2 
            0.132             0.534 
  ATE   q05   q95 
0.249 0.231 0.266 

The synthetic study has 10,391 students in 76 schools, a third of them treated. Because treatment depends on the covariates, the naive difference in mean outcomes, 0.304 standard deviations, is biased; the estimated propensity scores range from 0.13 to 0.53, so every kind of student has a reasonable chance of either arm, and overlap is good. The Bayesian causal forest estimates the average effect at 0.249 (90% interval 0.231 to 0.266). The forest also gives a posterior for the effect of every student, and averaging those within schools gives a posterior for each school's effect.

school <- ms |> mutate(row = row_number()) |> group_by(schoolid) |>
  summarise(X1 = first(X1), X2 = first(X2), rows = list(row), .groups = "drop") |>
  mutate(draws = map(rows, ~ colMeans(tau_ms[.x, , drop = FALSE])),
         mean = map_dbl(draws, mean), lo = map_dbl(draws, ~ quantile(.x, 0.05)), hi = map_dbl(draws, ~ quantile(.x, 0.95)))
school |> select(X1, X2, mean, lo, hi) |> pivot_longer(c(X1, X2), names_to = "moderator", values_to = "value") |>
  mutate(moderator = recode(moderator, X1 = "school fixed-mindset norms (X1)", X2 = "school achievement level (X2)")) |>
  mutate(moderator = factor(moderator, c("school achievement level (X2)", "school fixed-mindset norms (X1)"))) |>
  ggplot(aes(value, mean)) + geom_linerange(aes(ymin = lo, ymax = hi), colour = harbour[["muted"]], alpha = 0.35) +
  geom_point(colour = harbour[["sky"]], size = 1.3) + geom_smooth(method = "loess", se = FALSE, colour = harbour[["rose"]], linewidth = 0.8, formula = y ~ x) +
  facet_wrap(~ moderator, scales = "free_x") + labs(x = "school characteristic (standardised)", y = "average effect in school")

lowhigh <- function(v) { cut_v <- median(v); hi <- school$draws[v > cut_v]; lo <- school$draws[v <= cut_v]
  d <- Reduce(`+`, hi) / length(hi) - Reduce(`+`, lo) / length(lo); c(mean = mean(d), q05 = unname(quantile(d, 0.05)), q95 = unname(quantile(d, 0.95)), p_positive = mean(d > 0)) }
rbind(`high minus low achievement (X2)` = lowhigh(school$X2), `high minus low fixed mindset (X1)` = lowhigh(school$X1)) |> round(3)

                                    mean    q05    q95 p_positive
high minus low achievement (X2)    0.023 -0.005  0.056      0.909
high minus low fixed mindset (X1) -0.072 -0.105 -0.033      0.000

The schools' effects range from about 0.1 to 0.37 standard deviations, and the intervals overlap a good deal. Comparing schools above and below the median of each moderator, draw by draw, gives the posterior of the difference in average effects, which is what a claim of moderation needs. Schools with stronger fixed-mindset norms have smaller effects, by 0.072 (90% interval 0.033 to 0.105), with posterior probability essentially 1 that the difference is negative. Schools with higher achievement have slightly larger effects, by 0.023, but the interval (−0.005 to 0.056) includes zero and the probability of a positive difference is 0.91: suggestive, not established. A median split is a blunt summary; the figure shows the whole relationship, and the forest's prior, which favours homogeneous effects, means that the variation it does report has survived shrinkage.

Lecture 37 — Nonparametric is not assumption-free

15.9 A competition with known answers

How well do flexible Bayesian methods recover causal effects when nobody knows the right model? The 2016 Atlantic Causal Inference Conference ran a competition to find out (Dorie et al., 2019). The organisers took the covariates of a real study, 4,802 infants and 58 variables from the Collaborative Perinatal Project, and simulated treatment assignments and outcomes from 77 different data-generating processes that varied the degree of nonlinearity, the overlap between groups, the alignment between the assignment mechanism and the outcome, and the heterogeneity of effects, 100 data sets each. Thirty teams submitted methods, and because the truth was known, they could be scored. The simulation code is public as the R package aciccomp2016.

We repeat a small version of the competition for BART, the method that did best among the simplest submissions: ten of the 77 settings, one data set each, estimating the conditional average effect on the treated (the mean of $\mu (x,1)-\mu (x,0)$ over the treated units) by g-computation, with and without an estimated propensity score as an extra covariate.

# remotes::install_github("vdorie/aciccomp/2016")    # the competition's simulation package
library(aciccomp2016)
acic_path <- "fits/ch15_acic.rds"
settings <- c(1, 9, 17, 23, 31, 39, 47, 55, 63, 71)
if (file.exists(acic_path)) acic <- readRDS(acic_path) else {
  acic <- map_dfr(settings, function(s) {
    d <- dgp_2016(input_2016, s, 1); X <- model.matrix(~ . - 1, input_2016); tr <- d$z == 1
    truth <- mean(d$mu.1[tr] - d$mu.0[tr])
    ps <- colMeans(pnorm(bart(X, d$z, ndpost = 500, nskip = 250, verbose = FALSE)$yhat.train))     # BART propensity score
    est <- function(Xa) { f <- bart(cbind(Xa, z = d$z), d$y, x.test = cbind(Xa, z = 1 - d$z)[tr, ], ndpost = 500, nskip = 250, verbose = FALSE)
      att <- rowMeans(f$yhat.train[, tr] - f$yhat.test); c(est = mean(att), lo = unname(quantile(att, 0.025)), hi = unname(quantile(att, 0.975))) }
    bind_rows(tibble(setting = s, method = "BART", truth = truth, !!!est(X)),
              tibble(setting = s, method = "BART with propensity", truth = truth, !!!est(cbind(X, ps = ps))),
              tibble(setting = s, method = "difference in means", truth = truth, est = mean(d$y[tr]) - mean(d$y[!tr]), lo = NA, hi = NA)) })
  saveRDS(acic, acic_path) }
acic |> summarise(mean_abs_error = mean(abs(est - truth)), rmse = sqrt(mean((est - truth)^2)),
                  coverage_95 = mean(lo < truth & hi > truth), mean_width = mean(hi - lo), .by = method) |>
  mutate(across(where(is.double), ~ round(.x, 3)))

# A tibble: 3 x 5
  method               mean_abs_error  rmse coverage_95 mean_width
  <chr>                         <dbl> <dbl>       <dbl>      <dbl>
1 BART                          0.043 0.052         0.9      0.17 
2 BART with propensity          0.028 0.036         1        0.179
3 difference in means           0.886 1.06         NA       NA    

Both BART analyses are far better than the difference in means, whose error averages 0.89 because treatment in these simulations is strongly confounded. Plain BART misses the true effect by 0.043 on average, and 9 of its 10 intervals contain it; adding the propensity score cuts the average error to 0.028 and all 10 intervals cover, with intervals only slightly wider (0.179 against 0.170). Ten data sets are a small competition, and one miss in ten is consistent with 95% coverage, but the direction of the comparison is the one Dorie and colleagues found across 7,700 data sets: a flexible outcome model does most of the work, and giving it the propensity score makes its estimates and intervals more trustworthy.

15.10 Overlap, extrapolation and coverage

Lecture 34 drew a distinction between calibration averaged over a prior and calibration for every fixed truth. Nonparametric Bayesian methods are attractive because their priors are wide enough to contain almost any truth, so one might hope for both. Benchmarks such as ACIC show that this does not come free. A prior on functions is wide, but not uniform: BART prefers sums of shallow trees, Gaussian processes prefer smooth functions of a certain length-scale, and BCF prefers homogeneous effects. Where the data are plentiful the prior's preferences hardly matter; where they are thin, which in causal inference means where treated and control units do not overlap, the preferences become the answer, and posterior intervals inherit their confidence.

Three practices follow. First, examine overlap before modelling, by plotting propensity scores, and restrict the estimand to the region of overlap (the effect on the treated, or on a trimmed population) when overlap fails. Second, include the propensity score in the outcome model, or use BCF, which does so by design. Third, check coverage by simulation in settings that resemble your data, by generating outcomes from a model fitted to your data with known effects and refitting. It is the calibrated-Bayes recipe of Lecture 34, and it is the only way to know how much the posterior interval can be trusted.

15.11 Nonparametric Bayes in practice

15.12 Exercises

Each exercise has a hint and an answer button. Code for the answers is in R/ex15.R.

Further reading and references

1. Antoniak, C. E. (1974). Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems. Annals of Statistics, 2(6), 1152–1174.
2. Blackwell, D. and MacQueen, J. B. (1973). Ferguson distributions via Pólya urn schemes. Annals of Statistics, 1(2), 353–355.
3. Carvalho, C., Feller, A., Murray, J., Woody, S. and Yeager, D. (2019). Assessing treatment effect variation in observational studies: results from a data challenge. Observational Studies, 5(2), 21–35.
4. Chipman, H. A., George, E. I. and McCulloch, R. E. (2010). BART: Bayesian additive regression trees. Annals of Applied Statistics, 4(1), 266–298.
5. Dorie, V., Hill, J., Shalit, U., Scott, M. and Cervone, D. (2019). Automated versus do-it-yourself methods for causal inference: lessons learned from a data analysis competition. Statistical Science, 34(1), 43–68.
6. Escobar, M. D. and West, M. (1995). Bayesian density estimation and inference using mixtures. Journal of the American Statistical Association, 90(430), 577–588.
7. Ferguson, T. S. (1973). A Bayesian analysis of some nonparametric problems. Annals of Statistics, 1(2), 209–230.
8. Fong, E., Holmes, C. and Walker, S. G. (2023). Martingale posterior distributions. Journal of the Royal Statistical Society: Series B, 85(5), 1357–1391.
9. Ghosal, S. and van der Vaart, A. (2017). Fundamentals of Nonparametric Bayesian Inference. Cambridge University Press.
10. Hahn, P. R., Murray, J. S. and Carvalho, C. M. (2020). Bayesian regression tree models for causal inference: regularization, confounding, and heterogeneous effects. Bayesian Analysis, 15(3), 965–1056.
11. Hill, J. L. (2011). Bayesian nonparametric modeling for causal inference. Journal of Computational and Graphical Statistics, 20(1), 217–240.
12. Hill, J., Linero, A. and Murray, J. (2020). Bayesian additive regression trees: a review and look forward. Annual Review of Statistics and Its Application, 7, 251–278.
13. Miller, J. W. and Harrison, M. T. (2013). A simple example of Dirichlet process mixture inconsistency for the number of components. Advances in Neural Information Processing Systems, 26.
14. Miller, J. W. and Harrison, M. T. (2018). Mixture models with a prior on the number of components. Journal of the American Statistical Association, 113(521), 340–356.
15. Müller, P., Quintana, F. A., Jara, A. and Hanson, T. (2015). Bayesian Nonparametric Data Analysis. Springer.
16. Oganisian, A. and Roy, J. A. (2021). A practical introduction to Bayesian estimation of causal effects: parametric and nonparametric approaches. Statistics in Medicine, 40(2), 518–551.
17. Rubin, D. B. (1981). The Bayesian bootstrap. Annals of Statistics, 9(1), 130–134.
18. Sethuraman, J. (1994). A constructive definition of Dirichlet priors. Statistica Sinica, 4(2), 639–650.
19. Yeager, D. S., Hanselman, P., Walton, G. M., et al. (2019). A national experiment reveals where a growth mindset improves achievement. Nature, 573, 364–369.