Chapter 13
Causal inference: one framework

Lectures 26–30 · potential outcomes and structural models in one framework, 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, setting the two great notations for causation, Rubin's potential outcomes and Pearl's structural models, in one Bayesian framework; and it prepares the three chapters that follow, on experiments and decisions (Chapter 14), on flexible nonparametric models (Chapter 15) and on making artificial intelligence itself model causes rather than correlations (Chapter 16).

• Lecture 26 is about the framework. Potential outcomes and the fundamental problem of causal inference; Rubin's model, in which causal inference is the Bayesian imputation of missing potential outcomes; what the data can and cannot identify; structural causal models as generators of potential outcomes; and single-world intervention graphs, which unite the two notations.
• Lecture 27 is about the problem in observational data. A job-training programme whose observational evaluation gives the wrong sign, and the back-door criterion.
• Lecture 28 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 29 is about unmeasured confounding. Instrumental variables with Card's college-proximity data; mediation; the front-door criterion.
• Lecture 30 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 — One framework for causation

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 Rubin's model: causal inference as missing data

Write every unit's potential outcomes in one table, a row per unit and a column for each treatment, next to the unit's covariates $X_i$. Donald Rubin called this the science table, because it holds everything there is to know about the effect of the treatment on these units. An experiment or an observational study reveals exactly half of it: unit $i$ has treatment $T_i$ and we observe $Y_i^{\text{obs}}=Y_i(T_i)$, while the other entry, $Y_i^{\text{mis}}=Y_i(1-T_i)$, stays blank.

Seen this way, causal inference is a missing-data problem, and Chapter 7 has already told us how a Bayesian solves those: by the posterior predictive distribution of whatever is missing (Rubin, 1978). The full model has two parts. The science, $p(X,Y(0),Y(1)\mid \theta )$, describes the units and their potential outcomes. The assignment mechanism, $p(T\mid X,Y(0),Y(1),\varphi )$, describes how units came to be treated. The posterior of the missing potential outcomes is

and any causal estimand, a function of the completed science table, has its posterior induced by this one. The assignment mechanism appears in the formula, and in general it matters: if the workers who would gain most are the ones who enrol, the treated units' outcomes tell us about a special subset of $Y(1)$.

This is the condition of Chapter 7 in different clothes. There, missing values were missing at random if the chance of being missing depended only on what was observed; here, each potential outcome is missing at random given $X$ if assignment depends only on $X$. Randomisation makes assignment ignorable by design, whatever the covariates, and that is the precise sense in which it licenses causal conclusions. In an observational study, ignorability given the measured covariates is an assumption, and the rest of this chapter is about when it is credible, what to do when it is not, and how to tell.

Two more assumptions are built into the science table itself, and Rubin (1980) gathered them under the name SUTVA, the stable unit treatment value assumption. First, no interference: unit $i$'s outcome depends only on unit $i$'s treatment, not on anyone else's. Vaccination (my protection depends on whether you are vaccinated), classroom interventions and online marketplaces all break it, and each needs a bigger table indexed by everybody's treatments. Second, no hidden versions: "training" means one thing, so there is one $Y_i(1)$ per unit. If some trainees got a good instructor and some a bad one, there are two treatments, not one.

The table also makes clear that "the effect" is not one quantity. The finite-sample average effect is the average over the units we have,

a property of these 445 men. The super-population average effect, ${\tau }_{\text{sp}}=𝔼[Y(1)-Y(0)]$, is a property of the population they were drawn from. Both are legitimate; they answer different questions and have different uncertainty, because ${\tau }_{\text{fs}}$ has half its terms already observed. And there is a quantity that the data say nothing about at all. Each unit reveals one potential outcome, so the data inform the two margins, the distributions of $Y(0)$ and of $Y(1)$, but never their joint distribution: we never see the same unit under both treatments, so we cannot learn how $Y(0)$ and $Y(1)$ are related within a person. The super-population average depends only on the margins, because $𝔼[Y(1)-Y(0)]=𝔼[Y(1)]-𝔼[Y(0)]$. The finite-sample average, the variance of individual effects and the share of people who benefit all depend on the joint distribution, and so on an assumption the data cannot check.

13.3 Imputing what was never seen

