Chapter 1
The causal–Bayesian workflow

Lectures 1–2 · small worlds, graphs, priors and posteriors: the loop that every later chapter runs

In November 2020, a pharmaceutical company announced that its new COVID-19 vaccine was about 95% effective. The number came from 170 sick people. Eight of them had received the vaccine and 162 had received a placebo injection. From those two counts, and from a careful plan written down months earlier, the company's statisticians produced a statement that billions of people would act on: given these data, the probability that the vaccine's efficacy exceeds 30% is essentially 1.

That statement is Bayesian. It is a probability about something we cannot see (the vaccine's true efficacy), computed from something we can see (the counts), using rules stated in advance. You will reproduce it, exactly, in Section 1.2. But the more important lesson of this chapter is what stands behind such a number. Two kinds of assumption make it meaningful. One is a set of claims about what causes what: that the vaccine was assigned by a coin toss, so nothing else differs systematically between the two groups. The other is a set of claims about how the data were generated and what was plausible beforehand. Make the first kind of assumption wrongly and the most precise posterior in the world answers a question nobody asked. Make the second kind wrongly and the posterior answers the right question badly.

This course is about doing applied Bayesian research well enough to publish it and to rely on it. This chapter sets out the way of working that every later chapter follows. It has two halves that must be learned together:

• a causal half, in which we make our assumptions about the world explicit as a graph, and use the graph to decide which model to fit;
• a Bayesian half, in which we turn those assumptions, plus data, into a posterior distribution, and then check whether the result could plausibly have produced the data we saw.

Neither half is much use alone. A graph with no probability model cannot be fitted. A probability model with no graph can be fitted, but it cannot be interpreted. Here is the route through the chapter.

• Sections 1.1 to 1.3 set up the Bayesian half: what a model is, how a posterior is computed, and how to read one. The running example is the 2020 vaccine trial.
• Sections 1.4 to 1.6 set up the causal half: directed acyclic graphs, the four structures every graph is built from, and the back-door criterion that tells you which variables belong in a model.
• Sections 1.7 and 1.8 put the two halves together in one complete analysis of Antarctic penguins: simulate the priors, fit the model, check it, find what is wrong, and fix it.
• Section 1.9 draws the whole thing as a loop, and shows what the first step looks like for a real research project with ten variables.

Lecture 1 — Small worlds and the Bayesian update

1.1 Small worlds and large worlds

A statistical model is a small world: a self-contained, fully specified cartoon of reality inside which every event has a probability and every question has an answer. The large world is the messy reality the data actually came from. All inference happens in the small world. Its conclusions are only as trustworthy as the resemblance between the two.

That is why a model can be computed perfectly and still be wrong, and why most of the intellectual work in applied Bayesian analysis is not computation. It is construction (deciding what to put into the small world) and criticism (checking afterwards whether the small world could plausibly have produced what we saw). Computers have made the middle step, computing the posterior, almost free. They have done nothing for the first and last steps, which is where analyses succeed or fail.

Every model in this course has three layers, and we will write each one down explicitly.

1. A causal layer: a graph saying what influences what. It decides which variables enter the model and what their coefficients mean.
2. A stochastic layer: a likelihood, saying how the data are generated if we knew the unknowns, and a prior, saying which values of the unknowns are plausible before we see the data.
3. A computational layer: the algorithm that turns the first two layers into a posterior distribution.

Most failures of applied analysis are failures of layer one disguised as questions about layer three. "Should I run the sampler longer?" is rarely the question that matters. "Should this variable be in the model at all?" often is.

The posterior is the entire output of a Bayesian analysis. It is not an estimate; it is a distribution, a weighting over every version of the small world that is consistent with the data and the prior. Point estimates and intervals are summaries we extract from it at the end, and we will resist collapsing it to a single number until the last possible moment.

1.2 A first Bayesian analysis: did the vaccine work?

The data and the question

The trial of the Pfizer–BioNTech vaccine (Polack et al., 2020) randomised about 43,000 adults to two injections of either the vaccine or a saltwater placebo, in equal numbers. Then everyone went back to their lives, and the investigators waited for people to catch COVID-19. By the final analysis, 170 participants had developed symptomatic, laboratory-confirmed COVID-19 at least a week after their second dose. Eight were in the vaccine group. 162 were in the placebo group.

The question is: how effective is the vaccine? Vaccine efficacy is defined as the proportional reduction in the rate of disease:

A VE of 0 means the vaccine does nothing. A VE of 1 means it prevents every case. A negative VE means it makes disease more likely. Regulators had set the bar in advance: a vaccine would be considered effective if VE was convincingly above 30%.

The causal sketch

Before any probability, a graph. The question is causal: what would the rate of disease be if a person received the vaccine, compared with if they did not? Many things affect whether a person catches COVID-19: their age, their job, how many people they meet, how careful they are. Call all of that $U$. In an observational study, $U$ would also affect whether people chose to be vaccinated (cautious people get vaccinated and avoid crowds), and the comparison would be hopelessly confounded. Randomisation deletes that arrow.

That one missing arrow, from $U$ into $V$, is the entire reason the trial's numbers can be read causally. We will come back to this in Section 1.6. For now, note that the graph has already done its work: it has told us that a simple model comparing the two arms is enough.

A clever likelihood

The trial's statisticians used a neat trick. Because the two arms were the same size and followed for almost exactly the same length of time, we can condition on the total number of cases, 170, and ask a simpler question: of the people who got sick, what fraction were in the vaccine arm? Call that fraction $\theta$.

If the vaccine did nothing, a sick person would be equally likely to come from either arm, and $\theta$ would be 1/2. If the vaccine worked perfectly, nobody in the vaccine arm would get sick, and $\theta$ would be 0. In general, with equal-sized arms, the case split and the efficacy are linked by

Given $\theta$, each of the 170 cases is an independent "draw" that lands in the vaccine arm with probability $\theta$. So the number of vaccine-arm cases follows a binomial distribution:

A prior written down in advance

Here is the part that made headlines among statisticians. The trial's protocol, published before the data came in, specified the prior:

The strange first parameter was chosen for a reason. The mean of a Beta$(a,b)$ distribution is $a/(a+b)=0.700102/1.700102=0.4118$, and plugging $\theta =0.4118$ into the formula above gives VE = 30%. So the prior is centred on the regulators' minimum acceptable efficacy, and the total "weight" $a+b\approx 1.7$ makes it worth less than two observations. It is deliberately modest: the prior is centred on "just barely good enough", and is weak enough that the data will decide.

The update

The Beta distribution and the binomial likelihood fit together perfectly. If the prior is Beta$(a,b)$ and we observe $y$ successes in $n$ trials, the posterior is again a Beta distribution:

So the vaccine posterior is $\text{Beta}(0.700102+8,1+162)=\text{Beta}(8.700102,163)$. Every question about efficacy can now be answered by transforming this distribution. The easiest way to do that is by simulation: draw many values of $\theta$, convert each one to VE, and summarise the draws. Here is the first R code of the course. It starts by loading the packages we will use throughout and a plotting theme that matches these notes.

library(tidyverse)   # data wrangling and ggplot2
library(brms)        # Bayesian regression models, fitted with Stan
library(dagitty)     # causal graphs

# colours and a ggplot2 theme matching these notes
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())

