Chapter 6
Multilevel models

Lectures 11–12 · partial pooling, empirical Bayes, varying slopes, the geometry of hierarchies, and ambulance response times across New York

Data almost always come in groups: patients in hospitals, students in schools, neighbourhoods in boroughs, measurements of the same country over many years. The groups differ, and they differ in size. A few have hundreds of observations and some have one. How should we estimate each group's average? We could treat every group as unrelated to every other, trusting each group's own data entirely. We could ignore the groups and use one overall average for all of them. Or we could do something in between: estimate each group's average, but let the groups inform each other, so that groups with little data borrow strength from the rest.

The third option is the multilevel model (also called a hierarchical or mixed-effects model), and it is probably the single most useful technique in applied Bayesian statistics. It treats the group averages as themselves drawn from a population, learns from the data how much the groups differ, and pulls each group's estimate toward the population average by exactly the amount its own data warrant. Groups with lots of data barely move; groups with little data move a lot. This shrinkage sounds like a compromise, but it is not: shrunken estimates are, on average, more accurate than the raw group averages. That fact so surprised statisticians when it was proved in the 1950s that it became known as a paradox.

• Lecture 11 is about pooling. We will ask which country grows the best coffee, using ratings of 1,309 coffees from 36 countries, see why the raw rankings mislead, fit a multilevel model, check on held-out data that its estimates predict better, and meet empirical Bayes: the shortcut that estimates the population from the data and the history of Stein's paradox.
• Lecture 12 extends the idea. Groups can differ in their slopes as well as their averages: we will model how carbon emissions per person changed in 164 countries from 2000 to 2022, and ask whether richer countries cut them faster. Then we confront the computational price of hierarchies, the funnel, and the reparameterisation that removes it, and end with a case study from New York City's ambulance dispatch records: response times to cardiac arrests across nearly 200 ZIP codes.

Lecture 11 — Pooling

6.1 Many small groups: which country grows the best coffee?

The Coffee Quality Institute trains and certifies professional tasters ("Q graders") who score green coffee beans on a 100-point scale, adding up marks for aroma, flavour, aftertaste, acidity, body, balance and other attributes. In 2018 the data scientist James LeDoux collected the institute's public reviews of more than 1,300 coffees, and the data were later shared through the TidyTuesday project. We keep the 1,309 Arabica coffees that have a country of origin and a valid score (one record has a score of zero) and ask a simple question: which country produces the best coffee?

library(tidyverse)
library(brms)
library(posterior)
options(mc.cores = parallel::detectCores())

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

url <- "https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2020/2020-07-07/coffee_ratings.csv"
coffee <- read_csv(url, show_col_types = FALSE) |>
  filter(species == "Arabica", total_cup_points > 50, !is.na(country_of_origin)) |>
  transmute(country = country_of_origin, score = total_cup_points, method = processing_method)
nrow(coffee); n_distinct(coffee$country)

by_country <- coffee |> group_by(country) |> summarise(n = n(), mean = mean(score))
by_country |> arrange(-mean) |> head(6)          # the raw "top six"
by_country |> arrange(mean) |> head(4)           # and the bottom four

[1] 1309
[1] 36
# A tibble: 6 x 3
  country              n  mean
  <chr>            <int> <dbl>
1 United States        8  86.0
2 Papua New Guinea     1  85.8
3 Ethiopia            44  85.5
4 Japan                1  84.7
5 Kenya               25  84.3
6 Uganda              26  84.1
# A tibble: 4 x 3
  country           n  mean
  <chr>         <int> <dbl>
1 India             1  76.8
2 Haiti             6  77.2
3 Cote d?Ivoire     1  79.3
4 Nicaragua        26  80.5

The raw ranking puts the United States first (8 coffees, average 86.0), then Papua New Guinea, represented by a single coffee, then Ethiopia (44 coffees averaging 85.5) and Japan (again a single coffee). At the bottom are India (one coffee), Haiti (six) and Côte d'Ivoire (one). Would you bet that Papua New Guinea grows better coffee than Ethiopia on the strength of one bag? The top and bottom of a ranking of raw averages are dominated by small groups, for the simplest of reasons: the average of a few observations is noisy, and noise pushes some small groups far above the truth and others far below it.

6.2 Three ways to pool

Write $y_{ij}$ for the score of coffee $i$ from country $j$. Three models are possible, and they differ only in what they assume about the country averages ${\alpha }_j$.

• Complete pooling assumes all countries are the same: $y_{ij}\sim \text{Normal}(\alpha ,\sigma )$ with a single $\alpha$. It ignores real differences between countries.
• No pooling assumes the countries have nothing to do with each other: each ${\alpha }_j$ gets its own fixed prior, and its estimate is essentially the country's own average. It treats one coffee from Japan as if it told us all about Japanese coffee.
• Partial pooling assumes the countries are different but related: their averages are drawn from a common population, ${\alpha }_j\sim \text{Normal}(\mu ,\tau )$, where $\mu$ is the average of the country averages and $\tau$ how much they differ. Crucially, $\mu$ and $\tau$ are parameters, estimated from the data along with everything else.

The partial-pooling model is a multilevel model because it has two levels of variation: coffees vary around their country's average (with standard deviation $\sigma$), and countries vary around the population average (with standard deviation $\tau$). The prior for each ${\alpha }_j$ is not fixed in advance; it is learned from the other countries. In full:

When $\sigma$, $\mu$ and $\tau$ are known, each country's posterior mean has a simple and revealing form:

The estimate is a weighted average of the country's own mean and the population mean. The weight on the country's own data, $w_j$, compares how much countries genuinely differ (${\tau }^2$) with how noisy the country's own average is (${\sigma }^2/n_j$). With many coffees, ${\sigma }^2/n_j$ is small, $w_j$ is near 1, and the country keeps its own average. With one coffee, $w_j$ may be well below one-half and the estimate is pulled most of the way to the population mean. Complete pooling is the case $\tau =0$; no pooling is $\tau =\infty$. The widget lets you set $\tau$ by hand.

