Chapter 13
Causal inference in depth

Lectures 26–29 · potential outcomes, the back door, what not to adjust for, instruments, the front door, flexible structural equations, sensitivity analysis, and the road to causal AI

Every model in this course has been used to answer some kind of "what if" question, and most of them were careful about it: Chapter 1 drew the vaccine trial's causal graph, Chapter 3 separated total and direct effects of sex on weight, and Chapter 4 asked what a risk score's errors mean for the people it scores. But the machinery of the last three chapters (prediction with many features, clustering, classifying text and images) lives on what Judea Pearl calls the bottom rung of the ladder of causation (Pearl and Mackenzie, 2018). The first rung is association: seeing, and predicting what we will see. The second is intervention: what will happen if we do something. The third is counterfactuals: what would have happened to this person had things been different. No amount of data on the first rung answers questions on the second without assumptions, and the course has already met the consequences: the Table 2 fallacy of Chapter 10, and Chapter 12's detector of AI-written text, which learned how one data set was scraped rather than how ChatGPT writes.

This chapter climbs the ladder properly. It is the longest in the course because it carries two jobs. It completes the causal toolkit that Chapters 1 and 2 began, and it prepares the final chapter, which asks how the same ideas can make artificial intelligence itself model causes rather than correlations.

• Lecture 26 is about the problem. Potential outcomes and the fundamental problem of causal inference; a job-training programme whose observational evaluation gives the wrong sign; the back-door criterion.
• Lecture 27 is about adjustment. What not to condition on; Bayesian g-computation; the difference between the average effect and the effect on the treated; overlap.
• Lecture 28 is about unmeasured confounding. Instrumental variables with Card's college-proximity data; mediation; the front-door criterion.
• Lecture 29 is about flexibility and humility. Structural equations with Bayesian trees, scored against a benchmark whose true effects are known; sensitivity analysis; and what all of this means for AI.

Lecture 26 — The problem

13.1 Potential outcomes

Take one worker and one job-training programme. The worker has two potential outcomes: $Y_i(1)$, what she would earn next year if she took the training, and $Y_i(0)$, what she would earn if she did not. The causal effect of training for her is the difference $Y_i(1)-Y_i(0)$. The fundamental problem of causal inference (Holland, 1986) is that we only ever see one of the two: if she trains, $Y_i(0)$ is a counterfactual that never happens. Individual effects are not observed, for anyone, ever.

What can be estimated, under assumptions, are averages. The average treatment effect is $\text{ATE}=𝔼[Y(1)-Y(0)]$ over a population; the average treatment effect on the treated is $\text{ATT}=𝔼[Y(1)-Y(0)\mid T=1]$, the average over those who actually took the treatment. They differ whenever the effect differs between people and people with different effects choose differently. In Pearl's notation, $𝔼[Y\mid \text{do}(T=1)]$ is the average of $Y(1)$ over everyone, the outcome if we set $T=1$ by intervention, and it is generally not the same as $𝔼[Y\mid T=1]$, the average among those who happened to have $T=1$.

The naive comparison of treated and untreated averages estimates the ATE only if treatment is independent of the potential outcomes. Randomisation guarantees this; self-selection does not. The widget makes the problem concrete with twelve workers whose two potential outcomes we, unusually, know.

Under randomisation the naive difference scatters around the true effect, wider for twelve workers than for a thousand but centred correctly. When the poorest enrol, the treated group would have earned less anyway, and the naive difference is biased downward, sometimes enough to change its sign. When those who gain most enrol, the naive difference estimates neither the ATE nor the ATT reliably, and the ATT itself exceeds the ATE. These are the three worlds of observational data, and only the first comes with a guarantee.

13.2 A job-training programme and the wrong answer

In the mid-1970s the National Supported Work (NSW) demonstration offered disadvantaged workers, many of them recently released from prison, recovering from addiction or long unemployed, up to eighteen months of subsidised work. Eligible applicants were randomly assigned to the programme or to a control group, and their earnings in 1978 were recorded. Robert LaLonde (1986) used the experiment as a test of observational methods: he threw away the randomised control group, replaced it with men from national surveys who had not taken the programme, and asked whether the econometric methods of the day could recover the experimental answer from the resulting observational comparison. They could not, and his paper changed applied economics. Rajeev Dehejia and Sadek Wahba (1999) later reanalysed a subset with propensity scores. We use their subset: the 185 trained men, compared with 260 randomised controls (the experiment) or with 429 men from the Panel Study of Income Dynamics (PSID) who did not train (the observational study).

library(tidyverse)
library(brms)
library(posterior)
library(jsonlite)
library(dbarts)
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())

gh <- "https://raw.githubusercontent.com/cran/"
lalonde <- read.table(paste0(gh, "MatchIt/master/data/lalonde.tab"), header = TRUE) |>   # NSW trainees + PSID comparison men
  transmute(treat, age, educ, black = as.integer(race == "black"), hisp = as.integer(race == "hispan"),
            married, nodegree, re74 = re74 / 1000, re75 = re75 / 1000, re78 = re78 / 1000)   # earnings in $1,000s