We carry out Rubin's programme on the randomised part of the job-training study described in Section 13.6: 185 men assigned to subsidised work and 260 controls, with earnings in 1978 as the outcome. Randomisation makes assignment ignorable, so we need only a model for the science. Earnings have a spike at zero (a quarter of the trainees and a third of the controls earned nothing in 1978) and a long right tail, so for each arm we use a hurdle lognormal: with probability $h_t$ a man earns nothing, and otherwise his earnings are lognormal with parameters ${\mu }_t$ and ${\sigma }_t$.

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)          # 1978 earnings, $1,000s
nsw |> summarise(n = n(), mean_earnings = mean(re78), share_zero = mean(re78 == 0), .by = treat)
m_rubin <- brm(bf(re78 ~ treat, hu ~ treat, sigma ~ treat), data = nsw, family = hurdle_lognormal(),
  prior = c(prior(normal(1.5, 1), class = Intercept), prior(normal(0, 1), class = b),
            prior(normal(0, 1.5), class = Intercept, dpar = hu), prior(normal(0, 1), class = b, dpar = hu),
            prior(normal(0, 1), class = Intercept, dpar = sigma), prior(normal(0, 0.5), class = b, dpar = sigma)),
  seed = 13, refresh = 0, file = "fits/ch13_rubin")
d <- as_draws_df(m_rubin)
arm <- tibble(mu0 = d$b_Intercept, mu1 = d$b_Intercept + d$b_treat,
              s0 = exp(d$b_sigma_Intercept), s1 = exp(d$b_sigma_Intercept + d$b_sigma_treat),
              h0 = plogis(d$b_hu_Intercept), h1 = plogis(d$b_hu_Intercept + d$b_hu_treat))
tau_sp <- with(arm, (1 - h1) * exp(mu1 + s1^2 / 2) - (1 - h0) * exp(mu0 + s0^2 / 2))   # E[Y(1)] - E[Y(0)]
c(mean = mean(tau_sp), quantile(tau_sp, c(0.05, 0.95)), prob_positive = mean(tau_sp > 0)) |> round(3)

# A tibble: 2 x 4
  treat     n mean_earnings share_zero
  <dbl> <int>         <dbl>      <dbl>
1     1   185          6.35      0.243
2     0   260          4.55      0.354
         mean            5%           95% prob_positive 
        1.686         0.094         3.380         0.960 

The trainees earned $6,350 on average in 1978 and the controls $4,550, and fewer of the trainees earned nothing (24% against 35%). The posterior mean of the super-population effect, $𝔼[Y(1)]-𝔼[Y(0)]$, is $1,690, with a 90% interval from $90 to $3,380 and a posterior probability of 0.96 that the effect is positive. The model is fitted to each arm separately, so this is a comparison of two fitted distributions, and the difference in sample means ($1,794) is close to its centre.

For the finite-sample effect we must fill in the blanks. For each posterior draw, every trainee needs a value of $Y(0)$ and every control a value of $Y(1)$, drawn from the fitted distribution of the missing arm. How the imputed value relates to the man's observed one is exactly the joint distribution the data cannot see. We make the assumption explicit with the Gaussian copula of Chapter 11: transform the observed outcome to a normal score $z$ through its own arm's distribution, draw the missing outcome's score as $\rho z+\sqrt{1-{\rho }^2}\epsilon$, and transform back through the other arm's distribution. With $\rho =0$ the two potential outcomes are independent given the parameters; with $\rho$ near 1 a man who earned much in one world earns much in the other, the assumption of rank preservation.

F_hl <- function(y, h, mu, s) ifelse(y <= 0, runif(length(y), 0, h), h + (1 - h) * pnorm((log(y) - mu) / s))
Q_hl <- function(u, h, mu, s) ifelse(u <= h, 0, exp(mu + s * qnorm(pmin((u - h) / (1 - h), 1 - 1e-12))))
impute_table <- function(rho, draws = 1000) {
  set.seed(13); idx <- round(seq(1, nrow(arm), length.out = draws)); t1 <- nsw$treat == 1; y <- nsw$re78
  map_dfr(idx, function(k) { a <- arm[k, ]
    u_obs <- ifelse(t1, F_hl(y, a$h1, a$mu1, a$s1), F_hl(y, a$h0, a$mu0, a$s0))          # observed outcome -> uniform
    z_mis <- rho * qnorm(u_obs) + sqrt(1 - rho^2) * rnorm(length(y))                     # Gaussian copula with correlation rho
    y_mis <- ifelse(t1, Q_hl(pnorm(z_mis), a$h0, a$mu0, a$s0), Q_hl(pnorm(z_mis), a$h1, a$mu1, a$s1))
    y1 <- ifelse(t1, y, y_mis); y0 <- ifelse(t1, y_mis, y)
    tibble(tau_fs = mean(y1 - y0), share_gain = mean(y1 > y0), share_lose = mean(y1 < y0), sd_effect = sd(y1 - y0)) })
}
imputed <- map_dfr(c(0, 0.5, 0.9), ~ impute_table(.x) |> mutate(rho = .x))
imputed |> summarise(fs_mean = mean(tau_fs), fs_q05 = quantile(tau_fs, 0.05), fs_q95 = quantile(tau_fs, 0.95),
                     share_gain = mean(share_gain), share_lose = mean(share_lose), sd_effect = mean(sd_effect), .by = rho) |>
  mutate(across(-rho, ~ round(.x, 3)))

