Chapter 3
Bayesian linear models

Lectures 5–6 · the Gaussian model, the regression line, and what a coefficient means once the graph is drawn

a coefficient is an answer; the graph tells you which question it answers

More published science rests on the linear model than on any other statistical tool. It is the first model taught and the last one abandoned, and it is almost always the right place to begin. Its form is simple: the average of an outcome is a straight-line function of some predictors, and the outcomes scatter around that average with a bell-shaped spread. Every model in the rest of these notes, from logistic regression to neural networks, is a linear model with something changed.

The linear model's simplicity is also its trap. Fitting a line is easy; knowing what the fitted numbers mean is not. A regression of weight on height returns a slope of about 0.9 kg per centimetre. Is that how much heavier you would become if you grew a centimetre? How much heavier a taller person is, on average, among people like those in the sample? How much of the difference between men and women it explains? These are different questions, and the same coefficient answers at most one of them. Which one depends on the causal graph, not on the regression output. This chapter teaches the linear model and its Bayesian estimation, but its central lesson is the one in the cartoon: a coefficient is an answer, and the graph tells you which question it answers.

Lecture 5 builds the model. We will see why bell curves turn up everywhere, estimate the average height of American women from a national health survey (finding a trap in the data on the way), fit a line relating weight to height, separate the two kinds of uncertainty every prediction carries, and use a residual check to discover that the 19th-century formula for the body mass index is not quite right.
Lecture 6 asks what the coefficients mean. We will add a categorical predictor, separate a total effect from a direct one, watch a striking association between college majors and earnings shrink by half when we account for field of study, learn which "control variables" help and which do harm, and bend the line into a curve to measure when children grow fastest.

Think About It: Before You Start

Men in the United States are, on average, about 12 kg heavier than women. Men are also about 14 cm taller. Suppose you fit a regression of weight on sex and height. What do you expect the coefficient on sex to be, and what question would that number answer? Write down your guess; Section 3.7 gives the answer from data on 4,600 adults.

Lecture 5 — Lines, and the uncertainty around them

3.1 Why the bell curve keeps appearing

Measure the heights of a few thousand women and draw a histogram. You will get a bell: a single hump, roughly symmetric, falling away on both sides at a characteristic rate. Do the same for the lengths of beetles, the errors of a careful surveyor, the daily changes of a temperature reading, or the scores on a well-designed test, and you will keep getting bells. The normal (or Gaussian) distribution is so common that it is easy to forget it needs explaining. There are two explanations, and both matter for modelling.

The first is mechanical. A woman's adult height is the end result of many small influences: dozens of genetic variants each adding or subtracting a millimetre or two, nutrition in each year of childhood, illnesses, sleep, and luck. When an outcome is a sum of many small, roughly independent contributions, the contributions in each direction tend to cancel, the most common totals pile up in the middle, and extreme totals (which need nearly everything to push the same way) are rare. Whatever the shape of the individual contributions, the shape of the sum approaches the normal curve. This is the central limit theorem, and you can watch it work below.

Sums of small effects make bells; products make skewed bells

Each simulated person's outcome is built from $k$ small random influences. Choose how the influences are distributed and whether they add (each shifts the outcome) or multiply (each scales it by a few percent). The histogram shows 4,000 simulated people; the curve is the normal distribution with the same mean and standard deviation.

number of influences k = each influence is influences

Start with one flat influence: the histogram is flat, nothing like a bell. With two, it is a triangle. By five or six it is hard to tell from the curve, even when each influence is a coin flip or badly skewed. Now switch to multiply with many influences of ±5%. The result is a bell that leans to the right, with a long upper tail. That shape is the log-normal: when effects multiply, their logarithms add, so the logarithm of the outcome is normal. Body weight, incomes, and city sizes all grow multiplicatively and all have long right tails. Keep that in mind; it will matter in Section 3.5 and again in Section 3.8.

The second explanation is epistemic, and you met it in Section 2.8. If all you are prepared to say about a quantity is its mean and its variance, the normal distribution is the maximum-entropy distribution with that mean and variance: the one that assumes nothing else. Even when you have no idea what mechanism produced a variable, using a normal model for it commits you to nothing beyond a centre and a spread. That is why the normal is a sensible default likelihood for a continuous outcome, and also why it can fail: if the data have structure beyond a centre and a spread (two humps, a hard boundary, a long tail), the normal model will not see it. Section 3.5 shows how to look.

History: Galton's Board and the Origin of "Regression"

In the 1870s Francis Galton built a board studded with pins, the quincunx, and poured lead shot through it. Each ball bounced left or right at every pin, and the balls piled up at the bottom in a bell: a sum of coin flips, exactly as in the widget. Galton was fascinated by heredity, and in 1886 he published the heights of 928 adult children and their parents. Tall parents had tall children, but on average less extremely tall; short parents had children less extremely short. He called this "regression towards mediocrity", and the line he drew through the data became the first regression line. The phenomenon has nothing to do with heredity in particular: whenever two measurements are imperfectly correlated, extreme values of one are paired, on average, with less extreme values of the other. The name stuck to the whole method. Much earlier, in 1809, Carl Friedrich Gauss had derived the normal curve by asking which distribution of measurement errors would make the ordinary average the best estimate of a true value, an argument that was Bayesian in all but name, and had used it to justify least squares in astronomy.

Check Yourself

A biologist models the masses of individual fish in a lake with a normal distribution. A colleague objects that fish grow by a roughly constant percentage each season. Which model does the objection suggest?A normal distribution, because many small effects are involved.A uniform distribution between the smallest and largest fish.A normal distribution for the logarithm of mass (a log-normal for mass).

3.2 A model of one variable: how tall are American women?

Our data come from the National Health and Nutrition Examination Survey (NHANES), which the US Centers for Disease Control and Prevention runs continuously. Each year about 5,000 people are selected from across the country, interviewed at home, and then examined in a mobile clinic where trained staff measure height, weight, blood pressure and much more. We use the 2009–2012 cycles as packaged in the NHANES R package. First, the setup for the chapter.

library(tidyverse)
library(brms)
library(NHANES)     # US national health survey, 2009-2012
library(dagitty)
options(mc.cores = parallel::detectCores())   # run Markov chains in parallel

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())

Before modelling anything, look at how the data were made. Here is the first surprise.

nh <- NHANES |> select(ID, SurveyYr, Gender, Age, AgeMonths, Height, Weight)
c(rows = nrow(nh), people = n_distinct(nh$ID))
count(nh, ID, name = "copies") |> count(copies, name = "people")   # how often people repeat

adults <- nh |>
  distinct(ID, .keep_all = TRUE) |>        # keep one row per person
  filter(Age >= 20) |>
  drop_na(Height, Weight)
women <- filter(adults, Gender == "female")
c(adults = nrow(adults), women = nrow(women))

output

  rows people 
 10000   6779 
  copies people
1      1   4646
2      2   1413
3      3    461
4      4    177
5      5     60
6      6     18
7      7      3
8      8      1
adults  women 
  4609   2373

The package's data frame has 10,000 rows but only 6,779 different people: 2,133 people appear more than once, one of them eight times. This is not a mistake. NHANES deliberately oversamples some groups (older people, and Black, Hispanic and low-income Americans) so that it can estimate their health accurately, and every respondent carries a survey weight saying how many Americans they stand for. The package's authors turned the weighted survey into something that looks like a simple random sample by drawing 10,000 rows from it with replacement, in proportion to the weights. Heavily weighted people were drawn several times.

That is fine for drawing histograms, but it is poison for inference. A model that treats the rows as independent people believes it has 10,000 independent observations when it has 6,779, and its posterior will be too narrow. Repeating every person twice would make every standard deviation of the posterior too small by a factor of $\sqrt{2}$ . We keep one row per person, which leaves 4,609 adults with both height and weight measured, 2,373 of them women.

Research Practice: Know How Your Data Were Made

Reviewers of applied papers ask, before anything else, where the data came from and how the units were selected. Before you fit a model:

Count your units, not your rows. Check for duplicates (n_distinct(id) against nrow()), repeated measures, and multiple records per unit. Report the number of distinct units in the paper.
Read the documentation of every derived dataset. A teaching version of a survey, a Kaggle upload, or a colleague's "cleaned" file may differ from the source in ways that matter.
Know the design. National surveys such as NHANES, the American Community Survey and the British Household Panel use oversampling, clustering and weights. A population estimate that ignores them can be badly biased. The full answer for NHANES is to model the raw data with its weights and design, or to use multilevel regression and poststratification (Chapter 6). We will treat our de-duplicated sample as if it were a simple random sample of US adults, which is good enough to learn the methods and not good enough for a publication about the US population.

Check Yourself

A dataset accidentally contains every person exactly three times. By roughly what factor are the posterior standard deviations of a mean too small, if the duplicates are treated as independent people?

factor =

Now look at the heights. Adult heights of men and women together form a lumpy distribution, because they are a mixture of two groups with different averages. Within each sex, the bell shape is clear.