tmp <- tempfile(); download.file(paste0(gh, "causaldata/master/data/nsw_mixtape.rda"), tmp, quiet = TRUE, mode = "wb")
load(tmp)                                                                   # nsw_mixtape: the randomised experiment
experiment <- nsw_mixtape |> transmute(treat, re78 = re78 / 1000)
lalonde |> group_by(treat) |> summarise(n = n(), age = mean(age), educ = mean(educ), black = mean(black),
  married = mean(married), re74 = mean(re74), re75 = mean(re75), re78 = mean(re78)) |> mutate(across(age:re78, ~ round(.x, 2)))
c(naive_observational = mean(lalonde$re78[lalonde$treat == 1]) - mean(lalonde$re78[lalonde$treat == 0]),
  experimental_benchmark = mean(experiment$re78[experiment$treat == 1]) - mean(experiment$re78[experiment$treat == 0])) |> round(3)

# A tibble: 2 x 9
  treat     n   age  educ black married  re74  re75  re78
  <int> <int> <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl> <dbl>
1     0   429  28.0  10.2  0.2     0.51  5.62  2.47  6.98
2     1   185  25.8  10.4  0.84    0.19  2.1   1.53  6.35
   naive_observational experimental_benchmark 
                -0.635                  1.794 

The two groups hardly resemble each other. The PSID men are older (28.0 years against 25.8), far more often married (51% against 19%) and far less often Black (20% against 84%), and before the programme they earned much more: $5,620 in 1974 against $2,100, and $2,470 in 1975 against $1,530. Only education is similar (10.2 and 10.4 years). In 1978 the trainees earned $6,350 on average and the PSID men $6,980, so the naive observational comparison says that training cost each man $635. The randomised experiment says the opposite: trainees earned $1,794 more than the randomised controls. The naive comparison gets the sign wrong because the PSID men would have earned more than the trainees with or without training.

13.3 The back door

A back-door path from $T$ to $Y$ is a path that starts with an arrow into $T$. It carries association that is not causal. The back-door criterion (Pearl, 1993; Chapter 1) says that a set of variables $C$ is sufficient to identify the effect of $T$ on $Y$ if no variable in $C$ is a descendant of $T$ and $C$ blocks every back-door path. Then

the adjustment formula: estimate the outcome within each level of the confounders, then average over the confounders' distribution. In potential-outcome language, the same assumption is conditional ignorability, $\{Y(0),Y(1)\}⟂T\mid C$: within groups of identical $C$, treatment is as good as randomised. Two further conditions are needed in practice. Positivity (or overlap): every kind of person must have some chance of each treatment, $0<Pr(T=1\mid C=c)<1$, or the comparison within that kind of person is impossible. And consistency: the treatment must be well defined, so that "training" means the same thing for everyone.

The criterion is as much about what not to adjust for as about what to adjust for. The widget simulates the three structures of Chapter 1 (fork, pipe and collider) attached to a treatment-effect question.

With a confounder, adjustment removes the bias. With a mediator, the regression on $T$ alone recovers the total effect, and adjusting for $Z$ removes the part of the effect that flows through $Z$: the estimate is the direct effect, which is a different question. With a collider, the unadjusted regression is right and adjustment creates a bias, by the mechanism of Berkson's paradox (Chapter 1). Regression output looks identical in all three cases; only the graph says which number answers the question.

Lecture 27 — Adjustment

13.4 What not to condition on

Three classes of variables do harm when added to an adjustment set, and all three appear regularly in published work (Cinelli, Forney and Pearl, 2024).

• Post-treatment variables, including mediators. Anything affected by the treatment, such as employment during the programme or job satisfaction in 1978, absorbs part of the effect, or opens a collider path between treatment and an unmeasured cause of the outcome.
• Colliders and their descendants. Selecting the sample on a common effect of treatment and outcome (studying only people who were still in the area in 1978, or only patients admitted to hospital) creates association where there was none.
• Pure predictors of treatment that do not affect the outcome, the instruments of Section 13.6. Adjusting for them does not bias an estimate when there is no unmeasured confounding, but it removes useful variation in the treatment and inflates the variance; when there is unmeasured confounding, it amplifies the bias.

Pearl (2011) calls the last effect bias amplification. The opposite class deserves mention: variables that affect only the outcome are good controls. They are not needed for identification, but they reduce the residual variance and sharpen the estimate. Prior earnings in the job-training study are both a confounder and a strong predictor of later earnings, which is why they matter so much.

13.5 Bayesian g-computation

The adjustment formula suggests an algorithm, which James Robins (1986) called the g-formula or g-computation. Fit a model for the outcome given treatment and confounders. For every unit, predict the outcome twice: once with $T$ set to 1 and once with $T$ set to 0, keeping its own confounders. Average the difference over the units of interest: all units for the ATE, the treated units for the ATT. In a Bayesian fit, do this for every posterior draw, and the result is a posterior for the causal effect, with the outcome model's uncertainty carried through.

A note on the target. The NSW experiment randomised men who had applied to the programme, so the experimental benchmark is an average effect for the kind of men who take this programme: in effect, an ATT. The PSID men are a different population (older, more often married, with nearly three times the 1974 earnings), and the average effect of training over them, the ATE in our pooled sample, is a different and less meaningful quantity. The right observational estimate to compare with the benchmark is the ATT.