Now the analysis. Read it line by line: it is the whole of Bayesian inference in six lines.

a_prior <- 0.700102; b_prior <- 1          # the protocol's prior on theta
y <- 8; n <- 170                           # cases in the vaccine arm, total cases
a_post <- a_prior + y; b_post <- b_prior + n - y

ve <- function(theta) (1 - 2 * theta) / (1 - theta)   # theta -> efficacy

set.seed(2020)
ve_draws <- ve(rbeta(100000, a_post, b_post))    # 100,000 posterior draws of VE
round(quantile(ve_draws, c(0.025, 0.5, 0.975)), 3)
mean(ve_draws > 0.30)                            # Pr(VE > 30% | data)
pbeta(0.7 / 1.7, a_post, b_post, lower.tail = FALSE)   # exact Pr(VE <= 30% | data)

 2.5%   50% 97.5% 
0.903 0.949 0.976 
[1] 1
[1] 1.96e-28

The posterior median efficacy is 94.9%, and the central 95% of the posterior lies between 90.3% and 97.6%. The interval is exactly the 95% credible interval reported to the US Food and Drug Administration in December 2020 (FDA, 2020): 95.0% (90.3%, 97.6%). The point estimate of 95.0% there is not a posterior median but the observed efficacy, one minus the ratio of the two arms' case rates per person-year, $1-(8/2.214)/(162/2.222)$; the exact posterior median is 94.9%. None of our 100,000 draws fell below 30%. The exact calculation in the last line shows why: the posterior probability that VE is 30% or less is about $2\times {10}^{-28}$. The protocol required that probability to be below 1.4%.

Picture the prior and the posterior on the scale people care about, the VE scale.

draws <- tibble(
  prior     = ve(rbeta(100000, a_prior, b_prior)),
  posterior = ve_draws
) |>
  pivot_longer(everything(), names_to = "distribution", values_to = "VE") |>
  filter(VE > -1)

draws$distribution <- factor(draws$distribution, levels = c("prior", "posterior"))