stats <- adults |> group_by(Gender) |> summarise(mean = mean(Height), sd = sd(Height))
stats
curves <- reframe(stats, Height = seq(135, 205, by = 0.5),
                  density = dnorm(Height, mean, sd), .by = Gender)
ggplot(adults, aes(Height)) +
  geom_histogram(aes(y = after_stat(density)), binwidth = 1,
                 fill = harbour[["sky"]], alpha = 0.55) +
  geom_line(data = curves, aes(y = density), colour = harbour[["amber"]], linewidth = 0.9) +
  facet_wrap(~ Gender) +
  labs(x = "height (cm)", y = "density")

output

# A tibble: 2 x 3
  Gender  mean    sd
  <fct>  <dbl> <dbl>
1 female  162.  7.36
2 male    175.  7.58

Heights of 4,609 US adults (NHANES 2009–2012), with the normal curve that has each group's mean and standard deviation. Within each sex the normal describes the data closely.

We will start with the women. The question is simple: what is the average height of adult women in the United States, and how sure can we be? The model says that each woman's height $h_{i}$ is a draw from a normal distribution with unknown mean $μ$ and unknown standard deviation $σ$ :

h_{i} \sim Normal (μ, σ), μ \sim Normal (170, 20), σ \sim Exponential (0.1) .

Read the priors as statements of knowledge you had before seeing these data. The prior for $μ$ is centred at 170 cm with a standard deviation of 20 cm, so it considers average heights from about 130 cm to 210 cm plausible: wide enough to include every human population ever measured, but not absurd. The prior for $σ$ is the maximum-entropy distribution for a positive quantity with a known mean (Section 2.8); its mean is $1 / 0.1 = 10$ cm, which says "people differ in height by something like ten centimetres", and it gives very little weight to spreads above 40 cm. With 2,373 women, neither prior will matter much. In the small subgroups you meet in research, they will.

Definition: The Gaussian Model

A Gaussian model for an outcome $y$ says that each observation is drawn independently from a normal distribution, $y_{i} \sim Normal (μ, σ)$ . Its two parameters are the mean $μ$ , the centre of the distribution, and the standard deviation $σ$ , the typical distance of an observation from the centre. About 68% of observations fall within one $σ$ of $μ$ and about 95% within two. A Bayesian analysis adds priors for $μ$ and $σ$ and computes their joint posterior.

The posterior, computed by brute force

With two parameters, we can compute the posterior by the method of Chapter 1 with no cleverness at all. Lay a fine grid over all plausible pairs $(μ, σ)$ ; at each pair, compute the log-likelihood (the sum of the log normal densities of all 2,373 heights) and add the log priors; then exponentiate and normalise. This is grid approximation.

loglik <- function(mu, sigma, y) sum(dnorm(y, mu, sigma, log = TRUE))
grid <- expand_grid(mu    = seq(160.9, 162.3, length.out = 150),
                    sigma = seq(6.9, 7.85, length.out = 150)) |>
  mutate(loglik  = map2_dbl(mu, sigma, loglik, y = women$Height),
         logpost = loglik + dnorm(mu, 170, 20, log = TRUE) + dexp(sigma, 0.1, log = TRUE),
         post    = exp(logpost - max(logpost)),   # subtract the max first, to avoid underflow
         post    = post / sum(post))

set.seed(5)
grid_draws <- slice_sample(grid, n = 1e4, weight_by = post, replace = TRUE)
grid_draws |>
  summarise(across(c(mu, sigma), list(mean = mean,
                                      lo = ~ quantile(.x, 0.025), hi = ~ quantile(.x, 0.975))))

output

# A tibble: 1 x 6
  mu_mean mu_lo mu_hi sigma_mean sigma_lo sigma_hi
    <dbl> <dbl> <dbl>      <dbl>    <dbl>    <dbl>
1    162.  161.  162.       7.36     7.16     7.57

Two things in that code are worth copying into every grid you write. The likelihood of 2,373 observations is a product of 2,373 small numbers, about $10^{- 3500}$ , far below the smallest number a computer can store, so we work with log-likelihoods and subtract the maximum before exponentiating. And instead of summarising the grid directly, we sample 10,000 parameter pairs from it in proportion to their posterior probability. From then on, every question about the posterior is a question about a table of samples, exactly as in Chapter 1.

ggplot(grid, aes(mu, sigma)) +
  geom_point(data = slice_head(grid_draws, n = 600), colour = harbour[["sky"]],
             alpha = 0.25, size = 0.9,
             position = position_jitter(width = 0.004, height = 0.003, seed = 1)) +
  geom_contour(aes(z = post), colour = harbour[["amber"]], bins = 7, linewidth = 0.6) +
  labs(x = "mean height μ (cm)", y = "standard deviation σ (cm)")

The joint posterior of the mean and standard deviation of US women's height, computed on a 150 × 150 grid. Contours enclose regions of high posterior probability; dots are 600 of the 10,000 samples drawn from the grid.

The posterior is a small, almost circular hill. The average height of US women is between about 161.3 and 161.9 cm with 95% probability, and the standard deviation between about 7.16 and 7.57 cm. The hill is round rather than tilted, which says that what we learn about the centre tells us almost nothing about the spread: $μ$ and $σ$ are nearly uncorrelated in the posterior. The widget lets you watch the hill form as the data arrive.

With one woman measured, the posterior for $σ$ is still essentially the prior (you cannot learn about spread from one observation), and $μ$ is pinned only loosely near her height. By thirty, both are well located. Notice also that with few observations the hill is lopsided, with a longer tail toward large $σ$ : small samples cannot rule out a large spread. The normal approximation to a posterior, which many software packages report, is only trustworthy once the hill is round.

The same model in brms

Grids are perfect for two parameters and hopeless for twenty: a grid with 150 points per parameter needs $150^{20}$ points. Every serious model in this course is therefore fitted with brms, which translates a model formula into the Stan language and draws posterior samples by Hamiltonian Monte Carlo (Chapter 5 explains how). The formula Height ~ 1 means "height has a mean and nothing predicts it".

m_h <- brm(Height ~ 1, data = women,
           prior = c(prior(normal(170, 20), class = Intercept),
                     prior(exponential(0.1), class = sigma)),
           seed = 3, refresh = 0, file = "fits/ch3_height")
summary(m_h)

output

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: Height ~ 1 
   Data: women (Number of observations: 2373) 
  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   161.58      0.15   161.28   161.87 1.00     4056     2702

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     7.36      0.11     7.16     7.58 1.00     3598     2797

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

The two methods agree to within a few hundredths of a centimetre. In the brms summary, Estimate is the posterior mean, Est.Error the posterior standard deviation, and l-95% CI and u-95% CI the ends of the 95% credible interval. Rhat and the two ESS columns diagnose the sampler: an Rhat of 1.00 and effective sample sizes in the thousands mean the chains agree with each other and the summaries are precise. You met these in Section 1.8, and Chapter 5 explains them properly.

3.3 Drawing a line: weight and height

A mean is a model with no predictors. The linear model lets the mean depend on something. Taller people tend to be heavier; how much heavier? Before writing a model, draw the graph.

A causal sketch for adult weight. Height H affects weight W (a taller body has more bone and tissue to carry). Sex S affects both. Many other influences U on both (genes, childhood nutrition) are unmeasured. Our first model ignores S; Section 3.7 brings it back.

The first model is the simplest one that can answer "how does average weight change with height?":

W_{i}

\sim Normal (μ_{i}, σ)

μ_{i}

= α + β (H_{i} - 170)

α

\sim Normal (80, 20)

β

\sim Normal (0, 2)

σ

\sim Exponential (0.1)

The first line is the likelihood, as before, except that the mean now carries an index $i$ : each person has their own expected weight $μ_{i}$ . The second line is the linear model: $μ_{i}$ is a straight-line function of height. The equals sign is deliberate. The relationship between $μ_{i}$ and $H_{i}$ is deterministic; all the randomness lives in the first line.

Each prior is a statement you can defend before seeing the data:

$α$ is the expected weight of a person exactly 170 cm tall, because we centred height at 170 cm. A 170 cm adult weighs perhaps 60 to 100 kg, and Normal(80, 20) covers that comfortably.
$β$ is the change in expected weight, in kilograms, per centimetre of height. Normal(0, 2) allows anything from about −4 to 4 kg per cm. It does not even assume the slope is positive, and it rules out only the absurd (a slope of 10 would make a 190 cm person 200 kg heavier than a 170 cm one).
$σ$ is the typical deviation of a person's weight from the line, again with an exponential prior of mean 10 kg.

Definition: The Linear Model

A linear model for an outcome $y$ with predictors $x_{1}, \dots, x_{K}$ is

y_{i} \sim Normal (μ_{i}, σ), μ_{i} = α + β_{1} x_{1 i} + \dots + β_{K} x_{K i} .

"Linear" means linear in the parameters $α$ and $β_{k}$ . The predictors can be any fixed functions of the data: logarithms, squares, indicators of categories, or the spline basis functions of Section 3.10. A coefficient $β_{k}$ is the difference in the expected outcome between two units whose $x_{k}$ differ by one and whose other predictors are equal. Whether that difference is a causal effect depends on the graph.