We fit two outcome models in brms, with earnings in thousands of dollars. The first is the classical linear adjustment, with a single treatment effect for everyone. The second allows the effect of training to vary with every confounder (treatment-by-covariate interactions), so that ATE and ATT can differ.

f_lin <- re78 ~ treat + age + educ + black + hisp + married + nodegree + re74 + re75
f_int <- re78 ~ treat * (age + educ + black + hisp + married + nodegree + re74 + re75)
pri <- c(prior(normal(5, 10), class = Intercept), prior(normal(0, 5), class = b), prior(exponential(0.2), class = sigma))
m_lin <- brm(f_lin, data = lalonde, prior = pri, seed = 13, refresh = 0, file = "fits/ch13_lin")
m_int <- brm(f_int, data = lalonde, prior = pri, seed = 13, refresh = 0, file = "fits/ch13_int")
g_comp <- function(m, rows = TRUE) {           # posterior of the average of Y(1) - Y(0) over the chosen units
  d1 <- mutate(lalonde, treat = 1)[rows, ]; d0 <- mutate(lalonde, treat = 0)[rows, ]
  rowMeans(posterior_epred(m, newdata = d1) - posterior_epred(m, newdata = d0))
}
effects <- list(linear = g_comp(m_lin), interaction_ATE = g_comp(m_int), interaction_ATT = g_comp(m_int, lalonde$treat == 1))
map_dfr(effects, ~ c(mean = mean(.x), q05 = unname(quantile(.x, 0.05)), q95 = unname(quantile(.x, 0.95))), .id = "estimate") |>
  mutate(across(-estimate, ~ round(.x, 2)))

# A tibble: 3 x 4
  estimate         mean   q05   q95
  <chr>           <dbl> <dbl> <dbl>
1 linear           1.51  0.25  2.76
2 interaction_ATE  1.11 -0.51  2.77
3 interaction_ATT  1.63  0.34  2.97

Linear adjustment moves the estimate from −$635 to $1,510 (90% interval $250 to $2,760), and the interaction model's ATT is $1,630 ($340 to $2,970), both close to the experimental $1,794. The interaction model's ATE, averaged over all 614 men, is lower, $1,110, and its interval (−$510 to $2,770) includes zero: the model extrapolates the effect of training to well-paid PSID men, who would never take such a programme and for whom the data say little. All the intervals are wide, $2,500 to $3,300 from end to end. With 185 trainees and noisy earnings, even a correct adjustment leaves a great deal of uncertainty, and a single number such as "$1,510" would hide it.

Adjustment can only compare like with like where like exists. The propensity score, $e(c)=Pr(T=1\mid C=c)$ (Rosenbaum and Rubin, 1983), summarises how comparable the groups are: if the trainees and the PSID men overlap in their propensity scores, every trainee has comparable untrained men; if they do not, the outcome model is extrapolating. We fit a logistic propensity model and look.

ps_fit <- glm(treat ~ age + I(age^2) + educ + black + hisp + married + nodegree + re74 + re75 + I(re74 == 0) + I(re75 == 0),
              family = binomial, data = lalonde)
lalonde$ps <- fitted(ps_fit)
ggplot(lalonde, aes(ps, fill = factor(treat, labels = c("PSID comparison", "trainees")))) +
  geom_histogram(bins = 40, position = "identity", alpha = 0.6, boundary = 0) +
  scale_fill_manual(values = unname(harbour[c("sky", "rose")])) + labs(x = "propensity score", y = "men", fill = NULL)
lalonde |> group_by(treat) |> summarise(below_0.05 = sum(ps < 0.05), between = sum(ps >= 0.05 & ps <= 0.95), above_0.95 = sum(ps > 0.95))

# A tibble: 2 x 4
  treat below_0.05 between above_0.95
  <int>      <int>   <int>      <int>
1     0        216     213          0
2     1          6     172          7

Overlap is poor. Half the PSID men, 216 of 429, have a propensity score below 0.05: by their age, marriage and earnings, they are men for whom the model sees almost no chance of training, and there are almost no trainees like them. None has a score above 0.95. On the other side, 7 trainees have scores above 0.95 and 6 below 0.05. Any estimate that averages over the PSID men, such as the ATE, rests on extrapolation for half the sample; the ATT needs comparators only for the trainees, and most trainees have some. This is the second reason, after the meaning of the benchmark, to target the ATT here.

naive <- mean(lalonde$re78[lalonde$treat == 1]) - mean(lalonde$re78[lalonde$treat == 0])
bench <- mean(experiment$re78[experiment$treat == 1]) - mean(experiment$re78[experiment$treat == 0])
labels <- c(linear = "linear, constant effect", interaction_ATE = "interactions: ATE", interaction_ATT = "interactions: ATT")
map_dfr(effects, ~ tibble(effect = .x), .id = "estimate") |> mutate(estimate = labels[estimate]) |>
  ggplot(aes(effect, colour = estimate)) + geom_density(linewidth = 0.9) +
  geom_vline(xintercept = bench, linetype = 2, colour = harbour[["moss"]]) +
  geom_vline(xintercept = naive, linetype = 3, colour = harbour[["muted"]]) +
  scale_colour_manual(values = unname(harbour[c("rose", "amber", "sky")])) +
  labs(x = "effect of training on 1978 earnings ($1,000s)", y = NULL, colour = NULL)