At $\tau$ near zero, every country gets the overall average: no country is better than any other. At very large $\tau$, every country keeps its raw average, and the one-coffee countries sit at the extremes. In between, the large countries barely move while the small ones slide toward the centre, and the ranking changes: at the estimated $\tau$, Ethiopia, the United States and Kenya are the top three, Papua New Guinea drops from second to fifth, and Japan falls into the pack. Watch the held-out error as you move the slider. It is large at both ends and smallest somewhere in the middle, near the value of $\tau$ that the multilevel model estimates from the data.

6.3 Fitting a multilevel model

In brms, a varying intercept is written (1 | country): "an intercept that varies by country". The standard deviation of the country intercepts, $\tau$, is a parameter of class sd.

m_coffee <- brm(score ~ 1 + (1 | country), data = coffee,
                prior = c(prior(normal(82, 5), class = Intercept),
                          prior(exponential(1), class = sd),          # tau: spread between countries
                          prior(exponential(0.5), class = sigma)),    # spread within countries
                seed = 6, refresh = 0, file = "fits/ch6_coffee")
summary(m_coffee)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: score ~ 1 + (1 | country) 
   Data: coffee (Number of observations: 1309) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Group-Level Effects: 
~country (Number of levels: 36) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     1.44      0.24     1.02     1.97 1.00      927     1502

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    82.31      0.29    81.73    82.89 1.01      521     1091

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     2.44      0.05     2.35     2.54 1.00     4486     3296

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 summary has the familiar Intercept (now $\mu$, about 82.3 points) and sigma (2.4 points of variation between coffees from the same country), and a new block, Group-Level Effects (called Multilevel Hyperparameters in recent versions of brms), reporting sd(Intercept): $\tau$, the standard deviation of the true country averages, about 1.4 points with a 95% interval of roughly 1.0 to 2.0. Countries differ, but by less than individual coffees within a country do. A common summary is the intraclass correlation, ${\tau }^2/({\tau }^2+{\sigma }^2)$: the share of the variation in scores that lies between countries, here about a quarter.

est <- coef(m_coffee)$country[, , "Intercept"] |> as_tibble(rownames = "country") |>
  left_join(by_country, by = "country")
mu_hat <- fixef(m_coffee)["Intercept", "Estimate"]
ggplot(est, aes(n)) +
  geom_hline(yintercept = mu_hat, linetype = "dashed", colour = harbour[["muted"]]) +
  geom_segment(aes(xend = n, y = mean, yend = Estimate), colour = harbour[["muted"]], alpha = 0.6) +
  geom_point(aes(y = mean), shape = 1, colour = harbour[["rose"]], size = 2) +
  geom_point(aes(y = Estimate), colour = harbour[["sky"]], size = 2) +
  scale_x_log10() +
  labs(x = "number of coffees (log scale)", y = "average score")

The figure is the whole chapter in one picture. On the right, Mexico, Colombia and Guatemala, with well over a hundred coffees each, are estimated at essentially their own averages. On the left, the one-coffee countries are pulled about three-quarters of the way to the population mean (the weight from Check 6.1). Haiti, whose six coffees averaged 77.2, the lowest of any country, is moved up by more than a point; the United States, at the top with eight, is moved down by almost one.

A ranking should carry its uncertainty. With posterior draws of every country's average, we can compute, for each draw, which country comes first, and report the share of draws in which each country is best: the posterior probability of being the best.

draws_c <- coef(m_coffee, summary = FALSE)$country[, , "Intercept"]   # draws x countries
best <- colnames(draws_c)[apply(draws_c, 1, which.max)]
tibble(country = best) |> count(country, name = "draws") |>
  mutate(`Pr(best)` = draws / nrow(draws_c)) |> arrange(-draws) |> head(6)

# A tibble: 6 x 3
  country          draws `Pr(best)`
  <chr>            <int>      <dbl>
1 Ethiopia          2247     0.562 
2 United States     1224     0.306 
3 Papua New Guinea   181     0.0452
4 Japan              107     0.0268
5 Ecuador             63     0.0158
6 Rwanda              44     0.011 

Ethiopia is the best coffee-producing country in a little over half the posterior draws and the United States in about 30%; no other country reaches 5%. Papua New Guinea and Japan, second and fourth on raw averages, come first in under 5% and 3% of draws. "Ethiopia probably, possibly the United States (on the evidence of eight coffees, mostly from Hawaii), and we cannot tell beyond that" is a far more honest conclusion than a league table.

6.4 Does pooling predict better?

The claim that shrunken estimates are more accurate is testable. Split every country's coffees into two random halves. Estimate each country's average from the first half, in three ways (no pooling, complete pooling and partial pooling), and see which estimates come closest to the average of the other half.

set.seed(6)
coffee <- coffee |> group_by(country) |>
  mutate(half = sample(rep(1:2, length.out = n()))) |> ungroup()
first  <- filter(coffee, half == 1)
second <- coffee |> filter(half == 2) |> group_by(country) |> summarise(n2 = n(), target = mean(score))

m_half <- update(m_coffee, newdata = first, seed = 6, refresh = 0)
predictions <- first |> group_by(country) |> summarise(n1 = n(), `no pooling` = mean(score)) |>
  mutate(`complete pooling` = fixef(m_half)["Intercept", "Estimate"],
         `partial pooling` = coef(m_half)$country[country, "Estimate", "Intercept"])
predictions |> inner_join(second, by = "country") |>
  pivot_longer(c(`no pooling`, `complete pooling`, `partial pooling`), names_to = "method") |>
  group_by(method) |> summarise(countries = n(), mse = mean((value - target)^2)) |> arrange(mse)