ggplot(draws, aes(VE, fill = distribution, colour = distribution)) +
  geom_density(alpha = 0.35, adjust = 1.2, linewidth = 0.8) +
  geom_vline(xintercept = 0.3, linetype = "dotted", colour = harbour[["rose"]]) +
  facet_wrap(~ distribution, ncol = 1, scales = "free_y") +
  scale_fill_manual(values = c(prior = harbour[["amber"]], posterior = harbour[["sky"]])) +
  scale_colour_manual(values = c(prior = harbour[["amber"]], posterior = harbour[["sky"]])) +
  scale_x_continuous(labels = scales::percent) +
  labs(x = "vaccine efficacy", y = "density") +
  theme(legend.position = "none")

Notice the shape of the prior on this scale. Efficacy can be at most 100%, but it is unbounded below: a vaccine that tripled your risk would have VE = −200%. The protocol's prior gives such catastrophes some weight (its 2.5th percentile is around −2,600%), but it also gives just over half its weight to efficacy above 40%. A prior that looks modest on one scale can look strange on another, and that is a lesson we will meet again in Chapter 2.

Now experiment with the update yourself. The widget below performs exactly the calculation above, but lets you change the data and the prior.

Three things are worth trying. First, set both counts to zero: the posterior is the prior, because there are no data. Second, replay the trial and watch how quickly the posterior abandons the prior: after twenty or thirty cases, the choice of prior barely matters. Third, switch to the sceptical prior, which is centred on VE = 0 and worth 40 observations, and lower the case counts to something small like 2 and 15. With few data, the prior visibly holds the posterior back. Priors matter most when data are scarce, which is exactly when you should think hardest about them.

1.3 Reading a posterior

A posterior is a distribution, and a distribution answers many questions at once. Three summaries cover most needs.

• A central value, usually the posterior median or mean: a single "best guess". For the vaccine, the median VE is 94.9%.
• A credible interval: a range containing a stated share of the posterior. The 95% equal-tailed interval runs from the 2.5th to the 97.5th percentile of the draws: 90.3% to 97.6%.
• Posterior probabilities of statements: the share of draws for which a statement is true. Pr(VE > 30% | data) ≈ 1. Pr(VE > 95% | data) is about one half.

The phrase "given the model" carries a lot of weight. A credible interval is a statement in the small world. If the likelihood is wrong, or the prior is absurd, or (most often) the causal assumptions that make the parameter meaningful are false, the interval is a precise answer to the wrong question. That is not a weakness peculiar to Bayesian methods. A frequentist confidence interval rests on the same model; it just states its conclusion differently.

Two habits will serve you well. First, report the distribution, not just a point: a posterior median without an interval is a sentence with its verb missing. Second, prefer probabilities of substantive statements to tests against zero. "There is a 97% probability that the programme raised earnings by more than $500" is something a policymaker can use. "The effect is significantly different from zero" is not.

1.4 Drawing assumptions as a graph

The vaccine analysis was easy to interpret because of a coin toss. Most data are not like that. Most data come from observing the world as it is, where the things we want to compare differ in many ways at once. To reason about such data we need a precise language for causal assumptions. This course uses the directed acyclic graph, or DAG.

A DAG is a compact, checkable statement of a scientific theory. It is not a statistical model: it says nothing about the size or shape of any effect. It says only what can affect what. That turns out to be enough to decide which statistical model to fit, and that is the remarkable fact this course is built on.

An example: penguin bills

The palmerpenguins data (Horst, Hill and Gorman, 2020) record measurements of 344 penguins of three species (Adélie, Chinstrap and Gentoo) collected from 2007 to 2009 by Kristen Gorman at Palmer Station, Antarctica. Two of the measurements describe the bill: its length and its depth (top to bottom). Is a longer bill a shallower one?

library(palmerpenguins)
bills <- penguins |> drop_na(bill_length_mm, bill_depth_mm)

# one regression line through all 342 penguins
coef(lm(bill_depth_mm ~ bill_length_mm, data = bills))

# a separate line for each species
bills |>
  group_by(species) |>
  summarise(n = n(), slope = coef(lm(bill_depth_mm ~ bill_length_mm))[2])

   (Intercept) bill_length_mm 
      20.88547       -0.08502 
# A tibble: 3 x 3
  species       n slope
  <fct>     <int> <dbl>
1 Adelie      151 0.179
2 Chinstrap    68 0.222
3 Gentoo      123 0.205

Pooling all the penguins, each extra millimetre of bill length goes with a bill about 0.085 mm shallower. Within every one of the three species, each extra millimetre goes with a bill about 0.18 to 0.22 mm deeper. The sign flips. This is Simpson's paradox, and it is not a paradox at all once you draw the graph.