Lecture 28 — When the back door is closed

13.6 Instruments

Back-door adjustment needs every confounder to be measured. Often one is not: ability, motivation, health, the things that make people both seek a treatment and do well. An instrumental variable $Z$ offers a way round. It must satisfy three conditions: relevance (it affects the treatment), exclusion (it affects the outcome only through the treatment) and independence (it shares no common cause with the outcome). Then the instrument provides variation in the treatment that is unrelated to the unmeasured confounder, and the effect can be recovered from it.

In the simplest linear case, with one instrument and no covariates, the effect is the ratio of two regression slopes, the effect of $Z$ on $Y$ divided by the effect of $Z$ on $T$: $\beta =\text{Cov}(Z,Y)/\text{Cov}(Z,T)$. The Bayesian version models the two equations jointly,

with the correlation between $u$ and $v$ left free. That correlation is the footprint of the unmeasured confounder: modelling it, instead of assuming it is zero as ordinary regression does, is what lets the joint model separate the causal effect from the confounding.

David Card (1995) used this idea to estimate the return to schooling. Men with more education earn more, but ability and family background affect both. Card's instrument was whether a man grew up near a four-year college: proximity makes college cheaper and so increases schooling (relevance), and, Card argued, affects wages only through schooling, after accounting for region and urban residence (exclusion and independence). The data, from the US National Longitudinal Survey of Young Men, are 3,010 men interviewed in 1976. We fit ordinary regression and the joint instrumental-variable model, with the same controls: experience and its square, race, region and urban residence in 1966 and 1976.

tmp <- tempfile(); download.file(paste0(gh, "wooldridge/master/data/card.RData"), tmp, quiet = TRUE, mode = "wb")
load(tmp)                                                                    # card: NLS Young Men, 1976
card <- card |> mutate(across(c(exper, expersq), ~ as.numeric(scale(.x))))
controls <- "exper + expersq + black + south + smsa + smsa66 + reg662 + reg663 + reg664 + reg665 + reg666 + reg667 + reg668 + reg669"
m_ols <- brm(as.formula(paste("lwage ~ educ +", controls)), data = card, seed = 13, refresh = 0, file = "fits/ch13_ols")
f_iv <- bf(as.formula(paste("educ ~ nearc4 +", controls))) +               # first stage: schooling on the instrument
        bf(as.formula(paste("lwage ~ educ +", controls))) + set_rescor(TRUE) # outcome; correlated errors absorb U
m_iv <- brm(f_iv, data = card, seed = 13, refresh = 0, file = "fits/ch13_iv")   # about 9 minutes on two cores
rbind(ordinary = posterior_summary(m_ols, variable = "b_educ", probs = c(0.05, 0.95)),
      iv_return = posterior_summary(m_iv, variable = "b_lwage_educ", probs = c(0.05, 0.95)),
      first_stage = posterior_summary(m_iv, variable = "b_educ_nearc4", probs = c(0.05, 0.95)),
      error_correlation = posterior_summary(m_iv, variable = "rescor__educ__lwage", probs = c(0.05, 0.95))) |> round(3)

                    Estimate Est.Error     Q5   Q95
b_educ                 0.075     0.004  0.069 0.081
b_lwage_educ           0.132     0.060  0.041 0.231
b_educ_nearc4          0.301     0.088  0.156 0.444
rescor__educ__lwage   -0.257     0.246 -0.632 0.163

Ordinary regression puts the return to a year of schooling at 7.5% (90% interval 6.9% to 8.1%). The instrumental-variable model puts it at 13.2%, the same point estimate as Card's, but with a 90% interval from 4.1% to 23.1%, sixteen times as wide, and wide enough to include the ordinary-regression value. Growing up near a four-year college adds 0.30 years of schooling (0.16 to 0.44): a real but modest first stage, and the source of the imprecision. The residual correlation is the surprise. If unmeasured ability were the main confounder, it would raise both schooling and wages, the correlation would be positive and the IV estimate would be below the regression. Instead the IV estimate is above it and the correlation leans negative (−0.26, interval −0.63 to 0.16). Labour economists have offered two explanations: measurement error in reported schooling, which biases ordinary regression downward, and heterogeneous returns, with the compliers (men for whom cost was the obstacle) gaining more from a year of school than the average man. These data cannot tell them apart, and the honest summary is that the instrument is consistent with a return anywhere from modest to large.

The instrumental-variable estimate comes with two caveats that the posterior alone does not show. The first is weakness: when the instrument shifts the treatment only a little, the ratio divides by a small, noisy number, and the estimate is both imprecise and, in finite samples, biased toward the ordinary regression. The second is locality: with heterogeneous effects, an instrument identifies the average effect only for the compliers, the people whose treatment is changed by the instrument (Imbens and Angrist, 1994). Card's estimate is the return to schooling for men whose education depended on whether a college was nearby, plausibly men from families for whom cost was the obstacle, not the return for everyone. The widget explores the first caveat.