Why centre the predictor?

Without centring, $α$ would be the expected weight of a person of height zero, a quantity no one has prior knowledge about, and whose posterior is strongly correlated with $β$ (raise the slope a little and the intercept must fall a lot to keep the line through the cloud of data). Centring at a typical value makes $α$ interpretable, makes its prior easy to set, and makes $α$ and $β$ nearly independent in the posterior, which also helps the sampler. Centre at a round, meaningful number (170 cm, age 40, the year 2000) rather than at the sample mean, so that the intercept means the same thing in every analysis and every paper. Exercise 4 asks you to fit the model without centring and see what happens.

adults <- mutate(adults, h = Height - 170)      # centimetres above 170
m_wh <- brm(Weight ~ 1 + h, data = adults,
            prior = c(prior(normal(80, 20), class = Intercept),
                      prior(normal(0, 2), class = b),
                      prior(exponential(0.1), class = sigma)),
            seed = 3, refresh = 0, file = "fits/ch3_wh")
fixef(m_wh) |> round(3)
posterior_summary(m_wh, variable = "sigma") |> round(2)

output

          Estimate Est.Error  Q2.5  Q97.5
Intercept   83.602     0.285 83.03 84.166
h            0.903     0.027  0.85  0.958
      Estimate Est.Error  Q2.5 Q97.5
sigma    19.14       0.2 18.76 19.52

The posterior says: a 170 cm adult weighs 83.6 kg on average (95% interval 83.0 to 84.2), and each extra centimetre of height goes with 0.90 kg more weight on average (0.85 to 0.96). Individual weights scatter around the line with a standard deviation of about 19 kg, which is large: height explains only part of why people differ in weight.

Read the slope carefully. It compares the average weights of different people whose heights differ by a centimetre. It does not say what happens to one person who grows a centimetre, and, as the graph warns, it mixes the effect of height with the effect of sex, because men are both taller and heavier. Section 3.7 untangles this.

3.4 Two kinds of uncertainty

A fitted line is not one line but a posterior distribution of lines, one for each posterior sample of $(α, β)$ . Any question about the line is answered by asking it of every sample. There are two questions people ask most often, and they have very different answers.

Where is the average? For a given height, where is the mean weight $μ$ of people of that height? Compute $α + β (H - 170)$ for each posterior sample and summarise. This is the credible interval for the mean, and it reflects uncertainty about the parameters only.
Where is a person? For a given height, where will a new individual's weight fall? For each posterior sample, also draw a weight from Normal $(μ, σ)$ . This is the posterior predictive interval, and it adds the scatter of individuals around the mean.

grid_h <- tibble(h = seq(-30, 30, by = 1))
mu_band   <- fitted(m_wh, newdata = grid_h) |> as_tibble() |> bind_cols(grid_h)    # the mean
pred_band <- predict(m_wh, newdata = grid_h, seed = 4) |> as_tibble() |> bind_cols(grid_h)  # a person

set.seed(2)
ggplot(slice_sample(adults, n = 1000), aes(h + 170, Weight)) +     # plot 1,000 of the 4,609
  geom_point(colour = harbour[["muted"]], alpha = 0.25, size = 0.8) +
  geom_ribbon(data = pred_band, aes(x = h + 170, ymin = Q2.5, ymax = Q97.5),
              inherit.aes = FALSE, fill = harbour[["sky"]], alpha = 0.18) +
  geom_ribbon(data = mu_band, aes(x = h + 170, ymin = Q2.5, ymax = Q97.5),
              inherit.aes = FALSE, fill = harbour[["amber"]], alpha = 0.6) +
  geom_line(data = mu_band, aes(h + 170, Estimate), colour = harbour[["amber"]]) +
  labs(x = "height (cm)", y = "weight (kg)")

Weight against height (a random 1,000 of the 4,609 adults are shown; the model uses all of them), with the posterior mean line, the 95% credible interval for the mean weight at each height (narrow amber band), and the 95% posterior predictive interval for an individual's weight (wide blue band).

The amber band is barely wider than the line: with 4,609 people, we know the average weight at each height to within a kilogram or so. The blue band is about 75 kg wide: we cannot predict an individual's weight from height to better than that. Both statements are correct, and confusing them is one of the commonest errors in applied work. A paper that reports only the narrow band can make a weak predictor look precise; a paper that reports only the wide one can hide a well-estimated relationship.

Slide $n$ upward. The posterior lines converge and the amber band shrinks toward the line: more data always buys more certainty about the mean, roughly in proportion to $1 / \sqrt{n}$ . The blue band moves around at first, while $σ$ is poorly known (with a handful of people it can even be too narrow, if they happen to sit close to the line), and then settles: it can never be narrower than the spread of real people around the line. No amount of data makes height a good predictor of a particular person's weight. With very few points, notice how wide both bands flare at the edges: a line is pinned down best near the middle of the data and least well at the extremes, so extrapolation is where uncertainty grows fastest.

Research Practice: Say Which Interval You Mean