species_cols <- c(Adelie = harbour[["sky"]], Chinstrap = harbour[["amber"]],
                  Gentoo = harbour[["moss"]])
ggplot(bills, aes(bill_length_mm, bill_depth_mm)) +
  geom_point(aes(colour = species), alpha = 0.6, size = 1.6) +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE,
              colour = harbour[["muted"]], linetype = "dashed") +
  geom_smooth(aes(colour = species), method = "lm", formula = y ~ x, se = FALSE) +
  scale_colour_manual(values = species_cols) +
  labs(x = "bill length (mm)", y = "bill depth (mm)", colour = NULL)

Species affects both measurements. Gentoo penguins have long, shallow bills; Adélie penguins have short, deep ones. Within a species, bigger birds have bigger bills in every direction, so length and depth rise together. Here is that story as a graph:

There is no arrow between length and depth. Both associations, the negative pooled one and the positive within-species one, are produced by common causes. Pooling mixes the between-species pattern (driven by $S$) with the within-species pattern (driven by $B$). Separating by species removes the first and leaves the second. Neither line measures what would happen if you could stretch a penguin's bill.

Lecture 2 — Graphs, priors and the workflow

1.5 The four elementary structures

Every DAG, however large, is built from four local patterns. Knowing how each behaves is most of what causal reasoning requires in practice. Let $X$ and $Y$ be two variables whose association we are studying, and $Z$ a third.

The fork and the pipe behave identically in the data: in both, $X$ and $Y$ are associated, and the association disappears when you hold $Z$ fixed. The collider is the opposite: no association until you hold $Z$ fixed. The collider deserves special attention, because it is the structure that defeats intuition most reliably.

Suppose a city's restaurants survive if they are either genuinely good or cheaply located, and that food quality and rent are unrelated across all the restaurants that ever open. Among surviving restaurants, the only ones you can sample, a negative association appears: the expensive survivors must be good, and the mediocre survivors must be cheap. Survival is a collider, quality → survival ← cheap location, and sampling only survivors conditions on it. The association is created entirely by the selection. Explore all four structures in the widget.

The practical warning is blunt. "Control for more variables" is not a safe default. Conditioning on a variable is the right move for a fork and the wrong move for a collider, and only the graph, not the data, tells you which you are facing. The data look the same either way.

1.6 Paths, blocking and the back-door criterion

The four rules combine into one procedure for any graph. Between an exposure $X$ and an outcome $Y$ there may be many paths: sequences of arrows connecting them, followed in any direction. A path transmits association unless something blocks it.

Paths come in two kinds. A causal path follows the arrows forward from $X$ to $Y$; it carries the effect we want. A back-door path starts with an arrow pointing into $X$; it carries association that is not caused by $X$. The vaccine graph had no back-door paths because nothing pointed into $V$. That is precisely what randomisation buys.

Condition 1 rules out mediators and their descendants (Section 1.5). Condition 2 closes the confounding. You do not have to trace paths by hand: the dagitty package (Textor et al., 2016) reads a graph and returns the sufficient adjustment sets. Here it is on a small example. Suppose we want the effect of exercise on blood pressure, where age affects both, and both affect whether a person enrolled in the study through a health-screening programme.

g <- dagitty("dag {
  Age -> Exercise
  Age -> BP
  Exercise -> BP
  Exercise -> Enrolled
  BP -> Enrolled
}")
adjustmentSets(g, exposure = "Exercise", outcome = "BP")
isAdjustmentSet(g, c("Age", "Enrolled"), exposure = "Exercise", outcome = "BP")

{ Age }
[1] FALSE

The answer is to adjust for age and nothing else. Adding enrolment status, a collider, is declared invalid, and the data could never have told you that.

It is worth pausing on what has just happened, because it is the organising principle of the whole course. The back-door criterion is a statement about which variables enter the regression, and it was settled before any Bayesian machinery appeared. The causal layer chooses the model; the Bayesian layer fits it and quantifies uncertainty. Reversing the order, fitting first and interpreting causally later, is the most common way good computation produces a wrong conclusion. Collecting more data or running a longer chain cannot repair an adjustment set chosen by habit.

1.7 Priors are assumptions, so simulate them

Now we combine the two halves in a complete analysis. The question is practical. Weighing a penguin in the field means catching it and holding it on a scale in a bag, which is stressful for the bird and awkward for the researcher. Measuring a flipper with calipers is quick. How well can we predict a penguin's body mass from its flipper length?

First the graph. Flipper length does not cause body mass: stretching a penguin's flipper will not make it heavier. Both are consequences of the bird's overall size, which depends on its species and its sex.