13.7 Mediation and the front door

A treatment often works through intermediate steps. Training may raise earnings directly, through skills, and indirectly, by getting a man a first job whose record then helps him get the next. With a mediator $M$ on the path $T\to M\to Y$, the total effect splits into a natural direct effect, the effect of changing $T$ while holding $M$ at the value it would have taken without treatment, and a natural indirect effect, the effect of shifting $M$ from its untreated to its treated value while holding $T$ fixed (Robins and Greenland, 1992; Pearl, 2001). Both are computed from posterior draws of a model for $M$ and a model for $Y$ by the same g-computation as before. The price is an extra assumption: no unmeasured confounding of the mediator and the outcome, which randomising the treatment does not provide.

Mediators can also rescue a study. Suppose an unmeasured $U$ confounds $T$ and $Y$, so the back door cannot be closed, but the whole effect of $T$ on $Y$ passes through a measured mediator $M$, with $T\to M$ unconfounded and $M\to Y$ unconfounded once $T$ is controlled. Then the front-door criterion (Pearl, 1995) identifies the effect:

where

The first factor is the effect of $T$ on $M$, which is unconfounded. The second, $Q(m)$, is the effect of $M$ on $Y$, made unconfounded by adjusting for $T$, which blocks the back-door path $M\leftarrow T\leftarrow U\to Y$. Chaining the two gives the effect of $T$ on $Y$ without ever measuring $U$. Pearl's illustration was smoking ($T$), tar in the lungs ($M$) and cancer ($Y$), with a hypothetical gene ($U$) causing both smoking and cancer: R. A. Fisher's argument against the smoking–cancer link in the 1950s. We simulate such a world, where the truth is known, and estimate the effect three ways.

set.seed(13)
n <- 2000
sim <- tibble(U = rnorm(n), T = 0.8 * U + rnorm(n),            # U is never given to the models below
              M = 0.7 * T + rnorm(n, 0, 0.6), Y = 0.5 * M + 0.9 * U + rnorm(n))
truth <- 0.7 * 0.5                                              # effect of T on Y: through M only
m_naive <- brm(Y ~ T, data = sim, seed = 13, refresh = 0, file = "fits/ch13_fd_naive")
m_fd <- brm(bf(M ~ T) + bf(Y ~ M + T) + set_rescor(FALSE), data = sim, seed = 13, refresh = 0, file = "fits/ch13_fd")
dr <- as_draws_df(m_fd)
front_door <- dr$b_M_T * dr$b_Y_M                               # linear case of the front-door formula
rbind(naive = posterior_summary(as_draws_df(m_naive)$b_T, probs = c(0.05, 0.95)),
      front_door = posterior_summary(front_door, probs = c(0.05, 0.95))) |> round(3)
c(truth = truth)

     Estimate Est.Error    Q5   Q95
[1,]    0.778     0.022 0.741 0.814
[2,]    0.357     0.030 0.307 0.408
truth 
 0.35 

The naive regression of $Y$ on $T$ estimates 0.78 (90% interval 0.74 to 0.81), more than twice the true 0.35, because the unmeasured $U$ raises both. The front-door estimate, built only from the two unconfounded links, is 0.36 (0.31 to 0.41) and contains the truth. Neither model ever saw $U$.

With no direct path, the front-door estimate is centred on the truth however strong the confounding; the naive estimate is not. Add a direct path and the front-door estimate misses it entirely, because it only counts the effect that passes through $M$, so it now underestimates the total effect. Real front-door designs are rare for exactly this reason: few treatments act through a single measured mechanism. The criterion's importance is less in its everyday use than in what it shows, that identification depends on the graph, not on a list of approved methods.

Lecture 29 — Flexibility and humility

13.8 Flexible structural equations

Identification (the back door, an instrument, the front door) says which quantities the data can estimate and what to condition on. It says nothing about the shape of the relationships. Every estimate so far has assumed linear structural equations, and for the job-training data even the interaction model forced each covariate's modification of the effect to be linear.

Underneath every DAG is a structural causal model: one equation per variable, $X_j=f_j(\text{pa}(X_j),U_j)$, saying how each variable is generated from its parents and its own independent noise. An intervention $\text{do}(T=t)$ replaces $T$'s equation with the constant $t$ and leaves the others alone; that is the whole meaning of the do-operator. Nothing requires the $f_j$ to be linear, and the Bayesian toolkit of Chapters 8 and 10 supplies flexible ones with honest uncertainty: Gaussian processes, Bayesian additive regression trees (BART) and Bayesian neural networks. Jennifer Hill (2011) showed that BART, used as the outcome model in g-computation, is a strong default for causal effects, because it captures interactions and nonlinearity without the analyst specifying them and returns a posterior for every unit's effect.