# A tibble: 3 x 3
  method           countries   mse
  <chr>                <int> <dbl>
1 no pooling              28  1.49
2 partial pooling         28  1.63
3 complete pooling        28  3.38

A surprise: in this split, the raw averages (no pooling) predicted the other half slightly better than the multilevel estimates, although both beat complete pooling by a wide margin. Is shrinkage a fraud? No. Stein's result, which we meet in the next section, is about expected error, averaged over the noise; any single split is one draw of that noise, and here a couple of small countries happened to have first-half averages that their second halves agreed with. The honest test repeats the split many times. Refitting brms two hundred times would take a while, so the loop below uses lme4::lmer, the fast non-Bayesian fit of the same varying-intercept model, which with 36 countries gives nearly identical estimates.

library(lme4)
one_split <- function(seed) {
  set.seed(seed)
  d <- coffee |> group_by(country) |> mutate(half = sample(rep(1:2, length.out = n()))) |> ungroup()
  fit <- lmer(score ~ 1 + (1 | country), data = filter(d, half == 1))
  est <- coef(fit)$country |> as_tibble(rownames = "country") |> rename(partial = `(Intercept)`)
  d |> filter(half == 1) |> group_by(country) |> summarise(n1 = n(), raw = mean(score)) |>
    mutate(complete = fixef(fit)[[1]]) |> left_join(est, by = "country") |>
    inner_join(d |> filter(half == 2) |> group_by(country) |> summarise(target = mean(score)), by = "country") |>
    mutate(split = seed)
}
splits <- map_dfr(1:200, one_split)
mse <- splits |> group_by(split) |>
  summarise(across(c(raw, complete, partial), ~ mean((.x - target)^2)))
mse |> summarise(`Pr(partial beats raw)` = mean(partial < raw), across(c(raw, complete, partial), mean))

# A tibble: 1 x 4
  `Pr(partial beats raw)`   raw complete partial
                    <dbl> <dbl>    <dbl>   <dbl>
1                    0.75  2.52     3.38    2.01

Averaged over 200 random splits, partial pooling has the lowest error: about a fifth lower than the raw averages and less than two-thirds that of complete pooling, and it beats the raw averages in roughly three splits out of four. Neither extreme is a good idea; borrowing strength by the right amount is. And note what "the right amount" was: no one chose it. The model estimated $\tau$ from the data, and $\tau$ determined how much each country borrowed. The widget in Section 6.2 reports this same average over 50 of the splits, so you can watch the error change with $\tau$.

6.5 Empirical Bayes

The multilevel model treats $\mu$ and $\tau$ as unknown parameters with priors, and integrates over them. An older and simpler approach, called empirical Bayes, instead estimates them from the data, plugs the estimates in, and then applies Bayes' rule to each group with that estimated prior. The estimates come from the marginal likelihood: the probability of the observed group averages with the group effects integrated out. For our model, each country's observed average is normally distributed around $\mu$ with variance ${\tau }^2+{\sigma }^2/n_j$ (true variation plus sampling noise), so

This is also called type-II maximum likelihood, or evidence maximisation. It is a few lines of R.

sigma_w <- coffee |> group_by(country) |> mutate(dev = score - mean(score)) |> pull(dev) |> sd()
marg_ll <- function(par, d) {                    # log marginal likelihood of the country averages
  mu <- par[1]; tau <- exp(par[2])
  sum(dnorm(d$mean, mu, sqrt(tau^2 + sigma_w^2 / d$n), log = TRUE))
}
fit_eb <- optim(c(82, 0), marg_ll, d = by_country, control = list(fnscale = -1))
c(sigma_within = sigma_w, mu_hat = fit_eb$par[1], tau_hat = exp(fit_eb$par[2]))
posterior_summary(m_coffee, variable = c("b_Intercept", "sd_country__Intercept"))   # full Bayes

sigma_within       mu_hat      tau_hat 
       2.407       82.305        1.412 
                      Estimate Est.Error   Q2.5  Q97.5
b_Intercept             82.307    0.2942 81.733 82.890
sd_country__Intercept    1.443    0.2408  1.024  1.966

Empirical Bayes estimates $\widehat{\tau }\approx 1.4$ points, the same as the multilevel model's posterior mean, and its shrunken country averages are almost identical to the full Bayesian ones. With 36 countries, the data pin down $\tau$ well, and plugging in a single value loses little. Empirical Bayes is fast, transparent and widely used in genomics (estimating thousands of gene-level variances at once), in sports and insurance analytics, and in disease mapping.

Empirical Bayes has one weakness, and it matters most when groups are few. By plugging in $\widehat{\tau }$ as if it were known, it ignores the uncertainty in $\tau$, so its intervals for the group effects are too narrow. When the data pin $\tau$ down, that hardly matters. When they do not, a single number hides the fact that the data are consistent both with "these groups are identical" and with "they differ a lot"; and with few groups the maximum-likelihood estimate is often exactly zero, which forces complete pooling with complete confidence. The widget shows the marginal likelihood for $\tau$ with all 36 countries, and with only the eight countries represented by a single coffee.

With all 36 countries, the likelihood is a well-defined hump around 1.4. With only the eight one-coffee countries it still peaks near 1.4, but it is almost flat from zero to beyond 3: eight noisy scores cannot tell "these countries are identical" from "they differ by three points". Here is what the two approaches conclude from those eight countries alone.

one <- by_country |> filter(n == 1) |> mutate(se = sigma_w)        # eight countries, one coffee each
fit_eb8 <- optim(c(82, 0), marg_ll, d = one, control = list(fnscale = -1))
c(mu_hat = fit_eb8$par[1], tau_hat = exp(fit_eb8$par[2]))