For prediction we do not need an arrow from $F$ to $M$. Flipper length is useful because it is a descendant of size: it carries information about $B$, which drives $M$. The graph also warns us that species and sex will matter, and that a model ignoring them will mix several different relationships. We will start with the simplest model anyway, and let the checks tell us what is missing.

Measuring mass $M_i$ in kilograms and flipper length $F_i$ in millimetres, centred at 200 mm, the first model is a straight line with normal scatter:

Now we need priors for $\alpha$ (the expected mass of a penguin with 200 mm flippers), $\beta$ (the change in expected mass per millimetre of flipper) and $\sigma$ (the scatter). It is tempting to use "uninformative" priors, such as Normal(0, 10) for both. Before using any prior, we should look at what it says.

set.seed(11)
n_lines <- 50
prior_lines <- tibble(
  prior = rep(c("vague: Normal(0, 10) on both",
                "weakly informative: Normal(4, 1), Normal(0, 0.05)"), each = n_lines),
  alpha = c(rnorm(n_lines, 0, 10), rnorm(n_lines, 4, 1)),
  beta  = c(rnorm(n_lines, 0, 10), rnorm(n_lines, 0, 0.05))
)

ggplot(prior_lines) +
  geom_segment(aes(x = 170, xend = 235,
                   y = alpha + beta * (170 - 200), yend = alpha + beta * (235 - 200)),
               colour = harbour[["sky"]], alpha = 0.45) +
  geom_hline(yintercept = c(0, 22), linetype = "dashed", colour = harbour[["rose"]]) +
  facet_wrap(~ prior) +
  coord_cartesian(ylim = c(-15, 35)) +
  labs(x = "flipper length (mm)", y = "body mass (kg)")

# share of prior lines that imply an impossible penguin somewhere in 170-235 mm
prior_lines |>
  mutate(lo = pmin(alpha - 30 * beta, alpha + 35 * beta),
         hi = pmax(alpha - 30 * beta, alpha + 35 * beta)) |>
  group_by(prior) |>
  summarise(impossible = mean(lo < 0 | hi > 22))

# A tibble: 2 x 2
  prior                                             impossible
  <chr>                                                  <dbl>
1 vague: Normal(0, 10) on both                            1   
2 weakly informative: Normal(4, 1), Normal(0, 0.05)       0.02

Under the vague priors, every one of the 50 lines implies penguins that weigh less than nothing, or more than any bird alive, somewhere within the range of real flipper lengths. Under the weakly informative priors, one line in fifty does. The weakly informative priors were not chosen by looking at these data. They came from things anyone could look up: penguins of these species weigh a few kilograms (so $\alpha \sim \text{Normal}(4,1)$), and across the 60 mm range of flipper lengths the mass is unlikely to change by more than a few kilograms (so $\beta \sim \text{Normal}(0,0.05)$, allowing slopes up to about ±0.1 kg per mm). Notice that we did not even assume the slope is positive.

Experiment with the priors yourself.

1.8 Fit, diagnose, check, revise

Fit

With the graph drawn and the priors simulated, we fit the model. The brms package translates a model formula into a program in the probabilistic language Stan, compiles it, and draws samples from the posterior using Hamiltonian Monte Carlo (Chapter 5). Treat the sampler as a black box for now; we will open it in Chapter 5.

pen <- penguins |>
  drop_na(body_mass_g, flipper_length_mm, sex) |>
  mutate(mass = body_mass_g / 1000,              # kilograms
         flipper = flipper_length_mm - 200)      # mm, centred at 200

m1 <- brm(
  mass ~ 1 + flipper, data = pen, family = gaussian(),
  prior = c(prior(normal(4, 1), class = Intercept),
            prior(normal(0, 0.05), class = b),
            prior(exponential(1), class = sigma)),
  seed = 1, refresh = 0, file = "fits/ch1_m1"     # 'file' caches the fitted model
)
summary(m1)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: mass ~ 1 + flipper 
   Data: pen (Number of observations: 333) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     4.16      0.02     4.12     4.20 1.00     2006     1767
flipper       0.05      0.00     0.05     0.05 1.00     5705     3167

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.39      0.02     0.36     0.43 1.00     2157     2134

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Read the table from the bottom up. sigma, the typical deviation of a penguin from the line, is about 0.39 kg. The Intercept says a penguin with 200 mm flippers is expected to weigh 4.16 kg (95% interval 4.12 to 4.20). The flipper coefficient says expected mass rises about 0.05 kg, or 50 g, for each millimetre of flipper.

Diagnose

Before believing those numbers, check that the computation worked. Three columns of the summary do that job.