set.seed(13)
X <- as.matrix(select(lalonde, treat, age, educ, black, hisp, married, nodegree, re74, re75))
X1 <- X; X1[, "treat"] <- 1; X0 <- X; X0[, "treat"] <- 0
fit_bart <- bart(X, lalonde$re78, rbind(X1, X0), ndpost = 2000, nskip = 1000, verbose = FALSE)
n_l <- nrow(lalonde)
unit_effects <- fit_bart$yhat.test[, 1:n_l] - fit_bart$yhat.test[, n_l + 1:n_l]   # draws x units
effects$BART_ATT <- rowMeans(unit_effects[, lalonde$treat == 1])
effects$BART_ATE <- rowMeans(unit_effects)
map_dfr(effects, ~ c(mean = mean(.x), q05 = unname(quantile(.x, 0.05)), q95 = unname(quantile(.x, 0.95))), .id = "estimate") |>
  mutate(across(-estimate, ~ round(.x, 2)))

# A tibble: 5 x 4
  estimate         mean   q05   q95
  <chr>           <dbl> <dbl> <dbl>
1 linear           1.51  0.25  2.76
2 interaction_ATE  1.11 -0.51  2.77
3 interaction_ATT  1.63  0.34  2.97
4 BART_ATT         1.44  0.04  2.82
5 BART_ATE         1.21 -0.3   2.68

BART estimates the effect on the trainees at $1,440 (90% interval $40 to $2,820), and the ATE at $1,210 (−$300 to $2,680). These are close to the linear and interaction models: with 614 men and earnings that depend roughly linearly on past earnings, flexibility buys little here, and the honest conclusion is that three reasonable models, and their five estimates, agree on an effect of about $1,100 to $1,600, all below but within error of the experimental $1,794. The BART interval for the ATT only just excludes zero, which Section 13.9 will take seriously.

The job-training data cannot tell us whether a method captured individual effects correctly, because individual effects are never observed. To score methods on that, we need data where the truth is known, which means simulation. Hill (2011) built such a benchmark from a real study, the Infant Health and Development Program (IHDP), a randomised trial of intensive home visits and child care for low-birth-weight premature infants. She kept the 25 real covariates of the children and their mothers, made the treated group non-random by removing a subset of treated children (those with non-white mothers), and simulated the outcomes from known nonlinear functions, so that each child's true individual effect is known. The data set, 747 children of whom 139 are treated, has become a standard test for causal machine learning. We use the first of the replications distributed with the CEVAE paper (Louizos et al., 2017) and compare three outcome models: linear with a constant effect, linear with treatment interactions, and BART. The score is $\sqrt{\text{PEHE}}$, the root mean squared error of the estimated individual effects (precision in estimating heterogeneous effects), and the coverage of their 90% intervals.

ihdp <- read.csv("https://raw.githubusercontent.com/AMLab-Amsterdam/CEVAE/master/datasets/IHDP/csv/ihdp_npci_1.csv",
                 header = FALSE, col.names = c("treat", "y", "y_cf", "mu0", "mu1", paste0("x", 1:25)))
tau <- ihdp$mu1 - ihdp$mu0                                     # the true individual effects (known: simulated)
Xi <- as.matrix(select(ihdp, treat, starts_with("x"))); Xi1 <- Xi; Xi1[, "treat"] <- 1; Xi0 <- Xi; Xi0[, "treat"] <- 0
set.seed(13)
b_ihdp <- bart(Xi, ihdp$y, rbind(Xi1, Xi0), ndpost = 2000, nskip = 1000, verbose = FALSE)
n_i <- nrow(ihdp); td <- b_ihdp$yhat.test[, 1:n_i] - b_ihdp$yhat.test[, n_i + 1:n_i]
lm_add <- lm(y ~ ., data = select(ihdp, treat, y, starts_with("x")))
lm_int <- lm(y ~ treat * ., data = select(ihdp, treat, y, starts_with("x")))
tau_int <- predict(lm_int, mutate(ihdp, treat = 1)) - predict(lm_int, mutate(ihdp, treat = 0))
score <- function(est) c(ATE = mean(est), sqrt_PEHE = sqrt(mean((est - tau)^2)))
rbind(truth = c(ATE = mean(tau), sqrt_PEHE = 0), naive = c(mean(ihdp$y[ihdp$treat == 1]) - mean(ihdp$y[ihdp$treat == 0]), NA),
      linear_constant = score(rep(coef(lm_add)[["treat"]], n_i)), linear_interactions = score(tau_int), BART = score(colMeans(td))) |> round(3)
c(true_effect_sd = sd(tau), BART_coverage_90 = mean(tau > apply(td, 2, quantile, 0.05) & tau < apply(td, 2, quantile, 0.95))) |> round(3)

                      ATE sqrt_PEHE
truth               4.016     0.000
naive               4.021        NA
linear_constant     3.929     0.864
linear_interactions 3.961     0.583
BART                4.018     0.389
  true_effect_sd BART_coverage_90 
           0.860            0.971 

In this replication the true average effect is 4.02, and every method gets the average about right: even the naive difference (4.02) is on target. That is a coincidence of this replication, not a property of the design; Exercise 5 shows the naive error varying from −0.20 to 1.63 across the first ten replications. The individual effects are where the methods separate. The true effects vary (standard deviation 0.86), and a model with one constant effect has $\sqrt{\text{PEHE}}=0.86$, no better than ignoring the variation. Linear treatment interactions reduce the error to 0.58, and BART to 0.39, less than half the constant-effect error. BART's 90% intervals contain 97.1% of the true individual effects: slightly wide, which is the safe direction.