# A tibble: 3 x 7
    rho fs_mean fs_q05 fs_q95 share_gain share_lose sd_effect
  <dbl>   <dbl>  <dbl>  <dbl>      <dbl>      <dbl>     <dbl>
1   0      1.88  0.69    3.27      0.519      0.395     12.9 
2   0.5    1.80  0.499   3.28      0.508      0.342     10.3 
3   0.9    1.58  0.242   2.99      0.541      0.236      4.99

Three things stand out. First, the finite-sample effect is estimated more precisely than the super-population effect: at $\rho =0$ its 90% interval runs from $690 to $3,270, against $90 to $3,380, because half of the 445 terms in its average are observed rather than modelled. Second, its centre moves a little with $\rho$, from $1,880 at $\rho =0$ to $1,580 at $\rho =0.9$: with a high correlation each man's imputed outcome is tied to his rank in his own arm, so the answer also inherits how well the hurdle lognormal fits the tails of each arm. Third, and most important, the quantities that describe individuals change a great deal. If the two potential outcomes are independent, 52% of the men gain from training and 40% lose, and individual effects have a standard deviation of $12,900; with $\rho =0.9$, 24% lose and the standard deviation is $4,990. (The rest earn nothing either way.) The data are identical in every row of the table; only an assumption the data cannot check has changed.

bind_rows(tibble(estimand = "super-population", tau = tau_sp),
          imputed |> transmute(estimand = paste0("finite sample, \u03c1 = ", rho), tau = tau_fs)) |>
  ggplot(aes(tau, colour = estimand)) + geom_density(linewidth = 0.8, adjust = 1.2) +
  geom_vline(xintercept = 0, colour = harbour[["grid"]]) +
  scale_colour_manual(values = unname(harbour[c("ink", "sky", "moss", "rose")])) +
  labs(x = "average effect on 1978 earnings ($1,000s)", y = NULL, colour = NULL) +
  guides(colour = guide_legend(nrow = 2)) + theme(axis.text.y = element_blank())

All four posteriors agree that training raised earnings. The super-population posterior is the widest, since it generalises from 445 men to the population they came from, and the finite-sample posteriors are narrower. The share of men who were harmed by training, which a policy-maker might care about most, is not identified at all: it ranged from 24% to 40% across the assumptions in the table.

Move the slider. The average effect stays put, because it is a difference of two means; the spread of individual effects collapses as $\rho$ rises, and the share of workers who would be harmed by training falls from about a third at $\rho =0$ to none at $\rho =1$. A policy-maker who asks "how many people does this programme hurt?" is asking a question that no experiment can answer without an assumption about $\rho$, and an honest analysis reports the answer as a function of that assumption.

13.4 Structural causal models generate potential outcomes

Pearl's framework starts from mechanisms rather than tables. A structural causal model writes each variable as a function of its direct causes and an exogenous disturbance, and draws the graph of who-listens-to-whom.

Each value $u$ of the exogenous variables is a unit: it fixes everything about that individual that the model does not explain. Solving the submodels $M_0$ and $M_1$ for the same $u$ gives $Y_0(u)$ and $Y_1(u)$, the unit's row of the science table. So an SCM generates potential outcomes, and Rubin's assumptions become consequences of the model. Consistency, $Y=Y_T$, holds because if $T(u)=t$ then setting $T$ to $t$ changes nothing. No interference holds because each unit has its own $u$ and no equation takes another unit's treatment as an input. And ignorability can be read from the graph.