Every interval in a paper should say what it is an interval for: a parameter, the mean outcome at given predictor values, or a new observation. In captions write, for example, "shaded band: 95% credible interval for the expected weight" or "95% posterior predictive interval". When the purpose is prediction (forecasting response times, flagging a risky login, predicting a student's grade), report predictive intervals and their calibration (do 95% intervals contain about 95% of held-out cases?). When the purpose is explanation, report intervals for the parameters or contrasts that answer the scientific question.

Check Yourself

Doubling the number of people in the NHANES weight–height analysis would make the 95% predictive interval for a new person's weight

In the widget with

n = 10

, the bands are much wider at heights of 150 and 195 cm than near 170 cm. Why?The slope is uncertain, so lines that pivot around the centre of the data spread apart far from it.People at extreme heights have more variable weights.Height was centred at 170 cm, which forces the bands to be narrowest there.

3.5 Checking the model, and the shape of the human body

A posterior is a statement about the parameters of the model you wrote down. If the model is wrong in an important way, the posterior can be precise and misleading at once. Section 1.8 introduced the remedy: simulate data from the fitted model and compare them with the real data. For the weight model, compare the shape of the distribution of weights.

set.seed(6)
yrep <- posterior_predict(m_wh, ndraws = 200)                 # 200 simulated datasets
skew <- function(x) mean((x - mean(x))^3) / sd(x)^3
c(observed = skew(adults$Weight), simulated = mean(apply(yrep, 1, skew)),
  share_simulated_as_skewed = mean(apply(yrep, 1, skew) >= skew(adults$Weight)))

sims <- tibble(draw = rep(1:25, each = nrow(adults)), Weight = as.vector(t(yrep[1:25, ])))
ggplot() +
  geom_density(data = sims, aes(Weight, group = draw), colour = harbour[["sky"]],
               alpha = 0.3, linewidth = 0.3, n = 128) +
  geom_density(data = adults, aes(Weight), colour = harbour[["ink"]], linewidth = 1.1, n = 256) +
  labs(x = "weight (kg)", y = "density")

output

                 observed                 simulated share_simulated_as_skewed 
                 1.106532                  0.007017                  0.000000

Posterior predictive check for the linear model of weight. Thin blue curves: the distribution of weights in 25 datasets simulated from the fitted model. Thick curve: the real data, which have a sharper peak and a longer right tail than the model can produce.

The check fails clearly. The real weights are skewed: there are more people far above the average weight than far below it. Simulated datasets from the model are symmetric, with a skewness near 0, while the real data have a skewness of about 1.1, and not one of 200 simulated datasets comes close. The line itself may still be a fine summary of the mean, but the model's statements about individuals (such as the predictive band in Section 3.4) are wrong in the tails: it predicts too many very light people and too few very heavy ones.

Section 3.1 suggests the fix. Weight is built up multiplicatively: tissue grows in proportion to the tissue already there, and a 10% change in diet changes weight by a percentage, not by a fixed number of kilograms. So model log weight. If we also take the logarithm of height, the model becomes

\log W_{i} \sim Normal (μ_{i}, σ), μ_{i} = α + β \log (H_{i} / 170) .

On the original scale this says $W \approx e^{α} (H / 170)^{β}$ : weight is proportional to height raised to the power $β$ . Logarithms on both sides turn a power law into a straight line, and the slope $β$ becomes an elasticity: a 1% difference in height goes with a $β$ % difference in weight.

History: Quetelet's Index and the Exponent of Two

In 1835 the Belgian astronomer and statistician Adolphe Quetelet measured Belgian conscripts and claimed that adult weight increases roughly as the square of height. Dividing weight by height squared should therefore remove the effect of height, leaving a measure of how heavy someone is "for their size". In 1972 the physiologist Ancel Keys compared candidate indices and recommended Quetelet's, naming it the body mass index, BMI $= W / H^{2}$ . It now decides medical categories for billions of people. But why a power of two? If people were scaled-up copies of each other, weight would grow as the cube of height. In 2013 the Oxford mathematician Nick Trefethen argued in a letter to The Economist that an exponent of 2.5 fits real populations better, so that BMI misclassifies tall people as fatter, and short people as thinner, than they are. The exponent is an empirical question, and a linear model on the log scale answers it directly.

A model of $\log W$ on $\log H$ puts Quetelet's claim to the test: if BMI is the right adjustment, then $β = 2$ . We fit one model to all adults, and a second that allows the exponent to differ between women and men. The prior for $β$ is centred halfway between Quetelet's 2 and the 3 of geometric scaling, Normal(2.5, 1), so that it favours neither hypothesis.

adults <- mutate(adults, log_w = log(Weight), log_h = log(Height / 170))
m_log <- brm(log_w ~ 1 + log_h, data = adults,
             prior = c(prior(normal(4.4, 0.5), class = Intercept),   # log(80 kg) is 4.38
                       prior(normal(2.5, 1), class = b),
                       prior(exponential(5), class = sigma)),
             seed = 3, refresh = 0, file = "fits/ch3_log")
m_logs <- brm(log_w ~ 0 + Gender + Gender:log_h, data = adults,     # separate lines by sex
              prior = c(prior(normal(4.4, 0.5), class = b),
                        prior(normal(2.5, 1), class = b, coef = Genderfemale:log_h),
                        prior(normal(2.5, 1), class = b, coef = Gendermale:log_h),
                        prior(exponential(5), class = sigma)),
              seed = 3, refresh = 0, file = "fits/ch3_logs")

expo <- as_draws_df(m_logs) |>
  transmute(women = `b_Genderfemale:log_h`, men = `b_Gendermale:log_h`,
            `men - women` = men - women) |>
  mutate(`all adults` = as_draws_df(m_log)$b_log_h)
expo |>
  pivot_longer(everything(), names_to = "exponent") |>
  group_by(exponent) |>
  summarise(mean = mean(value), lo = quantile(value, 0.025), hi = quantile(value, 0.975),
            `Pr(< 2)` = mean(value < 2))

output

# A tibble: 4 x 5
  exponent     mean     lo    hi `Pr(< 2)`
  <chr>       <dbl>  <dbl> <dbl>     <dbl>
1 all adults  1.90  1.79   2.00      0.971
2 men         2.08  1.87   2.28      0.238
3 men - women 0.364 0.0777 0.648     1    
4 women       1.71  1.52   1.91      0.998

expo |>
  select(-`men - women`) |>
  pivot_longer(everything(), names_to = "group") |>
  ggplot(aes(value, colour = group)) +
  geom_density(linewidth = 1) +
  geom_vline(xintercept = c(2, 2.5), linetype = "dashed", colour = harbour[["muted"]]) +
  scale_colour_manual(values = c(women = harbour[["rose"]], men = harbour[["sky"]],
                                 `all adults` = harbour[["amber"]])) +
  labs(x = "exponent β  (weight ∝ height^β)", y = "density", colour = NULL)

Posterior distributions of the exponent β in weight ∝ height^β, for women, for men, and for all adults together. Dashed lines mark Quetelet's exponent of 2 (BMI) and Trefethen's proposed 2.5.

The answer is more interesting than either side of the argument. For all adults together, the exponent is about 1.90 (95% interval 1.79 to 2.00), slightly below Quetelet's 2, with posterior probability 0.97. But that pooled number hides a real difference: among women the exponent is about 1.71 (95% interval 1.52 to 1.91), among men about 2.08 (1.87 to 2.28), and the difference, 0.36 (0.08 to 0.65), is clearly positive. Neither sex is close to Trefethen's 2.5 in these data. A BMI built for men would be about right; for women, dividing by $H^{2}$ over-corrects, so that among women of the same build, taller women get slightly lower BMIs. (Whether that matters clinically is a question about how BMI predicts health, not one this model can answer.)

Notice also how the pooled exponent, 1.90, is not the average of the two sex-specific exponents in any simple way. Pooling men and women puts two clouds of points with different centres into one regression, and the line through the combined cloud reflects both the within-sex slopes and the difference between the groups. This is a gentle version of the confounding problem that dominates the next lecture.

Does the log scale fix the skew? The skewness of the residuals (each person's distance from the fitted line) falls from 1.35 on the kilogram scale to about 0.45 on the log scale. Not perfect (weight has a longer right tail than even a log-normal), but the model is now much closer to the data, and Exercise 5 asks you to run the posterior predictive check yourself. A model does not need to be perfect; it needs to be good enough for the question, and you need to know where it fails.

Research Practice: Report the Checks You Ran

Include at least one posterior predictive check for each main model in the paper or its supplement, chosen to probe the assumptions that matter for your conclusions: tails for risk predictions, group means for comparisons, skewness or zeros for counts. Say what failed and what you did about it. "The Gaussian model under-predicted heavy weights (posterior predictive skewness 0.0 against 1.1 observed), so we modelled log weight" is a sentence that builds a reviewer's trust. Choosing the outcome scale after looking at the data is a modelling decision like any other; report it, and if it could change the conclusion, report both analyses.

Check Yourself

In the model for men, the exponent is about 2.08. By about what percentage is the average weight of men who are 1% taller higher? (Use the elasticity interpretation.)

percent =%

The posterior probability that the exponent is below 2 is essentially 1 for women and about 0.24 for men. What is the most accurate summary?BMI is correct for men and wrong for women.The data favour an exponent below 2 for women; for men they are consistent with 2, though values somewhat above 2 are more likely.The exponent is 2 for both sexes, within error.

Key Takeaways: Lecture 5

Sums of many small effects are approximately normal; products are log-normal. The normal is also the maximum-entropy choice given a mean and a variance.
Count units, not rows, and know how the data were sampled before modelling.
A linear model lets the mean depend on predictors. Centre predictors at meaningful values and set priors you can defend on that scale.
The credible interval for the mean shrinks with data; the predictive interval for an individual does not. Say which one you report.
Check the model with posterior predictive checks. Logarithms turn multiplicative relationships and power laws into lines, and their slopes are elasticities.

Lecture 6 — What a coefficient means

3.6 Categories as predictors

Not every predictor is a number. Sex, field of study, treatment arm, borough and survey year are categorical. There are two ways to put a category into a linear model. The traditional way uses an indicator (or "dummy") variable: $M_{i} = 1$ for men and $0$ for women, and $μ_{i} = α + β M_{i}$ . Then $α$ is the women's mean and $β$ is the difference between men and women.

We will use the second way, an index variable: give each category its own parameter, and let the data's category label pick which one applies,

μ_{i} = α_{S [i]}, α_{female} \sim Normal (80, 20), α_{male} \sim Normal (80, 20) .

The two codings fit the same data equally well, but the index version is easier to think about. With indicators, the prior on $β$ is a prior on a difference, so the men's mean $α + β$ automatically gets a wider prior than the women's (it is the sum of two uncertain quantities), which is an assumption no one intended. With an index variable each category gets the same prior, the parameters mean what they say, and the extension to 16 categories (Section 3.8) needs no new ideas. In brms, the formula Weight ~ 0 + Gender asks for one parameter per category and no common intercept.

m_total <- brm(Weight ~ 0 + Gender, data = adults,
               prior = c(prior(normal(80, 20), class = b),
                         prior(exponential(0.1), class = sigma)),
               seed = 3, refresh = 0, file = "fits/ch3_total")
fixef(m_total) |> round(2)

output

             Estimate Est.Error  Q2.5 Q97.5
Genderfemale    76.21      0.41 75.40 77.03
Gendermale      88.21      0.44 87.35 89.09

How much heavier are men than women, on average? The model has no parameter for that, and it does not need one. A difference between categories is a function of the parameters, and as always we compute it from the posterior samples: subtract, sample by sample. This is called a contrast.

post <- as_draws_df(m_total)
diff_total <- post$b_Gendermale - post$b_Genderfemale     # one difference per posterior sample
c(mean = mean(diff_total), quantile(diff_total, c(0.025, 0.975)))

output

 mean  2.5% 97.5% 
12.01 10.79 13.17

Research Practice: Compare Groups with Contrasts, Not by Eye

A very common error is to compare two groups by checking whether their separate intervals overlap, or whether one estimate is "significant" and the other is not. Overlapping intervals do not mean there is no difference, and "the difference between significant and not significant is not itself statistically significant" (Gelman and Stern, 2006). The correct summary of a comparison is the posterior distribution of the contrast itself, computed draw by draw as above, with its interval. Every comparison your paper's argument depends on should appear as a contrast.

3.7 Total and direct effects

Men weigh about 12 kg more than women on average. Now return to the question in the opening box: what happens when we add height to the model? Look at the graph again, keeping only the parts we need.

Sex S affects weight W along two paths: directly, and through height H (the mediator). The total effect of S on W combines both paths; the direct effect is the part that does not go through H. U, unmeasured influences on both H and W, matters for the direct effect (see the warning below).

Sex affects weight in two ways: through height (men are taller, and taller bodies are heavier), and directly, by any other route (differences in muscle mass, fat distribution and bone density at the same height). These give two different questions, and two different estimands:

The total effect of sex on weight: how different are the average weights of men and women? It includes both paths. The graph has no back-door path from $S$ to $W$ (nothing causes sex in this graph), so the model with sex alone estimates it. That was m_total.
The direct effect: how different are the average weights of men and women of the same height? To estimate it, we block the path through $H$ by conditioning on height.

m_direct <- brm(Weight ~ 0 + Gender + h, data = adults,
                prior = c(prior(normal(80, 20), class = b),
                          prior(normal(0, 2), class = b, coef = h),
                          prior(exponential(0.1), class = sigma)),
                seed = 3, refresh = 0, file = "fits/ch3_direct")
post_d <- as_draws_df(m_direct)
diff_direct <- post_d$b_Gendermale - post_d$b_Genderfemale   # men minus women, same height

tibble(total = diff_total, direct = diff_direct) |>
  pivot_longer(everything(), names_to = "estimand") |>
  group_by(estimand) |>
  summarise(mean = mean(value), lo = quantile(value, 0.025), hi = quantile(value, 0.975))

output

# A tibble: 2 x 4
  estimand   mean    lo     hi
  <chr>     <dbl> <dbl>  <dbl>
1 direct   -0.790 -2.32  0.718
2 total    12.0   10.8  13.2

tibble(total = diff_total, direct = diff_direct) |>
  pivot_longer(everything(), names_to = "estimand") |>
  ggplot(aes(value, fill = estimand)) +
  geom_density(colour = NA, alpha = 0.7) +
  geom_vline(xintercept = 0, linetype = "dashed", colour = harbour[["muted"]]) +
  scale_fill_manual(values = c(total = harbour[["sky"]], direct = harbour[["amber"]])) +
  labs(x = "men minus women, average weight (kg)", y = "density", fill = NULL)

Posterior distributions of the difference in average weight between men and women. Total effect: about 12 kg. Direct effect at the same height: indistinguishable from zero. Essentially all of the weight difference between the sexes runs through height.

The total effect is 12.0 kg. The direct effect, comparing men and women of the same height, is about −0.8 kg, with an interval comfortably including zero. In these data, a man and a woman of the same height weigh, on average, about the same. All of the average weight difference between the sexes is accounted for by the difference in height.

Both numbers are correct. They answer different questions, and the question comes first. If you are a public-health planner estimating how much more food a male population needs, the total effect is the relevant one. If you are asking whether men are heavier "for their size", the direct effect is. What you must not do is fit the model with height, read off the sex coefficient, and report that "sex has no effect on weight". The coefficient in m_direct is a direct effect, and that is not what the sentence claims.

Warning: Direct Effects Need More Than a Graph Without Back Doors

The total effect needs only that nothing confounds $S$ and $W$ . The direct effect needs more: nothing may confound the mediator and the outcome. In the graph, unmeasured influences $U$ (genes, childhood nutrition) affect both height and weight. Conditioning on height then makes it a collider on the path $S \to H \leftarrow U \to W$ , and opens that path. Among people of the same height, a man is, on average, someone whose $U$ pushed his height less (he is "short for a man") and a woman someone whose $U$ pushed hers more; if the same $U$ also affects weight, the comparison is biased. The near-zero direct effect is therefore an estimate of a controlled comparison at equal height, which is a perfectly good descriptive quantity, but its causal reading depends on an assumption about $U$ . Chapter 13 develops the tools (mediation analysis and sensitivity analysis) for doing better.

Check Yourself

A university finds that students who use an AI tutor have higher exam scores. A reviewer asks the authors to "control for" the number of hours each student spent studying, which the tutor is thought to increase. What will the coefficient on tutor use then estimate?The total effect of the tutor, with less confounding.The effect of hours studied.The direct effect of the tutor at a fixed number of study hours, which leaves out the part of its effect that works by making students study more.

3.8 Confounding in the wild: majors and money

Now a harder example, with real stakes. In 2014 the data journalism site FiveThirtyEight published earnings data for 173 college majors, taken from the American Community Survey for 2010–2012: for each major, the median annual earnings of recent graduates working full time, and the share of graduates who are women. The question that leaps out is: do majors with more women pay less, and why?

url <- "https://raw.githubusercontent.com/fivethirtyeight/data/master/college-majors/recent-grads.csv"
majors <- read_csv(url, show_col_types = FALSE) |>
  drop_na(ShareWomen) |>                     # one major has no breakdown by sex
  mutate(log_earn = log(Median), field = Major_category)
nrow(majors)
majors |> select(Major, field, ShareWomen, Median, Sample_size) |> slice(c(1, 40, 90, 160))

output

[1] 172
# A tibble: 4 x 5
  Major                                  field       ShareWomen Median Sample_size
  <chr>                                  <chr>            <dbl>  <dbl>       <dbl>
1 PETROLEUM ENGINEERING                  Engineering      0.121 110000          36
2 ACCOUNTING                             Business         0.524  45000        2042
3 GEOSCIENCES                            Physical S~      0.591  36000          18
4 COSMETOLOGY SERVICES AND CULINARY ARTS Industrial~      0.585  29000         117

The raw association is strong: across majors, the correlation between the share of women and median earnings is −0.62. Petroleum engineering (12% women) pays $110,000; early childhood education (97% women) pays $28,000. We model the logarithm of earnings, because earnings are multiplicative and skewed (Section 3.5), and a coefficient on the log scale converts to a percentage difference.

Before fitting, the graph. Majors belong to broad fields (engineering, health, education, the arts, and so on). Field plausibly affects both the share of women (through decades of gendered sorting into fields, and the different histories of nursing and engineering) and earnings (through the jobs each field leads to). That makes field a common cause, a fork, and the raw association mixes any link between share of women and pay with the differences between fields.

A causal sketch for majors. Field F affects both the share of women SW and median earnings E. The question is about the arrow from SW to E (any pay penalty attached to majors with more women). U stands for unmeasured features of each major, such as the specific occupations it leads to.

We fit two models: one ignoring field, and one with a separate intercept for each of the 16 fields (an index variable, as in Section 3.6) and a common slope for the share of women. Formally, the second is

\log E_{i} \sim Normal (μ_{i}, σ), μ_{i} = α_{F [i]} + β S W_{i},

with Normal(10.6, 1) priors for the field intercepts ( $e^{10.6} \approx$ $40,000, and one unit on the log scale is a factor of 2.7, so this allows field averages from about $5,000 to $300,000) and a Normal(0, 1) prior for $β$ , which allows majors with no women and majors with only women to differ in earnings by up to a factor of about 7 in either direction.

m_maj1 <- brm(log_earn ~ 1 + ShareWomen, data = majors,
              prior = c(prior(normal(10.6, 1), class = Intercept),     # log($40,000)
                        prior(normal(0, 1), class = b),
                        prior(exponential(2), class = sigma)),
              seed = 3, refresh = 0, file = "fits/ch3_maj1")
m_maj2 <- brm(log_earn ~ 0 + field + ShareWomen, data = majors,
              prior = c(prior(normal(10.6, 1), class = b),
                        prior(normal(0, 1), class = b, coef = ShareWomen),
                        prior(exponential(2), class = sigma)),
              seed = 3, refresh = 0, file = "fits/ch3_maj2")

# percentage difference in median earnings, all-women major vs all-men major
pct <- function(fit) 100 * (exp(as_draws_df(fit)$b_ShareWomen) - 1)
tibble(`ignoring field` = pct(m_maj1), `within field` = pct(m_maj2)) |>
  pivot_longer(everything(), names_to = "model") |>
  group_by(model) |>
  summarise(mean = mean(value), lo = quantile(value, 0.025), hi = quantile(value, 0.975))

output

# A tibble: 2 x 4
  model           mean    lo    hi
  <chr>          <dbl> <dbl> <dbl>
1 ignoring field -50.3 -56.0 -44.0
2 within field   -30.5 -40.3 -18.7

Ignoring field, a major with only women graduates is predicted to pay about half (50% less, 95% interval 44% to 56%) what a major with only men pays. Within the same field, the gap is about 30% (95% interval 19% to 40%): still large, and with an interval well away from zero, but only about half as large on the log scale (about −0.36 against −0.70). Roughly half of the raw association is a difference between fields; the rest is a difference between majors within a field. The figure shows where each part comes from.

fe <- fitted(m_maj2) |> as_tibble() |> bind_cols(majors)          # within-field fits
overall <- tibble(ShareWomen = c(0, 1)) |>
  bind_cols(as_tibble(fitted(m_maj1, newdata = tibble(ShareWomen = c(0, 1)))))
centres <- majors |> group_by(field) |>
  summarise(ShareWomen = mean(ShareWomen), log_earn = mean(log_earn), n = n())

ggplot(majors, aes(ShareWomen, exp(log_earn))) +
  geom_point(colour = harbour[["muted"]], alpha = 0.5, size = 1.4) +
  geom_line(data = fe, aes(ShareWomen, exp(Estimate), group = field),
            colour = harbour[["sky"]], linewidth = 0.8) +
  geom_line(data = overall, aes(ShareWomen, exp(Estimate)),
            colour = harbour[["rose"]], linewidth = 1.2) +
  geom_point(data = centres, aes(size = n), colour = harbour[["amber"]], alpha = 0.9) +
  ggrepel::geom_text_repel(data = filter(centres, field %in% c("Engineering", "Education", "Health", "Arts")),
            aes(label = field), nudge_y = 0.12, size = 3, colour = harbour[["ink"]], min.segment.length = 0,
            segment.colour = harbour[["muted"]], seed = 1) +
  scale_y_log10(labels = scales::label_dollar(accuracy = 1)) +
  scale_x_continuous(labels = scales::percent) +
  scale_size_area(max_size = 5, guide = "none") +
  labs(x = "share of graduates who are women", y = "median earnings (log scale)")

Median earnings of recent graduates against the share of women, for 172 majors (grey). Amber points: field averages. The steep red line is the fit ignoring field; the short blue lines are the within-field fits, with a common, shallower slope. Much of the steep overall slope comes from differences between fields.

The steep red line is mostly drawn by a few clusters: engineering (few women, high pay) at the upper left, and education, psychology and health (many women, lower pay) at the lower right. Within each field, the blue lines still slope downward, but more gently. Conditioning on field has closed the fork $S W \leftarrow F \to E$ .

What "controlling for" a variable actually does

A coefficient in a multiple regression has a precise and illuminating meaning, discovered by Ragnar Frisch and Frederick Waugh in 1933. The coefficient on $S W$ in a model that also includes field equals the slope you get by (1) removing from log earnings everything that field predicts, (2) removing from the share of women everything that field predicts, and (3) regressing the first set of leftovers on the second. Multiple regression asks: once I know a major's field, does knowing its share of women tell me anything more about its earnings? You can check this numerically (with ordinary least squares, which matches the posterior mean closely here):

r_earn  <- resid(lm(log_earn ~ field, data = majors))     # earnings not predicted by field
r_women <- resid(lm(ShareWomen ~ field, data = majors))   # share of women not predicted by field
c(by_hand = coef(lm(r_earn ~ r_women))[["r_women"]],
  multiple_regression = coef(lm(log_earn ~ field + ShareWomen, data = majors))[["ShareWomen"]])

output

            by_hand multiple_regression 
            -0.3784             -0.3784

This is also why a regression coefficient can only be as good as the variables you condition on. If field is measured coarsely (16 broad categories), the leftovers in step (2) still contain differences between, say, nursing and health administration, and any confounding they carry remains.

So is the remaining within-field gap a causal effect of women's presence on pay? The data cannot say, and this is where the graph earns its keep. Two very different mechanisms produce exactly this pattern. In the sorting explanation, men and women choose different majors within each field, for many reasons, and the majors men choose happen to lead to better-paid occupations (the dashed $U \to S W$ and $U \to E$ arrows: unmeasured features of each major affect both). In the devaluation explanation, work becomes paid less because women do it (the $S W \to E$ arrow itself). Both predict a negative within-field slope. A cross-section of 172 majors cannot separate them, because it cannot block the unmeasured $U$ .

History: Does Pay Fall When Women Arrive?

Sociologists have argued about sorting versus devaluation for decades, and the debate is a model of how to make progress on a causal question with observational data. In 2009 Asaf Levanon, Paula England and Paul Allison used US census data from 1950 to 2000 to follow the same occupations over time. If devaluation is real, an occupation's pay should fall after women enter it; if sorting explains everything, women's entry should follow, not precede, changes in pay. Using models that compared each occupation with itself across decades (removing stable unmeasured differences between occupations), they found evidence that occupations' pay declined after the share of women rose, consistent with devaluation, though other researchers dispute the size of the effect. The design lesson is general: when a cross-section cannot block an unmeasured confounder, change the data: follow units over time, find a natural experiment, or measure the confounder. Chapters 8 and 13 develop these designs.

Research Practice: State the Estimand Before the Model

Lundberg, Johnson and Stewart (2021) argue that the most important sentence in a quantitative paper is the one that says what quantity, for which population, the analysis estimates: the estimand. "The association between share of women and earnings" is not an estimand; it does not say whether field is held fixed. Better: "the difference in median log earnings between majors in the same broad field whose shares of women differ by ten percentage points, among US college majors in 2010–2012". Writing the estimand first forces you to decide which variables to condition on, and the graph then tells you whether the model identifies it. Then state plainly what the estimate does not show: here, whether women's presence causes lower pay.

Two more lessons from these data. First, the Bayesian R² (the share of variance explained; Gelman et al., 2019) rises from about 0.42 to about 0.65 when field is added, but a high R² is not evidence that a model is causally correct: the model with field explains earnings better and still cannot tell sorting from devaluation. Second, the medians are themselves estimates: 10% of majors have fewer than 15 graduates in the survey sample, and Library Science's median is based on two people. A model that treats every median as exact gives these majors too much influence. Exercise 6 asks you to check whether this changes the conclusion; Chapters 6 and 7 show how to model measurement uncertainty directly.

rbind(`ignoring field` = bayes_R2(m_maj1), `with field` = bayes_R2(m_maj2)) |> round(2)

output

   Estimate Est.Error Q2.5 Q97.5
R2     0.42      0.04 0.33  0.50
R2     0.65      0.03 0.59  0.69

Check Yourself

A newspaper writes: "Within the same field, majors with more women pay about 30% less, so if more women chose a major its pay would fall by up to 30%." What is wrong with the second half?Nothing; the within-field estimate removes all confounding.Conditioning on field closes only the fork through field; unmeasured features of majors (

U

) can still produce the within-field association, so it may not be causal.The 30% is on the log scale and cannot be interpreted as a percentage.

Adding field raised R² from about 0.42 to about 0.65. This shows that

3.9 Good controls and bad controls

The last two sections conditioned on a variable and changed the answer: height turned a total effect into a direct one, and field removed a fork. It is tempting to conclude that adding variables makes estimates more trustworthy, and much applied research behaves as if it did, routinely "controlling for" every available variable. The graph says otherwise. A control variable $Z$ is good or bad depending on where it sits relative to the exposure $X$ and the outcome $Y$ , and there are only a handful of basic positions. The widget simulates each one.

What happens to the estimate of X → Y when you add Z?

Each setting simulates 300 datasets of 100 units from the graph shown, in which the true total effect of $X$ on $Y$ is marked by the dashed line. For each dataset we fit $Y$ on $X$ alone (blue) and $Y$ on $X$ and $Z$ (amber), and plot the histograms of the 300 estimated slopes on $X$ . The slider sets the strength of the arrows involving $Z$ .

where Z sits strength of Z's arrows =

Go through the six settings. Each corresponds to a rule.

Confounder (fork): without $Z$ the estimate is biased; with $Z$ it is centred on the truth. Good control.
Mediator (pipe): without $Z$ the estimate is centred on the total effect; with $Z$ it estimates only the direct part. Bad control for the total effect, good for the direct effect, as in Section 3.7.
Collider: without $Z$ the estimate is fine; with $Z$ it is biased, even though $Z$ is "just another variable". Bad control.
Cause of $Y$ only: both estimates are centred on the truth, but adding $Z$ makes the amber histogram narrower, because $Z$ explains some of the noise in $Y$ . Good control, for precision.
Cause of $X$ only: both are unbiased, but adding $Z$ makes the estimate less precise, because it removes variation in $X$ that was useful. Neutral to harmful.
Effect of $Y$ : conditioning on a consequence of the outcome biases the estimate, toward zero here. Bad control, and a common one: think of controlling for a diagnosis that the outcome itself triggers.

Research Practice: Choose Controls with the Graph, and Leave Their Coefficients Alone

Cinelli, Forney and Pearl (2022) catalogue these cases as "a crash course in good and bad controls", and reviewers increasingly expect authors to justify each covariate from a causal diagram rather than by habit. Two habits to adopt:

Justify every control. List each covariate with its role (confounder, precision variable, or otherwise) in a table or in the DAG figure. Never include post-treatment variables in an analysis of a total effect.
Avoid the "Table 2 fallacy" (Westreich and Greenland, 2013). A model built to estimate the effect of $X$ adjusts for the confounders of $X$ ; it does not, in general, estimate the effects of those confounders. In m_maj2, the field intercepts are not the causal effects of field on earnings, because nothing in the design removes the confounders of field. Report the coefficient you designed the model for, and leave the others in the supplement or out of the discussion.

Project Connection: AI and University Mathematics

A study of whether students' use of AI tools affects their grades in mathematics courses will be judged largely on its choice of controls. Prior attainment (earlier GPA, placement scores) plausibly affects both AI use and grades: a confounder, so adjust for it. Hours spent on homework or assignment completion may be affected by AI use: mediators, so do not adjust for them if the question is the total effect. Attending tutoring may be a collider if both struggling grades and AI enthusiasm lead students to it. Course section or instructor may affect grades only: precision variables, worth including. Draw the graph before collecting the data, so that you collect the confounders, and pre-register the adjustment set (Chapter 1).

Check Yourself

In the widget, for the structure "cause of

Y

only", the amber histogram is narrower than the blue one, and both are centred on the truth. Why does adding

Z

help?

Z

explains part of the variation in

Y

, so the remaining noise is smaller and the slope on

X

is estimated more precisely.

Z

removes confounding of the effect of

X

Z

makes the estimate biased but less variable.

A study of police body cameras (X) on use-of-force complaints (Y) controls for the number of arrests each officer makes (Z), which the cameras may change. In two or three sentences, explain what could go wrong.

3.10 Curves from lines: how fast do children grow?

Many relationships are not straight. Children grow quickly when they are small, slow down, speed up again at puberty, and stop. A linear model can still describe this, because "linear" refers to the parameters, not to the shape. The trick is to replace the single predictor "age" with several fixed functions of age, called basis functions, and let the model choose how much of each to use:

μ_{i} = α + \sum_{k = 1}^{K} w_{k} B_{k} ({age}_{i}) .

Each $B_{k}$ is a smooth bump that is non-zero only over a short range of ages. The weights $w_{k}$ are ordinary regression coefficients, so the model is still linear and fitted exactly like the others. Adding up the bumps, each scaled by its weight, produces a smooth curve. The bumps used in practice are B-splines: pieces of cubic polynomials joined smoothly at points called knots. The widget fits girls' heights from NHANES this way.

Building a growth curve from bumps

Points: heights of 451 girls aged 2 to 19 (NHANES 2009–10). The cubic B-spline basis has the number of knots you choose, evenly spaced in age. The weights are the posterior means under a Normal(0, $τ$ ) prior on each weight (with the noise standard deviation fixed at 6 cm), so $τ$ controls how freely the curve can bend. Thin curves: each basis function times its weight, drawn around the girls' average height. Thick curve: their sum.

interior knots = prior sd of weights τ = cm

With no interior knots, the basis is a single cubic polynomial, which cannot both rise steeply in childhood and level off sharply after puberty. With three or four knots, the curve follows the data well. With fifteen and a generous prior, it starts to wiggle, chasing noise in particular age ranges: that is overfitting. Shrink $τ$ and the wiggles disappear, because a tighter prior on the weights pulls them toward zero and the curve toward smoothness. A prior on spline weights is a smoothness assumption.

That observation leads to the method used in practice. Rather than choosing the number of knots by hand, use many basis functions and estimate how strongly to shrink their weights from the data, by giving $τ$ its own prior. This is a penalised spline, and in brms it is the s() term, borrowed from the mgcv package. The model below fits a separate smooth curve for girls and boys.

kids <- nh |>
  distinct(ID, .keep_all = TRUE) |>
  filter(Age >= 2, Age <= 19) |>
  drop_na(Height, AgeMonths) |>          # age in months was recorded in 2009-10
  mutate(age = AgeMonths / 12)
count(kids, Gender)

m_sp <- brm(Height ~ Gender + s(age, by = Gender, k = 10), data = kids,
            prior = c(prior(normal(130, 30), class = Intercept),
                      prior(normal(0, 50), class = b),
                      prior(exponential(0.2), class = sds),     # how wiggly the curves may be
                      prior(exponential(0.2), class = sigma)),
            seed = 3, refresh = 0, control = list(adapt_delta = 0.95),
            file = "fits/ch3_spline")

output

  Gender   n
1 female 451
2   male 490

The fitted curves are only the beginning. Because we have posterior samples of each whole curve, we can ask questions about features of the curve, with honest uncertainty. When do children grow fastest during puberty? For each posterior sample, compute the curve on a fine grid of ages, take differences to get the growth rate in cm per year, and find the age between 8 and 17 at which it peaks.

ages <- expand_grid(Gender = c("female", "male"), age = seq(2, 19.5, by = 0.1))
curve <- posterior_epred(m_sp, newdata = ages, ndraws = 1000)   # 1000 samples x 352 ages

band <- function(m) tibble(est = colMeans(m), lo = apply(m, 2, quantile, 0.025),
                           hi = apply(m, 2, quantile, 0.975))
heights <- bind_cols(ages, band(curve))
vel <- (curve[, -1] - curve[, -ncol(curve)]) / 0.1               # cm per year
same <- ages$Gender[-1] == ages$Gender[-nrow(ages)]              # don't difference across sexes
speed <- bind_cols(ages[-1, ][same, ], band(vel[, same]))

peak_age <- function(g) {
  cols <- which(ages$Gender[-1] == g & same & ages$age[-1] >= 8 & ages$age[-1] <= 17)
  ages$age[-1][cols][apply(vel[, cols], 1, which.max)]
}
peaks <- tibble(girls = peak_age("female"), boys = peak_age("male"), `boys - girls` = boys - girls)
peaks |> pivot_longer(everything()) |> group_by(name) |>
  summarise(mean = mean(value), lo = quantile(value, 0.025), hi = quantile(value, 0.975))

p1 <- ggplot(heights, aes(age, est, colour = Gender, fill = Gender)) +
  geom_point(data = kids, aes(age, Height), alpha = 0.25, size = 0.7, inherit.aes = FALSE,
             colour = harbour[["muted"]]) +
  geom_ribbon(aes(ymin = lo, ymax = hi), colour = NA, alpha = 0.3) + geom_line() +
  labs(x = "age (years)", y = "height (cm)", colour = NULL, fill = NULL)
p2 <- ggplot(speed, aes(age, est, colour = Gender, fill = Gender)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), colour = NA, alpha = 0.3) + geom_line() +
  coord_cartesian(ylim = c(0, 12)) +
  labs(x = "age (years)", y = "growth (cm per year)", colour = NULL, fill = NULL)
patchwork::wrap_plots(p1, p2) + patchwork::plot_layout(guides = "collect") &
  scale_colour_manual(values = c(female = harbour[["rose"]], male = harbour[["sky"]])) &
  scale_fill_manual(values = c(female = harbour[["rose"]], male = harbour[["sky"]])) &
  theme(legend.position = "top")

output

# A tibble: 3 x 4
  name          mean     lo    hi
  <chr>        <dbl>  <dbl> <dbl>
1 boys         12.0   9.3    13.2
2 boys - girls  2.73 -0.303   4.8
3 girls         9.26  8      10.5

Left: heights of 941 children (2009–10) with the posterior mean growth curve and 95% credible band for girls and boys. Right: growth velocity (cm per year) computed from the same posterior samples. The pubertal spurt is clear for boys and faint for girls, because a cross-section averages over children whose spurts come at different ages.

The left panel shows the familiar picture: girls and boys grow at similar rates until about age 10; girls then stop growing at around 15 and boys continue to 17 or 18, ending about 14 cm taller. The right panel shows something growth charts hide: growth slows through childhood and then speeds up again at puberty, clearly for boys and only faintly for girls. The posterior puts boys' fastest growth after age 8 at about 12 years (95% interval 9.3 to 13.2). For girls it cannot even separate a pubertal peak from the slow decline of childhood growth: the fastest growth after age 8 could be anywhere from age 8 to about 10.5.

Paediatric studies that measure the same children every few months find something much sharper: peak growth of 8 to 10 cm per year (our curves peak at about 7), at about 11.5 years for girls and 13.5 for boys. Is the model wrong? No: it answers a different question. Every child in NHANES was measured once, so our curve is the average height of children at each age, and its slope is the rate at which that average changes. Children enter puberty at different ages, so at any age the average mixes children in their spurt with children before and after it, and the population curve's spurt is lower and more spread out than any individual child's. With only about 25 girls per year of age, the faint hump is also hard to locate. To estimate a typical child's peak growth you need a different estimand and different data: repeated measurements of the same children, modelled with a curve for each child (Chapters 6 and 8). The cross-section was not useless. It said exactly as much as it could, and the posterior said it honestly, with intervals wide enough to show the ambiguity.

Why penalised splines are hierarchical models

In s(age), the wiggly part of the curve has weights $w_{k} \sim Normal (0, τ)$ , with $τ$ (reported by brms as sds) estimated from the data. A parameter that sets the spread of a whole family of other parameters is a hyperparameter, and a model with one is a multilevel (hierarchical) model. The data decide how much the curve may bend, just as in Chapter 6 the data decide how much tadpole tanks or survey subgroups should resemble each other. Penalised splines are hierarchical models in disguise, and the shrinkage you saw in the widget is the same partial pooling.

Check Yourself

In the spline widget with 15 knots, reducing the prior standard deviation

τ

of the weights makes the curve smoother. Why?It reduces the number of basis functions.It pulls the weights toward zero, so neighbouring bumps cannot take very different values, and sharp bends become improbable a priori.It makes the data less noisy.

Research Practice: A Complete Results Paragraph

Here is what a publishable results paragraph for the analysis of Section 3.7 looks like. It states the estimands, the model, the priors, the key contrasts with intervals, the checks, and the caveat.

"We estimated the total effect of sex on adult body weight, and the controlled direct effect at equal height, using 4,609 adults aged 20 and over in NHANES 2009–2012 (one record per person). Weight was modelled as normal with sex-specific means and, for the direct effect, a common linear term in height centred at 170 cm, with weakly informative priors (Normal(80, 20) for means in kg, Normal(0, 2) for the slope in kg per cm, Exponential(0.1) for the residual standard deviation). Men weighed 12.0 kg more than women on average (95% credible interval 10.8 to 13.2). At equal height the difference was −0.8 kg (−2.3 to 0.7), so the sex difference in weight is accounted for by height. A posterior predictive check showed that the normal model under-predicts the right tail of weight; a log-scale model gave the same substantive conclusion (Supplement S2). The direct effect assumes no unmeasured common causes of height and weight; Supplement S3 reports a sensitivity analysis. Code and data are available at [repository]."

Key Takeaways: Lecture 6

Use index variables for categories, and compute every comparison you care about as a contrast of posterior samples.
Conditioning on a mediator turns a total effect into a direct effect. Decide which you want before choosing the model.
Conditioning on a common cause removes confounding by that cause only; unmeasured confounders remain, and no fit statistic detects them.
Controls are good or bad according to their position in the graph: confounders and causes of the outcome help; mediators, colliders and effects of the outcome can hurt. Do not interpret the coefficients of controls.
Splines make curves from linear models; a prior on their weights is a smoothness assumption, and any feature of the curve has a posterior.

3.11 Exercises

Each exercise has a hint and an answer button.

Exercises

The cost of duplicates. Fit the model m_h to the women's heights with the duplicated rows kept (use filter(nh, Age >= 20, Gender == "female") |> drop_na(Height)). Compare the posterior standard deviation of the mean with that of m_h. Is the ratio what you expected?

Use update(m_h, newdata = ...) and compare fixef() outputs. The expected ratio is the square root of the ratio of the numbers of rows.

Keeping the duplicates gives 3,658 rows for 2,373 women. The posterior standard deviation of $μ$ falls from 0.151 to 0.121 cm, a ratio of 1.25, almost exactly the predicted $\sqrt{3658 / 2373} = 1.24$ : the model is 25% more confident with no new information. The posterior mean also moves, from 161.58 to 162.05 cm. The duplicates are not random: they encode the survey weights, so the duplicated data estimate the population average height (weighting by how many Americans each woman represents), while the de-duplicated data over-represent the oversampled groups. Neither is simply right. The duplicated version has the right weights and the wrong uncertainty; a proper survey analysis gets both right.
Men's heights. Fit the Gaussian model of Section 3.2 to the men, with the same priors. Report the posterior mean and 95% interval of their average height, and the posterior of the difference between men's and women's average heights.

Fit Height ~ 0 + Gender to adults and compute the contrast, or fit two separate models and subtract their draws.

Men's average height is 175.35 cm (95% interval 175.05 to 175.66), and women's 161.58 cm (161.29 to 161.88). The contrast, men minus women, has posterior mean 13.77 cm (13.34 to 14.21).
Predicting a person. Using m_wh, give the posterior mean and 95% interval for (a) the average weight of adults who are 185 cm tall, and (b) the weight of one new 185 cm adult. Explain in one sentence why (b) is so much wider than (a).

fitted(m_wh, newdata = tibble(h = 15)) and predict(m_wh, newdata = tibble(h = 15)).

(a) The average weight of 185 cm adults: 97.1 kg (95% interval 96.1 to 98.2). (b) One new 185 cm adult: about 97 kg (about 60 to 135). The predictive interval is about 75 kg wide, the credible interval about 2 kg, because a single person's weight varies around the average for their height with a standard deviation of about 19 kg, and more data cannot remove that variation.
Forgetting to centre. Refit m_wh with Weight ~ 1 + Height (height not centred) and the same priors for the slope and $σ$ , and a Normal(0, 100) prior for the intercept. What is the posterior mean of the intercept, and what does it mean? Compute the posterior correlation between the intercept and the slope, and compare it with the same correlation in m_wh.

as_draws_df(fit) |> with(cor(b_Intercept, b_Height)).

The intercept is about −69.9 kg (−79.0 to −60.8): the expected weight of a person of height zero, an impossible extrapolation with no meaning. The posterior correlation between intercept and slope is −0.998: almost any increase in the slope must be paid for by a decrease in the intercept to keep the line through the data. In the centred model, the correlation is only 0.17 (not exactly zero, because the centre, 170 cm, is a little above the sample mean of 168 cm). The slope itself is unchanged (0.90), but the intercept is uninterpretable and its prior impossible to set sensibly. (brms centres predictors internally for sampling, which is why this model still fits quickly; Stan without that trick would struggle.)
Checking the log model. Repeat the posterior predictive check of Section 3.5 for m_log, comparing the skewness of observed log weights with the skewness of simulated log weights. Does the log model pass this check?

Replace m_wh by m_log and adults$Weight by adults$log_w in the ppc_wh chunk.

The observed skewness of log weight is 0.18; simulated datasets have skewness near 0 (mean −0.01), and none of the 200 is as skewed as the data. So, strictly, the log model also fails this check. But the misfit is far smaller than on the kilogram scale (1.1 against 0), and with 4,609 people even a small discrepancy is detectable. Whether it matters depends on the use: for the exponent, which depends on the mean, it does not; for predicting the heaviest few percent of people, a model with a heavier right tail (a skew-normal or Student-t likelihood, Chapter 4) would be better.
Unreliable medians. Refit m_maj2 using only majors with at least 30 graduates in the survey sample (Sample_size >= 30). How many majors remain, and how does the within-field percentage difference change? What does this suggest about the influence of the small majors?

Filter, then update(m_maj2, newdata = ...), and apply the pct() function.

140 of the 172 majors have at least 30 graduates in the sample (all 16 fields remain). The within-field difference becomes −25% (95% interval −37% to −11%), against −31% (−40% to −19%) using all majors. The small majors pushed the estimate further from zero, but the conclusion that majors with more women pay less within the same field survives. Dropping data by an arbitrary threshold is crude; Chapter 7 shows how to weight each median by its precision instead.
Is age a confounder?. Age affects adult weight (people tend to gain weight into middle age) and also height (older adults are slightly shorter, partly because they were born earlier and partly because of spinal compression). Draw the graph with $A$ (age), $H$ and $W$ . Is age a confounder of the effect of height on weight? Fit Weight ~ 1 + h + Age and compare the slope on height with m_wh.

If $A \to H$ and $A \to W$ , there is a back-door path $H \leftarrow A \to W$ .

With $A \to H$ and $A \to W$ , age opens a back-door path $H \leftarrow A \to W$ , so it is a confounder of the effect of height on weight and should be adjusted for. In practice it matters little here: the slope on height moves from 0.903 to 0.915 kg per cm, because age is only weakly related to height among adults (correlation −0.16) and adds only 0.04 kg per year to weight (0.01 to 0.07). A confounder needs to be associated with both the exposure and the outcome to cause much bias.
Growth spurts. Using the posterior samples in the growth chunk, compute the posterior distribution of the peak growth velocity (cm per year at the peak, between ages 8 and 17) for girls and boys. Which sex grows faster at its peak, and by how much?

Replace which.max by max in the peak_age() function.

The posterior peak velocity is about 7.0 cm per year for girls (6.3 to 7.9) and 7.3 for boys (6.6 to 8.3), with a difference of 0.3 (−0.9 to 1.5): no clear difference. Longitudinal studies of individual children find peak velocities of 8 to 10 cm per year, with boys' clearly higher. The cross-sectional estimates are lower and similar, which is exactly the smearing described in Section 3.10: averaging growth spurts that happen at different ages for different children flattens them.

Chapter 3Bayesian linear models

Think About It: Before You Start

Lecture 5 — Lines, and the uncertainty around them

3.1 Why the bell curve keeps appearing

History: Galton's Board and the Origin of "Regression"

Check Yourself

3.2 A model of one variable: how tall are American women?

Research Practice: Know How Your Data Were Made

Check Yourself

Definition: The Gaussian Model

The posterior, computed by brute force

The same model in brms

3.3 Drawing a line: weight and height

Definition: The Linear Model

3.4 Two kinds of uncertainty

Research Practice: Say Which Interval You Mean

Check Yourself

3.5 Checking the model, and the shape of the human body

History: Quetelet's Index and the Exponent of Two

Research Practice: Report the Checks You Ran

Check Yourself

Key Takeaways: Lecture 5

Lecture 6 — What a coefficient means

3.6 Categories as predictors

Research Practice: Compare Groups with Contrasts, Not by Eye

3.7 Total and direct effects

Warning: Direct Effects Need More Than a Graph Without Back Doors

Check Yourself

3.8 Confounding in the wild: majors and money

History: Does Pay Fall When Women Arrive?

Research Practice: State the Estimand Before the Model

Check Yourself

3.9 Good controls and bad controls

Research Practice: Choose Controls with the Graph, and Leave Their Coefficients Alone

Project Connection: AI and University Mathematics

Check Yourself

3.10 Curves from lines: how fast do children grow?

Check Yourself

Research Practice: A Complete Results Paragraph

Key Takeaways: Lecture 6

3.11 Exercises

Exercises

Further reading and references

Further Reading

Chapter 3
Bayesian linear models