set.seed(1); show <- sample(n_i, 120)
bind_rows(tibble(model = "linear with interactions", true = tau, est = tau_int, lo = NA, hi = NA),
          tibble(model = "BART", true = tau, est = colMeans(td), lo = apply(td, 2, quantile, 0.05), hi = apply(td, 2, quantile, 0.95)) |>
            mutate(lo = replace(lo, -show, NA), hi = replace(hi, -show, NA))) |>
  mutate(model = factor(model, c("linear with interactions", "BART"))) |>
  ggplot(aes(true, est)) + geom_abline(colour = harbour[["muted"]], linetype = 2) +
  geom_linerange(aes(ymin = lo, ymax = hi), colour = harbour[["sky"]], alpha = 0.35, na.rm = TRUE) +
  geom_point(size = 0.8, colour = harbour[["rose"]], alpha = 0.6) + facet_wrap(~ model) +
  labs(x = "true individual effect", y = "estimated effect") + theme(panel.spacing = unit(1.5, "lines"))

The linear model with interactions is too straight for a curved truth: it overshoots for the children with the largest true effects and scatters widely for those with small ones. BART follows the diagonal much more closely, and its intervals are widest for the children with small effects, of whom there are few: the uncertainty is largest where the data are thinnest, as it should be.

13.9 Sensitivity analysis and the limits

No estimate in this chapter came with a guarantee. Back-door adjustment assumes every confounder is measured; instruments assume exclusion; the front door assumes a fully mediating, unconfounded mechanism. None of these can be checked from the data at hand. The honest response is sensitivity analysis: ask how strong an unmeasured confounder would have to be to change the conclusion, and report the answer.

For a linear outcome model there is a simple version. Suppose an unmeasured binary confounder $U$ raises 1978 earnings by $\gamma$ thousand dollars, and is more common among the trainees than among comparable controls by $\delta$ (a difference in proportions). Then the adjusted estimate is too large by about $\gamma \delta$, and the corrected effect is the estimate minus $\gamma \delta$. A Bayesian analysis can put priors on $\gamma$ and $\delta$ and carry them through; the simplest report is the tipping point, the combinations of $\gamma$ and $\delta$ at which the corrected effect's 90% interval first includes zero. (If the confounder lowers earnings, or is less common among trainees, the sign flips and the estimate is too small.) VanderWeele and Ding's (2017) E-value and the robustness value of Cinelli and Hazlett (2020) are refinements of the same idea.

att <- effects$BART_ATT
grid <- expand_grid(gamma = seq(0, 6, by = 0.1), delta = seq(0, 0.6, by = 0.01)) |>
  mutate(mean = mean(att) - gamma * delta, lower = quantile(att, 0.05) - gamma * delta)
ggplot(grid, aes(gamma, delta)) + geom_raster(aes(fill = mean)) +
  geom_contour(aes(z = lower), breaks = 0, colour = harbour[["ink"]], linewidth = 0.8) +
  geom_contour(aes(z = mean), breaks = 0, colour = harbour[["ink"]], linetype = 2) +
  scale_fill_gradient2(low = harbour[["rose"]], mid = harbour[["grid"]], high = harbour[["sky"]], name = "corrected\nATT ($1,000s)") +
  labs(x = "γ: effect of U on 1978 earnings ($1,000s)", y = "δ: excess prevalence of U among trainees")
c(posterior_mean = mean(att), q05 = unname(quantile(att, 0.05)),
  gamma_delta_to_reach_zero_lower = unname(quantile(att, 0.05)), gamma_delta_to_reach_zero_mean = mean(att)) |> round(2)

                 posterior_mean                             q05 
                           1.44                            0.04 
gamma_delta_to_reach_zero_lower  gamma_delta_to_reach_zero_mean 
                           0.04                            1.44 

The result is sobering. The 90% interval for the effect on the trainees already reaches down to $40, so an unmeasured confounder with $\gamma \delta$ of only 0.04 ($40), for example one that adds $200 to 1978 earnings and is 20 percentage points more common among trainees, would make it include zero. For the posterior mean to fall to zero, $\gamma \delta$ would need to reach 1.44: for instance, a factor that adds $3,000 to 1978 earnings and is 48 percentage points more common among trainees than among comparable PSID men. The observational analysis agrees with the experiment in its centre, but on its own it could not have ruled out no effect, and a careful report would say exactly that.

Three other families of designs avoid the unmeasured-confounder argument by exploiting how the data arose, and each deserves a chapter of its own. Difference-in-differences compares the change in outcomes before and after a treatment between treated and untreated groups, assuming parallel trends in the absence of treatment. Regression discontinuity compares units just above and just below a threshold that determines treatment, such as a test score for a scholarship, assuming nothing else jumps at the threshold. Synthetic control builds a weighted combination of untreated units whose past matches the treated unit, and reads the effect from the post-treatment gap. All three can be fitted as Bayesian models with the tools of this course: multilevel models (Chapter 6), time-series models (Chapter 8) and Gaussian processes.