The SCM also answers the question Section 13.3 had to leave open. In an SCM, a unit's two potential outcomes are computed from the same $u$, so the model determines their joint distribution. If earnings are $Y=f(T,U_Y)$ with $f$ increasing in $U_Y$, a man with a high $U_Y$ earns more in both worlds, and the two potential outcomes are perfectly rank-correlated: the $\rho =1$ case. Pearl's third rung, the individual counterfactual ("what would this trainee have earned without training?"), is answered by abduction, action and prediction: infer $u$ from what was observed about the unit, set $T$ to the other value, and recompute $Y$. With a probabilistic model and a prior on $u$, that is a posterior predictive imputation, Rubin's procedure, with the cross-world correlation fixed by the structural equations rather than left as a sensitivity parameter. The two frameworks agree on everything the data can identify; they differ in where they put the assumptions the data cannot check.

13.5 One framework, two notations

The graph and the potential outcomes meet in a single picture. Thomas Richardson and James Robins (2013) proposed the single-world intervention graph (SWIG). To represent the intervention $\text{do}(T=t)$, split the node $T$ in two: the random variable $T$ keeps the arrows coming into it, and a fixed value $t$ sends the arrows going out. Every descendant of the split node is relabelled as a potential outcome, $Y$ becoming $Y(t)$.

The SWIG is an ordinary graph, so d-separation applies to it, and it contains $T$ and $Y(t)$ as nodes. Here the only path between them is $T\leftarrow C\to Y(t)$, which conditioning on $C$ blocks: the graph displays the ignorability statement $Y(t)⟂T\mid C$ directly, rather than deriving it. If there were an unmeasured common cause of $T$ and $Y$, the SWIG would show an open path and ignorability would fail. SWIGs also show what cannot be represented in a single world. The natural direct effect of Section 13.11 involves $Y(t,M(t^{\prime }))$, the outcome under one treatment with the mediator at the value it would take under another. No single intervention produces both values, so no SWIG contains them, and the independence that identifies natural effects is a cross-world assumption that no experiment can test. It is the same lesson as the correlation $\rho$ of Section 13.3, at a larger scale.

The table is a dictionary, not a scorecard; each tradition is stronger where the other is weaker. The potential-outcome tradition keeps attention on the design: how units came to be treated, whether there is overlap, which estimand a policy needs, and how to analyse experiments and their finite populations. The graphical tradition keeps attention on assumptions: it says exactly which variables to adjust for and which to leave alone (Section 13.8), finds identification where no ignorability statement is natural (the front door, Section 13.11), derives testable implications (Chapter 2) and transports results between populations. Guido Imbens (2020) reviews both from an economist's point of view and argues that each community should borrow the other's strengths. This course uses graphs to choose models and potential outcomes to state and compute estimands, and the rest of this chapter moves between them freely.

Lecture 27 — The problem in observational data

13.6 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.7 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 28 — Adjustment

13.8 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.10. 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.9 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 29 — When the back door is closed

13.10 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.11 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 30 — Flexibility and humility

13.12 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.13 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. Chapter 15 explains why flexible outcome models work for causal effects, when they fail (regularisation-induced confounding and poor overlap) and how Bayesian causal forests repair the first failure.

13.13 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.

f_pl <- re75 ~ treat * (age + educ + black + hisp + married + nodegree + re74)
m_pl <- brm(f_pl, data = lalonde, seed = 13, refresh = 0, file = "fits/ch13_placebo",
            prior = c(prior(normal(3, 10), class = Intercept), prior(normal(0, 5), class = b), prior(exponential(0.2), class = sigma)))
tr <- lalonde$treat == 1
att_pl <- rowMeans(posterior_epred(m_pl, newdata = mutate(lalonde, treat = 1)[tr, ]) -
                   posterior_epred(m_pl, newdata = mutate(lalonde, treat = 0)[tr, ]))
exp_pl <- with(nsw_mixtape, mean(re75[treat == 1]) - mean(re75[treat == 0])) / 1000
c(observational_placebo = mean(att_pl), q05 = unname(quantile(att_pl, 0.05)), q95 = unname(quantile(att_pl, 0.95)),
  experimental_placebo = exp_pl) |> round(3)

observational_placebo                   q05                   q95 
                0.131                -0.380                 0.663 
 experimental_placebo 
                0.265 