• Rhat compares the four independent chains with each other. If they have all explored the same distribution, Rhat is 1.00. Values above 1.01 mean the chains disagree, and the summary should not be trusted.
• Bulk_ESS and Tail_ESS are effective sample sizes: how many independent draws the correlated samples are worth, for the centre and for the tails of the distribution. A few hundred is enough for most summaries.
• Any warnings about divergent transitions (there were none) mean the sampler met geometry it could not handle (Chapter 5).

Here every Rhat is 1.00 and every ESS is in the thousands, so the computation is sound. An algorithm that fails these checks has not given you a slightly worse answer. It has given you no answer.

Check

The computation worked. Does the model? A posterior predictive check draws parameters from the posterior, simulates new datasets, and compares them with the data actually observed. If the model is a good small world, its simulated penguins should look like the real ones in every way we care about. Our graph told us that species and sex matter, so a natural check is: does the model reproduce the average mass of each species and sex?

We also fit the model the graph suggests, adding species and sex, so that the check can compare the two.

m2 <- brm(
  mass ~ 1 + flipper + species + sex, data = pen, family = gaussian(),
  prior = c(prior(normal(4, 1), class = Intercept),
            prior(normal(0, 1), class = b),
            prior(normal(0, 0.05), class = b, coef = flipper),
            prior(exponential(1), class = sigma)),
  seed = 1, refresh = 0, file = "fits/ch1_m2"
)
round(fixef(m2), 3)

                 Estimate Est.Error   Q2.5 Q97.5
Intercept           3.643     0.047  3.548 3.732
flipper             0.020     0.003  0.014 0.026
speciesChinstrap   -0.089     0.046 -0.178 0.001
speciesGentoo       0.829     0.084  0.667 0.997
sexmale             0.528     0.038  0.452 0.602

group_means <- function(fit) {
  yrep <- posterior_predict(fit, ndraws = 500)      # 500 simulated datasets
  grp  <- paste(pen$species, pen$sex)
  sims <- apply(yrep, 1, function(y) tapply(y, grp, mean))
  tibble(group = rownames(sims),
         lo = apply(sims, 1, quantile, 0.05),
         hi = apply(sims, 1, quantile, 0.95))
}
set.seed(3)
checks <- bind_rows("flipper only" = group_means(m1),
                    "flipper + species + sex" = group_means(m2), .id = "model")
observed <- pen |> group_by(group = paste(species, sex)) |> summarise(mass = mean(mass))

ggplot(checks, aes(y = group)) +
  geom_linerange(aes(xmin = lo, xmax = hi, colour = model), linewidth = 3,
                 position = position_dodge(width = 0.55)) +
  geom_point(data = observed, aes(x = mass), shape = 18, size = 3.6,
             colour = harbour[["ink"]]) +
  scale_colour_manual(values = c("flipper only" = harbour[["amber"]],
                                 "flipper + species + sex" = harbour[["sky"]])) +
  labs(x = "average body mass (kg)", y = NULL, colour = NULL)

# the same comparison as a table
checks |>
  left_join(observed, by = "group") |>
  mutate(inside = mass >= lo & mass <= hi) |>
  select(model, group, observed = mass, lo, hi, inside) |>
  arrange(group, model) |>
  print(n = 12, width = 82)

# A tibble: 12 x 6
   model                   group            observed    lo    hi inside
   <chr>                   <chr>               <dbl> <dbl> <dbl> <lgl> 
 1 flipper + species + sex Adelie female        3.37  3.32  3.47 TRUE  
 2 flipper only            Adelie female        3.37  3.47  3.64 FALSE 
 3 flipper + species + sex Adelie male          4.04  3.94  4.09 TRUE  
 4 flipper only            Adelie male          4.04  3.69  3.88 FALSE 
 5 flipper + species + sex Chinstrap female     3.53  3.29  3.50 FALSE 
 6 flipper only            Chinstrap female     3.53  3.63  3.88 FALSE 
 7 flipper + species + sex Chinstrap male       3.94  3.98  4.18 FALSE 
 8 flipper only            Chinstrap male       3.94  4.05  4.27 FALSE 
 9 flipper + species + sex Gentoo female        4.68  4.65  4.82 TRUE  
10 flipper only            Gentoo female        4.68  4.70  4.89 FALSE 
11 flipper + species + sex Gentoo male          5.48  5.36  5.52 TRUE  
12 flipper only            Gentoo male          5.48  5.13  5.34 FALSE 