13.10 From causal inference to causal AI

Everything in this chapter was done with a few thousand rows and a handful of variables. Modern AI systems are trained on billions of examples with billions of parameters, and they predict superbly within the distribution they were trained on. Yet the failures that matter most, a detector that learned the scraper (Chapter 12), a clinical model that learned which hospital took the X-ray, a language model that answers "what if" questions by reciting associations, are failures on the second and third rungs of the ladder. The lesson of this chapter is that climbing the ladder takes three things that more data from the same source cannot supply:

1. Causal assumptions, written down as a graph or as structural equations, which say what may be adjusted for and what may not.
2. Interventional variation, from experiments, instruments, natural experiments or data from several environments in which the mechanisms are the same but the associations differ.
3. Honest uncertainty, which widens where the data cannot speak, as the overlap check and the sensitivity analysis did here.

Chapter 14 takes these one at a time into deep learning. It fits Bayesian neural networks as the structural equations of this section and asks whether their uncertainty widens where it should; it looks at transformers as machines that can learn Bayesian inference, and causal inference, from simulated worlds; and it asks what it would take for an AI system to learn causes instead of shortcuts.

13.11 Exercises

Each exercise has a hint and an answer button.

Further reading and references

1. Card, D. (1995). Using geographic variation in college proximity to estimate the return to schooling. In L. N. Christofides, E. K. Grant and R. Swidinsky (eds), Aspects of Labour Market Behaviour: Essays in Honour of John Vanderkamp, 201–222. University of Toronto Press.
2. Cinelli, C., Forney, A. and Pearl, J. (2024). A crash course in good and bad controls. Sociological Methods & Research, 53(3), 1071–1104.
3. Cinelli, C. and Hazlett, C. (2020). Making sense of sensitivity: extending omitted variable bias. Journal of the Royal Statistical Society, Series B, 82(1), 39–67.
4. Dehejia, R. H. and Wahba, S. (1999). Causal effects in nonexperimental studies: reevaluating the evaluation of training programs. Journal of the American Statistical Association, 94(448), 1053–1062.
5. Fisher, R. A. (1958). Cancer and smoking. Nature, 182, 596.
6. Hill, J. L. (2011). Bayesian nonparametric modeling for causal inference. Journal of Computational and Graphical Statistics, 20(1), 217–240.
7. Holland, P. W. (1986). Statistics and causal inference. Journal of the American Statistical Association, 81(396), 945–960.
8. Imbens, G. W. and Angrist, J. D. (1994). Identification and estimation of local average treatment effects. Econometrica, 62(2), 467–475.
9. LaLonde, R. J. (1986). Evaluating the econometric evaluations of training programs with experimental data. American Economic Review, 76(4), 604–620.
10. Louizos, C., Shalit, U., Mooij, J., Sontag, D., Zemel, R. and Welling, M. (2017). Causal effect inference with deep latent-variable models. Advances in Neural Information Processing Systems, 30.
11. Lundberg, I., Johnson, R. and Stewart, B. M. (2021). What is your estimand? Defining the target quantity connects statistical evidence to theory. American Sociological Review, 86(3), 532–565.
12. Neyman, J. (1923/1990). On the application of probability theory to agricultural experiments. Essay on principles. Section 9. Translated by D. M. Dabrowska and T. P. Speed. Statistical Science, 5(4), 465–472.
13. Pearl, J. (1993). Comment: graphical models, causality and intervention. Statistical Science, 8(3), 266–269.
14. Pearl, J. (1995). Causal diagrams for empirical research. Biometrika, 82(4), 669–688.
15. Pearl, J. (2001). Direct and indirect effects. Proceedings of the Seventeenth Conference on Uncertainty in Artificial Intelligence, 411–420.
16. Pearl, J. (2011). Invited commentary: understanding bias amplification. American Journal of Epidemiology, 174(11), 1223–1227.
17. Pearl, J. and Mackenzie, D. (2018). The Book of Why: The New Science of Cause and Effect. Basic Books.
18. Robins, J. (1986). A new approach to causal inference in mortality studies with a sustained exposure period: application to control of the healthy worker survivor effect. Mathematical Modelling, 7(9–12), 1393–1512.
19. Robins, J. M. and Greenland, S. (1992). Identifiability and exchangeability for direct and indirect effects. Epidemiology, 3(2), 143–155.
20. Rosenbaum, P. R. and Rubin, D. B. (1983). The central role of the propensity score in observational studies for causal effects. Biometrika, 70(1), 41–55.
21. Rubin, D. B. (1974). Estimating causal effects of treatments in randomized and nonrandomized studies. Journal of Educational Psychology, 66(5), 688–701.
22. Stock, J. H. and Trebbi, F. (2003). Retrospectives: who invented instrumental variable regression? Journal of Economic Perspectives, 17(3), 177–194.
23. VanderWeele, T. J. and Ding, P. (2017). Sensitivity analysis in observational research: introducing the E-value. Annals of Internal Medicine, 167(4), 268–274.
24. Wright, P. G. (1928). The Tariff on Animal and Vegetable Oils. Macmillan.