In our data the placebo test passes: the estimated "effect" of training on 1975 earnings is $130, with a 90% interval from −$380 to $660, against $265 in the experiment, where randomisation guarantees no real effect. This differs from Imbens and Xu's result for an instructive reason. They use the full PSID comparison group of 2,490 men; this chapter uses the 429-man comparison group distributed with the MatchIt package, a subset that already resembles the trainees much more closely. Selecting a comparison group is itself part of the design, and a placebo test that passes in one sample can fail in another. Passing is necessary for credibility, not sufficient: 1974 earnings are in the model, so 1975 earnings are predicted well, but the real question is whether anything that shaped 1978 earnings and enrolment is missing, which no placebo on an earlier outcome can settle.

13.14 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.

The rest of the course takes these one at a time. Chapter 14 turns to the cleanest source of interventional variation, the randomised experiment, and to the decisions experiments are run to inform. Chapter 15 replaces the fixed functional forms of this chapter with priors on whole distributions and functions, and asks when their flexibility helps and when it only hides an assumption. Chapter 16 takes all three requirements 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.15 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. and Hazlett, C. (2020). Making sense of sensitivity: extending omitted variable bias. Journal of the Royal Statistical Society, Series B, 82(1), 39–67.
3. Cinelli, C., Forney, A. and Pearl, J. (2024). A crash course in good and bad controls. Sociological Methods & Research, 53(3), 1071–1104.
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. (2020). Potential outcome and directed acyclic graph approaches to causality: relevance for empirical practice in economics. Journal of Economic Literature, 58(4), 1129–1179.
9. Imbens, G. W. and Angrist, J. D. (1994). Identification and estimation of local average treatment effects. Econometrica, 62(2), 467–475.
10. Imbens, G. W. and Xu, Y. (2025). Comparing experimental and nonexperimental methods: what lessons have we learned four decades after LaLonde (1986)? Journal of Economic Perspectives, 39(4), 173–202.
11. LaLonde, R. J. (1986). Evaluating the econometric evaluations of training programs with experimental data. American Economic Review, 76(4), 604–620.
12. 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.
13. 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.
14. Pearl, J. (1993). Comment: graphical models, causality and intervention. Statistical Science, 8(3), 266–269.
15. Pearl, J. (1995). Causal diagrams for empirical research. Biometrika, 82(4), 669–688.
16. Pearl, J. (2001). Direct and indirect effects. Proceedings of the Seventeenth Conference on Uncertainty in Artificial Intelligence, 411–420.
17. Pearl, J. (2009). Causality: Models, Reasoning, and Inference, 2nd edition. Cambridge University Press.
18. Pearl, J. (2011). Invited commentary: understanding bias amplification. American Journal of Epidemiology, 174(11), 1223–1227.
19. Pearl, J. and Mackenzie, D. (2018). The Book of Why: The New Science of Cause and Effect. Basic Books.
20. Richardson, T. S. and Robins, J. M. (2013). Single world intervention graphs (SWIGs): a unification of the counterfactual and graphical approaches to causality. Working Paper 128, Center for Statistics and the Social Sciences, University of Washington.
21. 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.
22. Robins, J. M. and Greenland, S. (1992). Identifiability and exchangeability for direct and indirect effects. Epidemiology, 3(2), 143–155.
23. 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.
24. Rubin, D. B. (1974). Estimating causal effects of treatments in randomized and nonrandomized studies. Journal of Educational Psychology, 66(5), 688–701.
25. Rubin, D. B. (1978). Bayesian inference for causal effects: the role of randomization. Annals of Statistics, 6(1), 34–58.
26. Rubin, D. B. (1980). Comment on "Randomization analysis of experimental data: the Fisher randomization test" by D. Basu. Journal of the American Statistical Association, 75(371), 591–593.
27. Splawa-Neyman, J. (1990). On the application of probability theory to agricultural experiments. Essay on principles. Section 9 (translated by D. M. Dabrowska and T. P. Speed; original work 1923). Statistical Science, 5(4), 465–472.
28. Stock, J. H. and Trebbi, F. (2003). Retrospectives: who invented instrumental variable regression? Journal of Economic Perspectives, 17(3), 177–194.
29. VanderWeele, T. J. and Ding, P. (2017). Sensitivity analysis in observational research: introducing the E-value. Annals of Internal Medicine, 167(4), 268–274.
30. Wright, P. G. (1928). The Tariff on Animal and Vegetable Oils. Macmillan.