The flipper-only model fails the check badly. It predicts that female Adélies average 3.47 to 3.64 kg (they actually average 3.37 kg), that male Adélies average 3.69 to 3.88 kg (they average 4.04 kg), and that male Gentoos average 5.13 to 5.34 kg (they average 5.48 kg). Look at the pattern: every female group is over-predicted, and the heavier Adélie and Gentoo males are under-predicted. The model is not merely noisy; it is systematically wrong in exactly the way the graph warned, because it has no way to know that males are heavier than females with the same flippers. The model with species and sex fits four of the six groups, and its residual misses are informative: it over-predicts Chinstrap males and under-predicts Chinstrap females, which suggests that the sex difference in mass is smaller in Chinstraps than in the other species. That calls for a species-by-sex interaction, and the loop goes round again.

Adding species and sex also changes the flipper coefficient, from 0.050 to about 0.020 kg per mm. Within a species and sex, flipper length is a weaker guide to mass than across all penguins, because much of the pooled relationship was between-group differences (the same mechanism as Simpson's paradox in Section 1.4). Neither number is "the effect of flippers on mass"; the graph told us there is no such effect. They are two different predictive summaries, answering different questions.

For prediction, we can also ask directly which model would predict new penguins better. Approximate leave-one-out cross-validation, computed by the loo package from the posterior draws without refitting the model, estimates each model's expected log predictive density for new data (Section 9.7 explains the importance-sampling idea behind it):

loo_compare(loo(m1), loo(m2))

   elpd_diff se_diff
m2   0.0       0.0  
m1 -93.6      11.2  

The richer model is better by about 94 units of expected log predictive density, with a standard error of 11: an overwhelming difference.

1.9 The workflow loop, and graphs that are not small

Putting the pieces together gives the loop below, which follows the "Bayesian workflow" of Gelman et al. (2020). It is a loop, not a pipeline: a failed check sends you back to the graph or the likelihood, never onward to the write-up.

Three principles travel with the loop. Build models you understand: start simple and add structure when a check demands it, as we did with species and sex. Computation is part of the model: a sampler that fails to converge has given you no answer. And criticism is not optional: the posterior predictive check is where the small world is held up against the large world, and skipping it skips the only honest test the method offers.

When the graph is not small

The graphs so far have three to five nodes. That is deliberate: a small graph isolates one structure so that its rule can be seen. But a DAG drawn for a question somebody actually needs answered routinely has ten or twenty nodes, several of them unobserved, mixing all four structures. The reasoning does not change; what changes is that you can no longer trace the paths in your head.

Here is a causal sketch for a real research question: do New York City neighbourhoods receive equitable emergency medical service? Specifically, does a neighbourhood's racial and ethnic composition $R$ affect the clinical outcome $O$ of residents after a 911 call, and does it do so by way of ambulance response time $RT$? The graph encodes that hypothesis and every rival explanation a sceptical reviewer would raise.

Read the graph as a set of claims. The latent node $L$ is the analyst's honesty made visible. Historical segregation and disinvestment shaped, together, the racial composition of neighbourhoods ($L\to R$), their present wealth ($L\to E$), and the siting of emergency infrastructure ($L\to W$). Population density $P$ drives traffic, distance to units and how busy ambulances are. Socioeconomic status $E$ shapes how severe calls are, through baseline health, and affects outcomes directly. Response time is produced by four operational variables and, possibly, by the exposure itself: the arrow $R\to RT$ is the disparity the project exists to detect.

The target is the total effect of $R$ on $O$, which flows along the causal paths from $R$ through $RT$ to $O$. Before reaching for software, try to find an adjustment set by hand. Four open back-door paths leave $R$ through the latent $L$: one through wealth ($R\leftarrow L\to E\to O$), two through wealth and call acuity ($R\leftarrow L\to E\to K\to O$ and $R\leftarrow L\to E\to K\to RT\to O$), and one through where stations were built ($R\leftarrow L\to W\to RT\to O$). Conditioning on $E$ blocks the first three and conditioning on $W$ blocks the fourth, so $\{E,W\}$ looks like the answer. Press condition on {E, W} in the widget.

It fails. Distance to the nearest unit has two parents, structural history and population density: $L\to W\leftarrow P$. So $W$ is a collider on the path $R\leftarrow L\to W\leftarrow P\to G\to RT\to O$, and conditioning on it opens that path (and a similar one through $V$). Blocking one path opened two others. This is exactly the kind of mistake that is easy to make by hand and impossible to make with software:

ems <- dagitty("dag {
  L [latent]
  L -> R   L -> E   L -> W
  P -> W   P -> G   P -> V
  E -> K   E -> O
  R -> RT
  W -> RT  K -> RT  G -> RT  V -> RT
  K -> O   RT -> O
}")
adjustmentSets(ems, exposure = "R", outcome = "O")
isAdjustmentSet(ems, c("E", "W"), exposure = "R", outcome = "O")
isAdjustmentSet(ems, c("E", "W", "P", "K", "G", "V", "RT"), exposure = "R", outcome = "O")

{ E, G, V, W }
{ E, P, W }
[1] FALSE
[1] FALSE

There are two minimal sufficient adjustment sets: $\{E,P,W\}$ (also condition on density, which closes the opened paths at their source) and $\{E,G,V,W\}$ (condition on traffic and availability instead). Which to use is a practical question: population density is measured precisely from census data, while traffic and ambulance availability at the moment of each call are harder to obtain. Neither set is the one that "control for everything available" would produce. That habit would condition on response time, the mediator, which erases the very effect under study, and because $RT$ is a collider with five parents, it would open still more non-causal paths. Try it in the widget. The graph turns a vague instinct into a precise instruction, and the instruction is different from both the instinct and the first careful guess.

One feature of the graph carries the deepest lesson. The set $\{E,P,W\}$ works only if structural history $L$ affects outcomes through no other route. Suppose it also raises residents' chronic exposure to environmental hazards, and so directly affects health: an arrow $L\to O$. Press add arrow L → O in the widget. A back-door path $R\leftarrow L\to O$ opens that no observed variable can block, because $L$ is unobserved.

ems_lo <- dagitty("dag {
  L [latent]
  L -> R   L -> E   L -> W   L -> O
  P -> W   P -> G   P -> V
  E -> K   E -> O
  R -> RT
  W -> RT  K -> RT  G -> RT  V -> RT
  K -> O   RT -> O
}")
adjustmentSets(ems_lo, exposure = "R", outcome = "O")   # prints nothing: no valid set
length(adjustmentSets(ems_lo, exposure = "R", outcome = "O"))

[1] 0

No adjustment set exists. The total effect is then not identified from these observational data, however large the sample. The right response is not to delete the arrow so that the graph looks solvable. It is to draw the arrow if you believe it belongs, and to say honestly that the effect can only be bounded, or reached through other designs such as instruments, natural experiments or sensitivity analysis (Chapter 13). A graph is allowed to tell you that your question cannot be answered with the data in hand. That message, delivered before any model is fitted, is one of the most valuable things the causal layer produces, and no Bayesian computation downstream can substitute for it.

1.10 Exercises

Each exercise has a hint and an answer button. Try the exercise first, open the hint if you are stuck, and look at the answer only to check your own.

Further reading and references

1. Berkson, J. (1946). Limitations of the application of fourfold table analysis to hospital data. Biometrics Bulletin, 2(3), 47–53.
2. FDA (2020). Vaccines and Related Biological Products Advisory Committee Briefing Document: Pfizer-BioNTech COVID-19 Vaccine, 10 December 2020.
3. Gelman, A., Vehtari, A., Simpson, D., et al. (2020). Bayesian workflow. arXiv:2011.01808.
4. Griffith, G. J., Morris, T. T., Tudball, M. J., et al. (2020). Collider bias undermines our understanding of COVID-19 disease risk and severity. Nature Communications, 11, 5749.
5. Horst, A. M., Hill, A. P. and Gorman, K. B. (2020). palmerpenguins: Palmer Archipelago (Antarctica) penguin data. R package. Data originally from Gorman, Williams and Fraser (2014), PLoS ONE 9(3), e90081.
6. Kruschke, J. K. (2021). Bayesian analysis reporting guidelines. Nature Human Behaviour, 5, 1282–1291.
7. McElreath, R. (2020). Statistical Rethinking: A Bayesian Course with Examples in R and Stan (2nd ed.). CRC Press.
8. Pearl, J. (1993). Comment: graphical models, causality and intervention. Statistical Science, 8(3), 266–269.
9. Pearl, J. (2009). Causality: Models, Reasoning, and Inference (2nd ed.). Cambridge University Press.
10. Polack, F. P., Thomas, S. J., Kitchin, N., et al. (2020). Safety and efficacy of the BNT162b2 mRNA Covid-19 vaccine. New England Journal of Medicine, 383, 2603–2615.
11. Tennant, P. W. G., Murray, E. J., Arnold, K. F., et al. (2021). Use of directed acyclic graphs (DAGs) to identify confounders in applied health research: review and recommendations. International Journal of Epidemiology, 50(2), 620–632.
12. Textor, J., van der Zander, B., Gilthorpe, M. S., Liśkiewicz, M. and Ellison, G. T. H. (2016). Robust causal inference using directed acyclic graphs: the R package 'dagitty'. International Journal of Epidemiology, 45(6), 1887–1894.
13. Vehtari, A., Gelman, A., Simpson, D., Carpenter, B. and Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: an improved $\widehat{R}$ for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667–718.
14. Wright, S. (1921). Correlation and causation. Journal of Agricultural Research, 20(7), 557–585.