m_one <- brm(mean | se(se) ~ 1 + (1 | country), data = one,            # known measurement error
             prior = c(prior(normal(82, 5), class = Intercept), prior(exponential(1), class = sd)),
             seed = 6, refresh = 0, control = list(adapt_delta = 0.99), file = "fits/ch6_one")
posterior_summary(m_one, variable = "sd_country__Intercept")
coef(m_one)$country[, , "Intercept"] |> round(2)

 mu_hat tau_hat 
 81.958   1.369 
                      Estimate Est.Error    Q2.5 Q97.5
sd_country__Intercept   0.9511    0.8266 0.02546  3.03
                 Estimate Est.Error  Q2.5 Q97.5
Cote d?Ivoire       81.54      1.34 78.55 83.92
Ecuador             82.25      1.27 80.01 85.17
India               81.13      1.56 77.36 83.63
Japan               82.38      1.34 80.04 85.53
Mauritius           81.70      1.25 78.95 83.98
Papua New Guinea    82.57      1.39 80.23 86.03
Rwanda              82.10      1.26 79.76 84.72
Zambia              81.97      1.23 79.58 84.50

Empirical Bayes reports $\widehat{\tau }\approx 1.4$ and proceeds as if that were known. The full Bayesian 95% interval for $\tau$ runs from almost zero to about 3, with posterior mean about 0.95 (the exponential prior pulls gently toward smaller values). The eight country estimates are shrunk hard toward about 82, each with an interval 5 to 6 points wide that averages over every plausible $\tau$. With a slightly different eight countries the empirical Bayes estimate could as easily have been 0 or 2.5, and every downstream interval would have inherited that accident. The formula mean | se(se) tells brms that each observation is itself an estimate with a known standard error: the syntax of a Bayesian meta-analysis, which is exactly this model applied to the results of several studies.

Lecture 12 — Slopes, levels and geometry

6.6 Varying slopes: carbon emissions in 164 countries

Groups can differ in more than their averages. When the same units are measured repeatedly over time, each unit has its own trend, and the trends differ too. Our example is a question at the centre of climate policy. Between 2000 and 2022, did emissions of carbon dioxide per person fall in rich countries and rise in poorer ones, and how much do countries at the same income level differ? The data come from Our World in Data, which compiles territorial (production-based) emissions from the Global Carbon Project (Friedlingstein et al., 2023) and national accounts.

url_co2 <- "https://raw.githubusercontent.com/owid/co2-data/master/owid-co2-data.csv"
owid <- read_csv(url_co2, show_col_types = FALSE) |>
  filter(nchar(iso_code) == 3, !startsWith(iso_code, "OWID"), between(year, 2000, 2022))  # countries only
gdp2000 <- owid |> filter(year == 2000, !is.na(gdp)) |>
  transmute(country, log_gdp = log(gdp / population))                                     # GDP per person in 2000
co2 <- owid |> filter(co2_per_capita > 0) |>
  inner_join(gdp2000, by = "country") |>
  group_by(country) |> filter(n() == 23) |> ungroup() |>                                  # complete series only
  mutate(y = log(co2_per_capita),                 # log tonnes of CO2 per person
         t = (year - 2010) / 10,                  # decades from 2010
         g = log_gdp - mean(gdp2000$log_gdp))     # centred log GDP per person in 2000
c(countries = n_distinct(co2$country), rows = nrow(co2))

countries      rows 
      164      3772 

Each country $j$ now gets its own intercept $a_j$ (its log emissions per person in 2010) and its own slope $b_j$ (the change in log emissions per decade; a slope of −0.1 means about a 10% fall per decade). Both depend on the country's income in 2000, and each country deviates from that dependence by its own amount:

Two things are new. First, log GDP per person is a group-level predictor: it varies between countries, not within them, and it explains part of the differences in both intercepts and slopes. Second, the country deviations in intercept and slope are drawn jointly from a two-dimensional normal distribution with a correlation $\rho$, so the model can learn, for example, whether countries that emit more than expected also tend to cut faster. The prior for a correlation matrix is the LKJ distribution (Lewandowski, Kurowicka and Joe, 2009); LKJ(2) gently favours weak correlations. In brms, (1 + t | country) asks for a correlated varying intercept and slope, and t * g gives the population slope, the effect of income and their interaction (income's effect on the slope).

m_co2 <- brm(y ~ 1 + t * g + (1 + t | country), data = co2,
             prior = c(prior(normal(0, 2), class = Intercept),
                       prior(normal(0, 1), class = b),
                       prior(exponential(1), class = sd),
                       prior(lkj(2), class = cor),
                       prior(exponential(5), class = sigma)),
             seed = 6, refresh = 0, file = "fits/ch6_co2")
fixef(m_co2) |> round(3)
posterior_summary(m_co2, variable = c("sd_country__Intercept", "sd_country__t",
                                      "cor_country__Intercept__t", "sigma")) |> round(3)

          Estimate Est.Error   Q2.5  Q97.5
Intercept    0.672     0.055  0.570  0.789
t            0.105     0.020  0.066  0.144
g            1.205     0.045  1.115  1.291
t:g         -0.163     0.015 -0.192 -0.132
                          Estimate Est.Error   Q2.5 Q97.5
sd_country__Intercept        0.689     0.037  0.623 0.764
sd_country__t                0.235     0.014  0.211 0.264
cor_country__Intercept__t   -0.039     0.075 -0.187 0.106
sigma                        0.120     0.002  0.117 0.123

The population-level results answer the policy question. A country with average income in 2000 emitted about $e^{0.67}\approx 2.0$ tonnes of CO₂ per person in 2010, and its emissions per person were rising, by about 11% per decade. Each unit of log GDP per person (a factor of 2.7 in income) multiplied emissions per person by about $e^{1.21}\approx 3.4$, and lowered the decade-on-decade slope by about 0.16. So emissions per person were rising in poorer countries and falling in the richest: for a country about 2.7 times richer than the average (one log unit), the slope is roughly $0.11-0.16\approx -0.06$ per decade.

But countries at the same income level vary enormously around those trends. The standard deviation of country slopes, ${\tau }_b$, is about 0.23 per decade: larger than the effect of income itself. And once income is in the model, the correlation between a country's intercept and slope deviations is close to zero. Countries that emit more than their income predicts are not systematically cutting faster or slower; income explained the correlation one would see without it (Exercise 6).

show <- c("China", "India", "United States", "United Kingdom", "Brazil", "Nigeria")
nd <- co2 |> filter(country %in% show)
lines_c <- nd |> mutate(fit = fitted(m_co2, newdata = nd)[, "Estimate"],
                        pop = fitted(m_co2, newdata = nd, re_formula = NA)[, "Estimate"])
p1 <- ggplot(lines_c, aes(year, y)) +
  geom_point(colour = harbour[["muted"]], size = 0.8) +
  geom_line(aes(y = fit), colour = harbour[["sky"]], linewidth = 0.9) +
  geom_line(aes(y = pop), colour = harbour[["amber"]], linetype = "dashed") +
  facet_wrap(~ country, ncol = 3) + scale_x_continuous(breaks = c(2005, 2015)) +
  labs(x = NULL, y = "log tonnes CO2 per person")
fe <- fixef(m_co2)[, "Estimate"]
re_t <- ranef(m_co2)$country[, "Estimate", "t"]                      # country deviations in slope
slopes <- co2 |> distinct(country, g) |>
  mutate(slope = fe["t"] + fe["t:g"] * g + re_t[country])
p2 <- ggplot(slopes, aes(g, slope)) +
  geom_hline(yintercept = 0, colour = harbour[["muted"]]) +
  geom_abline(intercept = fixef(m_co2)["t", 1], slope = fixef(m_co2)["t:g", 1], colour = harbour[["rose"]]) +
  geom_point(colour = harbour[["sky"]], alpha = 0.7, size = 1.3) +
  labs(x = "log GDP per person, 2000 (centred)", y = "slope (log change per decade)")
patchwork::wrap_plots(p1, p2, widths = c(1.5, 1))

The left panel shows partial pooling working on slopes. Each country's solid line uses its own 23 years of data, which dominate here, so the lines follow the points; the dashed line is what the model would predict from income alone, for a new country it had never seen. China's emissions per person rose far faster than its income in 2000 predicted; Nigeria's fell although its income predicted a rise; the United Kingdom fell faster than predicted. With 23 observations per country, shrinkage is mild. With a few observations per country (a new survey of a few hundred firms over three years, say), the same model would pull the noisy country slopes strongly toward the income-predicted line.

Try four years for China, the United Kingdom and Nigeria. The multilevel line beats the country's own line on the unseen years for the United Kingdom, but not for China, whose emissions rose far faster than countries of its income usually did, or for Nigeria. Pooling helps on average, not in every case. Across all 164 countries, with four years observed, the multilevel trend predicts the remaining 19 years better than the country's own line for about three countries in four, and its root mean squared error is less than half as large (0.30 against 0.68 on the log scale). With ten years observed, the advantage is smaller but still clear (0.31 against 0.37).

6.7 The funnel and the non-centred parameterisation

Multilevel models have a computational weak spot that Chapter 5 promised to explain. Consider again the eight countries represented by a single coffee each, but now write the model as Stan code, exactly as the mathematics states it. We draw each country's true average ${\theta }_j$ directly from Normal$(\mu ,\tau )$: the centred parameterisation.

library(rstan)
centred <- "
data { int<lower=1> J; vector[J] ybar; vector<lower=0>[J] se; }
parameters { real mu; real<lower=0> tau; vector[J] theta; }
model {
  mu ~ normal(82, 5); tau ~ exponential(1);
  theta ~ normal(mu, tau);          // centred: country means drawn directly around mu
  ybar ~ normal(theta, se);
}"
noncentred <- "
data { int<lower=1> J; vector[J] ybar; vector<lower=0>[J] se; }
parameters { real mu; real<lower=0> tau; vector[J] z; }
transformed parameters { vector[J] theta = mu + tau * z; }   // non-centred: theta built from z
model {
  mu ~ normal(82, 5); tau ~ exponential(1);
  z ~ normal(0, 1);
  ybar ~ normal(theta, se);
}"
dl <- list(J = nrow(one), ybar = one$mean, se = one$se)
fit_c <- stan(model_code = centred, data = dl, seed = 6, refresh = 0)
fit_n <- stan(model_code = noncentred, data = dl, seed = 6, refresh = 0, control = list(adapt_delta = 0.95))
divergent <- function(f) sum(sapply(get_sampler_params(f, inc_warmup = FALSE), function(x) sum(x[, "divergent__"])))
tibble(model = c("centred", "non-centred"), divergences = c(divergent(fit_c), divergent(fit_n))) |>
  bind_cols(bind_rows(summary(fit_c, pars = "tau")$summary[, c("mean", "n_eff", "Rhat"), drop = FALSE] |> as_tibble(),
                      summary(fit_n, pars = "tau")$summary[, c("mean", "n_eff", "Rhat"), drop = FALSE] |> as_tibble()))

# A tibble: 2 x 5
  model       divergences  mean n_eff  Rhat
  <chr>             <dbl> <dbl> <dbl> <dbl>
1 centred             117 1.07   103.  1.05
2 non-centred           0 0.974 2045.  1.00

The centred model produces over a hundred divergent transitions, an $\widehat{R}$ of 1.05 for $\tau$ and an effective sample size of about a hundred from 4,000 draws. The non-centred model, which is mathematically identical, produces no divergences (with a slightly smaller step size, adapt_delta = 0.95) and an effective sample size twenty times larger. The figure shows why.

pc <- as.data.frame(fit_c) |> transmute(log_tau = log(tau), effect = `theta[1]` - mu) |>
  mutate(divergent = unlist(lapply(get_sampler_params(fit_c, inc_warmup = FALSE), function(x) x[, "divergent__"])) == 1,
         model = "centred: theta[1] - mu")
pn <- as.data.frame(fit_n) |> transmute(log_tau = log(tau), effect = `z[1]`, divergent = FALSE,
                                        model = "non-centred: z[1]")
bind_rows(pc, pn) |> group_by(model, divergent) |> slice_head(n = 1500) |> ungroup() |>   # thin for display
  arrange(divergent) |>
  ggplot(aes(effect, log_tau, colour = divergent)) +
  geom_point(alpha = 0.4, size = 0.8) +
  facet_wrap(~ model, scales = "free_x") +
  scale_colour_manual(values = c(`FALSE` = harbour[["sky"]], `TRUE` = harbour[["rose"]]), labels = c("no", "yes")) +
  labs(x = "country effect", y = "log τ", colour = "divergent") + guides(colour = guide_legend(override.aes = list(size = 2, alpha = 1)))

In the centred parameterisation, the scale of each country effect ${\theta }_j-\mu$ is set by $\tau$. When $\tau$ is large the effects can spread widely; when $\tau$ is near zero they must all be within a hair of $\mu$. Plotted against $\log \tau$, the posterior is a funnel, with a wide mouth at the top and a narrow neck at the bottom. HMC must use one step size for the whole posterior. A step size suited to the wide mouth is far too large for the neck, where the curvature is extreme, so trajectories that try to enter it diverge, and the sampler under-explores small values of $\tau$. This was first described by Radford Neal (2003), and Betancourt and Girolami (2015) analyse its geometry in detail. It is why the centred model's posterior mean for $\tau$ came out too high (1.07, against 0.97 from the non-centred model).

The non-centred parameterisation writes each effect as ${\theta }_j=\mu +\tau z_j$ with $z_j\sim \text{Normal}(0,1)$. The model is identical, but the sampler now explores $z_j$, whose scale does not depend on $\tau$, so the funnel becomes a round hill. It is the same trick as centring a predictor in Chapter 5: change coordinates so that the posterior has no extreme curvature. brms uses the non-centred form for every ( | ) term automatically, which is why the coffee and CO₂ models had no divergences. When you write your own Stan models, or when brms still reports divergences, think about the geometry first.

6.8 A case study: ambulance response times across New York

The techniques of this chapter are made for data like New York City's ambulance dispatch records. The Fire Department publishes every emergency medical incident since 2005 (Fire Department of the City of New York, 2026), more than 30 million of them, on the city's open-data portal, with the call type, a severity code, the borough and ZIP code, and timestamps from the moment the incident is created in the dispatch system to the moment the first unit arrives. We take the most urgent calls of 2024 (severity level 1, which is almost entirely cardiac arrests and choking) and ask a descriptive question: how much does the typical response time vary between the city's ZIP codes, and which ZIP codes are slowest? Some ZIP codes have hundreds of such calls a year and some only a handful, so this is the coffee problem of Section 6.1 with ambulances.

url <- paste0("https://data.cityofnewyork.us/resource/76xm-jjuj.json?",
  "$select=incident_datetime,initial_call_type,borough,zipcode,incident_response_seconds_qy,valid_incident_rspns_time_indc",
  "&$where=initial_severity_level_code='1' AND incident_datetime between '2024-01-01T00:00:00' and '2024-12-31T23:59:59'",
  "&$limit=50000")
ems <- jsonlite::fromJSON(URLencode(url)) |> as_tibble() |>
  filter(valid_incident_rspns_time_indc == "Y", !is.na(zipcode), as.numeric(incident_response_seconds_qy) > 0) |>
  mutate(minutes = as.numeric(incident_response_seconds_qy) / 60, log_min = log(minutes),
         borough = str_to_title(str_remove(borough, "^RICHMOND / ")), date = as.Date(substr(incident_datetime, 1, 10)))
zip_n <- count(ems, zipcode)
c(calls = nrow(ems), zip_codes = nrow(zip_n), min_calls = min(zip_n$n), median_calls = median(zip_n$n), max_calls = max(zip_n$n),
  cardiac_arrest_share = round(mean(ems$initial_call_type == "ARREST"), 3), median_minutes = round(median(ems$minutes), 2))
ems |> group_by(borough) |> summarise(calls = n(), zip_codes = n_distinct(zipcode), median_minutes = round(median(minutes), 2),
                                      over_10_minutes = round(mean(minutes > 10), 3))

               calls            zip_codes            min_calls 
           24313.000              192.000                1.000 
        median_calls            max_calls cardiac_arrest_share 
             102.000              404.000                0.911 
      median_minutes 
               5.570 
# A tibble: 5 x 5
  borough       calls zip_codes median_minutes over_10_minutes
  <chr>         <int>     <int>          <dbl>           <dbl>
1 Bronx          5466        25           5.75           0.077
2 Brooklyn       6522        38           5.45           0.058
3 Manhattan      5568        54           5.1            0.059
4 Queens         5503        65           5.97           0.105
5 Staten Island  1254        12           5.68           0.059

After removing the 2% of calls without a valid response time or a ZIP code, there are 24,313 calls in 192 ZIP codes; 91% are cardiac arrests and most of the rest are choking. The citywide median response time is 5.6 minutes. The ZIP codes differ enormously in how much data they provide: the median ZIP code had 102 such calls, the busiest 404 and the smallest one. The boroughs differ modestly: Manhattan has the fastest median (5.1 minutes) and Queens the slowest (6.0 minutes), and in Queens 10.5% of responses took longer than ten minutes, against 6% to 8% elsewhere.

Response times are skewed, with a long right tail of delayed responses, so we model the log of the time in minutes. ZIP codes are nested in boroughs, and both get varying intercepts. Because a few calls in any ZIP code can take far longer than the rest, we use a Student-$t$ likelihood on the log scale (Chapter 3), which lets those calls be unusual without dragging the ZIP estimates:

The prior on $\alpha$ is Normal(1.7, 0.5), centred on about 5.5 minutes, and the group standard deviations have Exponential(5) priors, which put most of their mass on differences between places of less than about 30%.

ems_priors <- c(prior(normal(1.7, 0.5), class = Intercept), prior(exponential(5), class = sd),
                prior(exponential(2), class = sigma), prior(gamma(2, 0.1), class = nu))
m_ems <- brm(log_min ~ 1 + (1 | borough) + (1 | zipcode), data = ems, family = student(), prior = ems_priors,
             seed = 6, refresh = 0, file = "fits/ch6_ems")
summary(m_ems)$fixed[, c("Estimate", "l-95% CI", "u-95% CI", "Rhat")] |> round(3)
VarCorr(m_ems) |> map(~ round(.x$sd[, c("Estimate", "Q2.5", "Q97.5")], 3))
summary(m_ems)$spec_pars[, c("Estimate", "l-95% CI", "u-95% CI")] |> round(2)

          Estimate l-95% CI u-95% CI  Rhat
Intercept    1.721    1.618    1.822 1.004
$borough
Estimate     Q2.5    Q97.5 
   0.097    0.041    0.220 

$zipcode
Estimate     Q2.5    Q97.5 
   0.129    0.113    0.148 

$residual__
Estimate     Q2.5    Q97.5 
   0.355    0.350    0.360 

      Estimate l-95% CI u-95% CI
sigma     0.35     0.35     0.36
nu        4.49     4.25     4.74

The typical response to a life-threatening call is $e^{1.72}\approx 5.6$ minutes (95% interval 5.0 to 6.2). ZIP codes within a borough vary with a standard deviation of 0.13 on the log scale (0.11 to 0.15): a ZIP code one standard deviation slower than its borough's typical one takes about 14% longer, roughly 45 seconds. The borough standard deviation, 0.10, is estimated from only five boroughs and its interval is wide (0.04 to 0.22), the situation Section 6.9's practice box warns about. The Student-$t$ degrees of freedom, $\nu \approx 4.5$, confirm that the tails are heavy: a normal model on the log scale would have been pulled around by a minority of very long responses. All $\widehat{R}$ values are at most 1.004.

zip_borough <- distinct(ems, zipcode, borough)
newd <- zip_borough
post_zip <- posterior_linpred(m_ems, newdata = newd)                          # draws x ZIP codes, log minutes
zips <- newd |> mutate(post = exp(apply(post_zip, 2, median)), lo = exp(apply(post_zip, 2, quantile, 0.05)),
                       hi = exp(apply(post_zip, 2, quantile, 0.95))) |>
  left_join(ems |> group_by(zipcode) |> summarise(n = n(), raw = median(minutes)), by = "zipcode")
city <- exp(fixef(m_ems)["Intercept", "Estimate"])
ggplot(zips, aes(n)) + geom_hline(yintercept = city, linetype = 2, colour = harbour[["muted"]]) +
  geom_segment(aes(xend = n, y = raw, yend = post), colour = harbour[["grid"]]) +
  geom_point(aes(y = raw), shape = 1, colour = harbour[["muted"]], size = 1.5) +
  geom_point(aes(y = post, colour = borough), size = 1.6) + scale_x_log10() +
  scale_colour_manual(values = unname(harbour[c("sky", "rose", "amber", "moss", "lilac")])) +
  labs(x = "life-threatening calls in 2024 (log scale)", y = "typical response time (minutes)", colour = NULL)
bind_rows(zips |> arrange(desc(raw)) |> head(5) |> mutate(ranking = "raw median"),
          zips |> arrange(desc(post)) |> head(5) |> mutate(ranking = "posterior median")) |>
  select(ranking, zipcode, borough, n, raw, post, lo, hi) |> mutate(across(c(raw, post, lo, hi), ~ round(.x, 2))) |> as.data.frame()

            ranking zipcode   borough  n   raw post   lo    hi
1        raw median   10112 Manhattan  1 19.07 5.35 4.33  6.68
2        raw median   11697    Queens 19 11.35 9.52 8.43 10.67
3        raw median   10020 Manhattan  6 10.57 6.43 5.29  7.68
4        raw median   10178 Manhattan  1 10.22 5.39 4.36  6.71
5        raw median   10464     Bronx 15  8.90 7.51 6.57  8.55
6  posterior median   11697    Queens 19 11.35 9.52 8.43 10.67
7  posterior median   11360    Queens 44  8.03 7.67 7.01  8.40
8  posterior median   11414    Queens 73  8.13 7.62 7.03  8.22
9  posterior median   11356    Queens 41  7.73 7.51 6.76  8.32
10 posterior median   10464     Bronx 15  8.90 7.51 6.57  8.55

The figure is Section 6.3's shrinkage plot again, with ambulances. Ranked by their raw medians, three of the five slowest ZIP codes have six calls or fewer, and two have just one: ZIP code 10112 in Midtown Manhattan, which had a single response of 19 minutes, is estimated at 5.4 minutes once it is pooled with the rest of Manhattan. Ranked by posterior medians, the slowest ZIP codes are different and make geographic sense: Breezy Point at the western tip of the Rockaway peninsula (11697: 9.5 minutes, 90% interval 8.4 to 10.7, from 19 calls), followed by Bayside (11360), Howard Beach (11414) and College Point (11356) in Queens and City Island (10464) in the Bronx, all at the edges of the city with limited road access. These estimates rest on enough calls to survive pooling, and they are where an analysis of response-time equity should start looking.

As in Section 6.4, the case for pooling is predictive, and it can be tested: fit the model to the first half of 2024 and predict each ZIP code's typical response in the second half.

h1 <- filter(ems, date < as.Date("2024-07-01")); h2 <- filter(ems, date >= as.Date("2024-07-01"))
m_h1 <- update(m_ems, newdata = h1, seed = 6, refresh = 0, file = "fits/ch6_ems_h1")
truth <- h2 |> group_by(zipcode) |> summarise(n2 = n(), target = median(log_min))
raw1 <- h1 |> group_by(zipcode) |> summarise(n1 = n(), raw = median(log_min))
pooled <- zip_borough |> filter(zipcode %in% raw1$zipcode) |>
  mutate(pooled = colMeans(posterior_linpred(m_h1, newdata = zip_borough |> filter(zipcode %in% raw1$zipcode))))
cmp <- truth |> inner_join(raw1, by = "zipcode") |> inner_join(pooled, by = "zipcode")
err <- function(d) c(zip_codes = nrow(d), rmse_raw = sqrt(weighted.mean((d$raw - d$target)^2, d$n2)),
                     rmse_pooled = sqrt(weighted.mean((d$pooled - d$target)^2, d$n2)))
rbind(all = err(cmp), fewer_than_30_calls = err(filter(cmp, n1 < 30)), at_least_30_calls = err(filter(cmp, n1 >= 30))) |> round(3)

                    zip_codes rmse_raw rmse_pooled
all                       187    0.097       0.096
fewer_than_30_calls        43    0.184       0.176
at_least_30_calls         144    0.088       0.088

For the 187 ZIP codes with calls in both halves of the year, the pooled estimates predict the second half's typical response slightly better than the raw first-half medians (weighted RMSE on the log scale 0.096 against 0.097). The gain is concentrated where it should be: for the 43 ZIP codes with fewer than 30 calls in the first half, the error falls from 0.184 to 0.176, while for the 144 larger ones the two are identical (0.088). Most ZIP codes have enough calls for their own median to be reliable, so pooling changes little for them; its value here is in protecting rankings, maps and league tables from ZIP codes with a handful of calls, as the raw top five show.

Everything here is descriptive. The model says where responses are slow, not why. Whether neighbourhood composition causes slower responses, the equity question behind this data, is a causal question with the structure of Section 1.9: distance to the nearest ambulance station, traffic, call volume and hospital locations all lie on paths between a neighbourhood and its response times, and some of them are mediators that an equity analysis should not adjust away. A multilevel model like this one is the right outcome model for that analysis; the graph decides what goes into it. Note also what the clock measures: it starts when the incident is created in the dispatch system, after the call has been answered and classified, so any delay before that point is invisible in these data.

6.9 Multilevel models in practice

McElreath (2020) argues that multilevel regression deserves to be the default form of regression, and that a paper should have to justify not using it. Whenever data have clusters (repeated measures on people, students in classes, events in places, measurements over time), single-level models either ignore the clustering, which makes standard errors far too small, or treat each cluster as unrelated, which wastes information and overfits. The models of this chapter extend naturally:

• Several grouping factors. Students belong to classes and are also assessed by markers; a model with (1 | class) + (1 | marker) has crossed varying effects. Students within classes within schools give nested effects, (1 | school/class).
• Any likelihood. Varying effects work identically in logistic and Poisson models: recid ~ age + (1 | county) or calls ~ offset(log(pop)) + (1 | neighbourhood).
• Multilevel regression and poststratification (MRP). To estimate a population quantity from a non-representative sample, such as an opinion poll or the NHANES data of Chapter 3, fit a multilevel model for the outcome by demographic cell and reweight the cell estimates by census counts. MRP is now standard in public opinion research (Gelman and Hill, 2007).

6.10 Exercises

Each exercise has a hint and an answer button.

Further reading and references

1. Betancourt, M. and Girolami, M. (2015). Hamiltonian Monte Carlo for hierarchical models. In S. K. Upadhyay et al. (eds), Current Trends in Bayesian Methodology with Applications, 79–101. CRC Press.
2. Efron, B. and Morris, C. (1975). Data analysis using Stein's estimator and its generalizations. Journal of the American Statistical Association, 70(350), 311–319.
3. Efron, B. and Morris, C. (1977). Stein's paradox in statistics. Scientific American, 236(5), 119–127.
4. Fire Department of the City of New York (2026). EMS Incident Dispatch Data. NYC Open Data, dataset 76xm-jjuj, data.cityofnewyork.us.
5. Friedlingstein, P. et al. (2023). Global carbon budget 2023. Earth System Science Data, 15(12), 5301–5369. Data via Our World in Data: github.com/owid/co2-data.
6. Gelman, A. and Hill, J. (2007). Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.
7. Goldstein, H. and Spiegelhalter, D. J. (1996). League tables and their limitations: statistical issues in comparisons of institutional performance. Journal of the Royal Statistical Society A, 159(3), 385–443.
8. James, W. and Stein, C. (1961). Estimation with quadratic loss. Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability, 1, 361–379.
9. LeDoux, J. (2018). Coffee Quality Institute database. Distributed via TidyTuesday, 7 July 2020: github.com/rfordatascience/tidytuesday.
10. Lewandowski, D., Kurowicka, D. and Joe, H. (2009). Generating random correlation matrices based on vines and extended onion method. Journal of Multivariate Analysis, 100(9), 1989–2001.
11. McElreath, R. (2020). Statistical Rethinking (2nd ed.). CRC Press.
12. Neal, R. M. (2003). Slice sampling. Annals of Statistics, 31(3), 705–767.
13. Robbins, H. (1956). An empirical Bayes approach to statistics. Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, 1, 157–163.
14. Stein, C. (1956). Inadmissibility of the usual estimator for the mean of a multivariate normal distribution. Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, 1, 197–206.
15. Wainer, H. (2007). The most dangerous equation. American Scientist, 95(3), 249–256.