Chapter 4
Generalized linear models

Lectures 7–8 · probabilities and counts: links, logistic and Poisson regression, and the fairness of a risk score

The linear model of Chapter 3 says that an outcome scatters in a bell shape around a straight line. That is a fine description of heights and log earnings. It is a poor description of whether a defendant is re-arrested (a yes or a no), how many incidents an airline has in fifteen years (a count, never negative, often zero), or which of several options a voter chooses. For outcomes like these the normal distribution is the wrong likelihood, and a straight line for the mean will sooner or later predict impossible values: probabilities above one, negative counts.

The generalized linear model (GLM) keeps everything that made the linear model useful (a linear predictor built from the variables you choose, priors on its coefficients, a causal graph that says which variables belong) and changes two things: the likelihood, which is chosen to respect the kind of outcome you have, and a link function, which bends the linear predictor so that it always lands in the range the outcome allows. Logistic regression for yes-or-no outcomes and Poisson regression for counts are the two most important members of the family, and between them they cover a very large share of applied research.

The price of the link function is that coefficients no longer mean what they meant in Chapter 3. A coefficient of 0.8 is no longer "0.8 more of the outcome"; it is 0.8 on the scale of log-odds or log-rates, and what that means for the outcome depends on everything else. Learning to translate coefficients back to the outcome scale, and to report results there, is most of the skill of using GLMs well.

• Lecture 7 is about probabilities. We will build logistic regression from the link function up, set priors on the log-odds scale, and fit it to the data at the centre of the most famous recent controversy in applied statistics: the COMPAS criminal risk score, analysed by ProPublica in 2016. We will learn to report effects as differences in probability, and to check whether a model's probabilities can be trusted.
• Lecture 8 is about fairness and counts. We will see why the COMPAS debate could not be settled by better statistics, and what that teaches about choosing an estimand. Then we turn to counts, with fifteen years of incidents at 56 airlines, and meet exposure, overdispersion, the negative binomial and regression to the mean.

Lecture 7 — Probabilities

4.1 When a straight line breaks

Suppose $y_i$ is 1 if defendant $i$ is re-arrested within two years and 0 otherwise. The quantity we want to model is the probability $p_i=Pr(y_i=1)$. The simplest idea is to make $p_i$ a straight line in a predictor such as the number of prior arrests. This linear probability model is common in economics and often harmless near the middle of the data, but it has a fatal flaw: lines go on forever, and probabilities do not. A line steep enough to describe the difference between zero and five prior arrests will predict a probability above 1 for someone with thirty.

The GLM solution has three ingredients, and you have already met the reasoning behind each.

1. A likelihood that matches the outcome. A yes-or-no outcome with probability $p_i$ is a Bernoulli variable; a count of yeses out of $n$ independent tries is binomial. Section 2.8 showed these are the maximum-entropy distributions for binary outcomes with a given mean, so they assume nothing beyond the probability itself.
2. A linear predictor, exactly as in Chapter 3: ${\eta }_i=\alpha +{\beta }_1{x}_{1i}+\cdots +{\beta }_K{x}_{Ki}$, which can be any real number.
3. A link function connecting the two. For probabilities, the standard choice is the logit, the logarithm of the odds: $$logit(p_i)=\log \frac{p_i}{1-p_i}={\eta }_i⟺p_i=\frac{1}{1+e^{-{\eta }_i}}.$$ The right-hand function, the inverse logit or logistic function, squeezes the whole real line into the interval from 0 to 1.

The logistic function is S-shaped: nearly flat when $\eta$ is very negative (the probability is close to 0 and hard to move), steepest at $\eta =0$ (where $p=1/2$), and nearly flat again when $\eta$ is large. So a coefficient no longer has a single meaning on the probability scale. A change of one unit in $x$ always changes the log-odds by $\beta$, but it changes the probability a lot when $p$ is near one-half and very little when $p$ is near 0 or 1. Explore this below.

Two facts from the widget are worth memorising. First, the steepest slope of the probability curve, reached where $p=1/2$, is exactly $\beta /4$. This gives the divide-by-4 rule (Gelman, Hill and Vehtari, 2020): a logistic coefficient divided by 4 is an upper bound on the change in probability per unit of $x$. A coefficient of 0.8 means "at most about 20 percentage points per unit". Second, far from $p=1/2$ the effect on probability is much smaller, so the same coefficient implies different changes in probability for different people. Logistic regression has interactions built in: the effect of one predictor on the probability depends on the values of all the others.

4.2 A risk score and its data

Courts across the United States use risk assessment instruments to inform decisions about bail, sentencing and parole. One of the most widely used is COMPAS, sold by the company Northpointe (now Equivant). A defendant answers a questionnaire of 137 items, and the software produces a "decile score" from 1 (lowest risk) to 10 (highest) for the risk of committing another crime. In 2016 the news organisation ProPublica obtained the scores of more than 7,000 people arrested in Broward County, Florida, in 2013 and 2014, linked them to criminal records, and recorded whether each person was arrested again within two years. Their article, "Machine Bias" (Angwin et al., 2016), argued that the scores were biased against Black defendants. The data and code were published on GitHub, and they have since become a standard test case for fairness in prediction.

We load the data with exactly the filters ProPublica used: keep people whose score was computed within 30 days of their arrest, whose two-year outcome is known, whose charge was an ordinary felony or misdemeanour, and whose score is valid.

library(tidyverse)
library(brms)
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/propublica/compas-analysis/master/compas-scores-two-years.csv"
compas <- read_csv(url, show_col_types = FALSE, name_repair = make.unique) |>
  filter(between(days_b_screening_arrest, -30, 30),   # scored near the time of arrest
         is_recid != -1,                               # outcome known
         c_charge_degree != "O",                       # ordinary charges
         score_text != "N/A") |>                       # valid score
  transmute(id, sex, age, race, priors = priors_count, decile = decile_score,
            recid = two_year_recid)                    # 1 = arrested again within 2 years
nrow(compas)
compas |> group_by(race) |>
  summarise(n = n(), recid = mean(recid), age = mean(age), priors = mean(priors)) |>
  arrange(-n)

[1] 6172
# A tibble: 6 x 5
  race                 n recid   age priors
  <chr>            <int> <dbl> <dbl>  <dbl>
1 African-American  3175 0.523  32.4   4.24
2 Caucasian         2103 0.391  37.5   2.29
3 Hispanic           509 0.371  35.0   2.10
4 Other              343 0.362  34.8   1.75
5 Asian               31 0.258  38.2   1.35
6 Native American     11 0.455  33.2   5.18

We have 6,172 people. Over two years, 45% were arrested again. The two largest groups differ substantially: 52% of the 3,175 Black defendants and 39% of the 2,103 white defendants were re-arrested. They also differ in age (the Black defendants are on average five years younger) and in recorded prior arrests (4.2 against 2.3 on average). Those differences will matter a great deal.

4.3 Priors on the log-odds scale

Priors for a logistic regression live on the log-odds scale, where intuition is weak. Section 2.9 warned that a "vague" prior on one scale can be strongly informative on another. Here is the classic example: the prior Normal(0, 10) on a log-odds intercept, which looks agnostic, and the prior Normal(0, 1.5), which looks restrictive. What do they say about the probability?

set.seed(1)
prior_p <- tibble(sd = rep(c(10, 1.5), each = 1e4)) |>
  mutate(alpha = rnorm(n(), 0, sd), p = plogis(alpha),
         prior = factor(paste0("alpha ~ Normal(0, ", sd, ")"),
                        levels = c("alpha ~ Normal(0, 10)", "alpha ~ Normal(0, 1.5)")))
prior_p |> group_by(prior) |> summarise(`share below 0.01 or above 0.99` = mean(p < 0.01 | p > 0.99))
ggplot(prior_p, aes(p)) +
  geom_histogram(binwidth = 0.025, boundary = 0, fill = harbour[["sky"]], alpha = 0.7) +
  facet_wrap(~ prior) + labs(x = "implied prior probability p", y = "draws")

# A tibble: 2 x 2
  prior                  `share below 0.01 or above 0.99`
  <fct>                                             <dbl>
1 alpha ~ Normal(0, 10)                            0.644 
2 alpha ~ Normal(0, 1.5)                           0.0023

The "vague" prior puts about two-thirds of its weight on probabilities below 1% or above 99%: a very strong and silly opinion about re-arrest. The Normal(0, 1.5) prior is close to flat on the probability scale, which is what most people mean by "I have no idea". The lesson generalises: in a logistic regression, priors with standard deviations of about 1 to 1.5 on the intercept and on coefficients of standardised predictors are weakly informative, and anything much wider pushes predictions toward 0 and 1.

4.4 Fitting a logistic regression

Our first model predicts re-arrest from three things a court knows about every defendant: age, number of prior arrests and sex. Before fitting, think about the scales. Age enters in decades from 35, so the coefficient is the change in log-odds per ten years. Prior arrests are very skewed (most people have none or a few; some have thirty), and the difference between 0 and 3 priors plausibly matters more than the difference between 20 and 23, so we use $\log (1+\text{priors})$. The model is

compas <- compas |>
  mutate(age10 = (age - 35) / 10, lpriors = log1p(priors), male = as.integer(sex == "Male"))
m_lr <- brm(recid ~ 1 + age10 + lpriors + male, data = compas, family = bernoulli(),
            prior = c(prior(normal(0, 1.5), class = Intercept),
                      prior(normal(0, 1), class = b)),
            seed = 4, refresh = 0, file = "fits/ch4_logit")
fixef(m_lr) |> round(2)
exp(fixef(m_lr)[-1, -2]) |> round(2)        # odds ratios with 95% intervals

          Estimate Est.Error  Q2.5 Q97.5
Intercept    -1.34      0.07 -1.48 -1.20
age10        -0.46      0.03 -0.52 -0.41
lpriors       0.85      0.03  0.78  0.91
male          0.32      0.07  0.18  0.46
        Estimate Q2.5 Q97.5
age10       0.63 0.60  0.66
lpriors     2.33 2.18  2.49
male        1.38 1.20  1.59

Every coefficient is far from zero, with narrow intervals. Read them on three scales.

• Log-odds (the raw coefficients): ten years of age lower the log-odds of re-arrest by about 0.46; each unit of $\log (1+\text{priors})$ raises it by about 0.85; being male raises it by about 0.32. Direction and relative size are visible here, but not much else.
• Odds ratios ($e^{\beta }$): ten more years of age multiply the odds of re-arrest by about 0.63; being male multiplies them by about 1.38. Odds ratios are standard in medicine and epidemiology, and Section 4.5 explains why they often mislead.
• Probabilities: the scale people care about. To get there, compute predicted probabilities from every posterior draw and summarise them.

grid_p <- expand_grid(priors = 0:20, age = c(20, 35, 50)) |>
  mutate(age10 = (age - 35) / 10, lpriors = log1p(priors), male = 1)
pred <- posterior_epred(m_lr, newdata = grid_p)            # draws x rows of grid_p
grid_p <- grid_p |>
  mutate(p = colMeans(pred), lo = apply(pred, 2, quantile, 0.025),
         hi = apply(pred, 2, quantile, 0.975), age = factor(age))
ggplot(grid_p, aes(priors, p, colour = age, fill = age)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), colour = NA, alpha = 0.25) +
  geom_line(linewidth = 1) +
  scale_colour_manual(values = c(`20` = harbour[["rose"]], `35` = harbour[["amber"]], `50` = harbour[["sky"]])) +
  scale_fill_manual(values = c(`20` = harbour[["rose"]], `35` = harbour[["amber"]], `50` = harbour[["sky"]])) +
  scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
  labs(x = "prior arrests", y = "Pr(re-arrest within 2 years)", colour = "age", fill = "age")

A man of 20 with no prior arrests has about a 42% chance of re-arrest within two years; a man of 50 with none, about 15%. The first few prior arrests raise the probability steeply; after ten or so, extra priors change it little, because of the logarithm and because probabilities near the top of the S-curve are hard to move. These are the statements a reader can act on, and none of them can be read off a coefficient table.

4.5 Reporting on the scale people care about

"How much does age matter?" now has many answers, one for each defendant. The standard summary is the average comparison (also called the average marginal effect): for every defendant in the data, predict the probability of re-arrest as they are and as they would be ten years older, take the difference and average over defendants. Doing this for every posterior draw gives a posterior distribution for the average difference.

set.seed(4)
p_now   <- posterior_epred(m_lr, ndraws = 1000)                                     # as observed
p_older <- posterior_epred(m_lr, ndraws = 1000, newdata = mutate(compas, age10 = age10 + 1))
avg_diff  <- rowMeans(p_older - p_now)                     # average change in probability
avg_ratio <- rowMeans(p_older) / rowMeans(p_now)           # ratio of average probabilities
rbind(`risk difference` = quantile(avg_diff, c(0.5, 0.025, 0.975)),
      `risk ratio`      = quantile(avg_ratio, c(0.5, 0.025, 0.975))) |> round(3)

                   50%   2.5%  97.5%
risk difference -0.094 -0.112 -0.075
risk ratio       0.793  0.757  0.833

Being ten years older lowers the probability of re-arrest by about 9 percentage points on average (95% interval 7.5 to 11.2), or multiplies it by about 0.79. Compare that with the odds ratio of 0.63. The odds ratio suggests a 37% reduction; the risk actually falls by about 21%. This gap is not an error. When an outcome is common (here 45% of people are re-arrested), odds and probabilities are very different things, and odds ratios exaggerate the corresponding ratios of probabilities. Only when the outcome is rare (below about 10%) are the two close.

4.6 Can the probabilities be trusted?

A model that outputs probabilities makes a promise: among people given a probability of 30%, about 30% should experience the outcome. A model that keeps this promise is calibrated. Calibration is the first thing to check in any probability model, and it is easy to check: group people by their predicted probability and compare the predicted and observed rates in each group.

We can also ask how well the model discriminates: whether people who were re-arrested got higher probabilities than people who were not. The standard summary is the area under the ROC curve (AUC): the probability that a randomly chosen person who was re-arrested has a higher predicted probability than a randomly chosen person who was not. An AUC of 0.5 is a coin toss; 1 is perfect.

This makes possible a revealing comparison. COMPAS uses 137 items and a proprietary algorithm. Our model uses three facts available in any court record. Which predicts better? We fit a second logistic regression using only the COMPAS decile score, so that both are on the same probability scale.

m_dec <- update(m_lr, formula = recid ~ 1 + decile, newdata = compas,
                seed = 4, refresh = 0, file = "fits/ch4_decile")
auc <- function(p, y) {                       # probability that a random positive outranks a random negative
  r <- rank(p); n1 <- sum(y); n0 <- sum(1 - y)
  (sum(r[y == 1]) - n1 * (n1 + 1) / 2) / (n1 * n0)
}
compas <- compas |>
  mutate(p_model = colMeans(posterior_epred(m_lr)), p_compas = colMeans(posterior_epred(m_dec)))
compas |> summarise(auc_model = auc(p_model, recid), auc_compas = auc(p_compas, recid),
                    brier_model = mean((p_model - recid)^2), brier_compas = mean((p_compas - recid)^2))

# A tibble: 1 x 4
  auc_model auc_compas brier_model brier_compas
      <dbl>      <dbl>       <dbl>        <dbl>
1     0.732      0.710       0.208        0.215

calib <- compas |>
  pivot_longer(c(p_model, p_compas), names_to = "model", values_to = "p") |>
  mutate(model = recode(model, p_model = "age + priors + sex", p_compas = "COMPAS decile only")) |>
  group_by(model) |> mutate(bin = ntile(p, 10)) |>
  group_by(model, bin) |>
  summarise(pred = mean(p), k = sum(recid), n = n(), .groups = "drop") |>
  mutate(obs = k / n, lo = qbeta(0.025, k + 1, n - k + 1), hi = qbeta(0.975, k + 1, n - k + 1))
ggplot(calib, aes(pred, obs, colour = model)) +
  geom_abline(linetype = "dashed", colour = harbour[["muted"]]) +
  geom_pointrange(aes(ymin = lo, ymax = hi), position = position_dodge(width = 0.015)) +
  scale_colour_manual(values = c(harbour[["sky"]], harbour[["amber"]])) +
  coord_equal(xlim = c(0, 1), ylim = c(0, 1)) +
  labs(x = "predicted probability", y = "observed re-arrest rate", colour = NULL)

Both models are well calibrated: within each bin, the observed rate is close to the predicted probability. And the three-variable model discriminates slightly better than the COMPAS score (AUC 0.73 against 0.71), with a slightly lower Brier score (0.208 against 0.215; the Brier score is the mean squared difference between probability and outcome). An AUC of about 0.7 means that if you pick one person who was re-arrested and one who was not, the model ranks them correctly about 70% of the time. That is useful but far from certain, and both approaches are doing it with essentially the same information: how young a person is and how many times they have been arrested before.

Lecture 8 — Fairness, counts and exposure

4.7 What would a fair risk score be?

ProPublica's central finding is easy to state. Call a defendant "flagged" if their COMPAS decile score is 5 or more (medium or high risk). Among defendants who were not re-arrested, a much larger share of Black than white defendants had been flagged: the tool's false positive rate was higher for Black defendants. Among those who were re-arrested, a larger share of white defendants had been missed: the false negative rate was higher for white defendants. Northpointe replied (Dieterich, Mendoza and Brennan, 2016) that the score was fair in a different sense: at every score, Black and white defendants were re-arrested at about the same rate, so a score meant the same thing whoever received it. That is calibration within groups, the property of Section 4.6, checked separately for each group. Both claims can be checked in the data, with Bayesian uncertainty.

two <- compas |> filter(race %in% c("African-American", "Caucasian")) |>
  mutate(race = recode(race, `African-American` = "Black", Caucasian = "white"),
         flagged = decile >= 5)
counts <- two |> group_by(race) |>
  summarise(fp = sum(flagged & recid == 0), neg = sum(recid == 0),     # false positives / all negatives
            fn = sum(!flagged & recid == 1), pos = sum(recid == 1),    # false negatives / all positives
            tp = sum(flagged & recid == 1), flag = sum(flagged))       # true positives / all flagged
set.seed(7)
draws <- counts |> rowwise() |>                                        # Beta(1, 1) priors on each rate
  reframe(race, FPR = rbeta(4000, fp + 1, neg - fp + 1), FNR = rbeta(4000, fn + 1, pos - fn + 1),
          PPV = rbeta(4000, tp + 1, flag - tp + 1))
draws |> pivot_longer(-race, names_to = "rate") |> group_by(rate, race) |>
  summarise(mean = mean(value), lo = quantile(value, 0.025), hi = quantile(value, 0.975), .groups = "drop")
draws |> group_by(race) |> mutate(i = row_number()) |> ungroup() |>
  pivot_wider(names_from = race, values_from = c(FPR, FNR, PPV)) |>
  summarise(`FPR: Black - white` = mean(FPR_Black - FPR_white),
            `FNR: Black - white` = mean(FNR_Black - FNR_white),
            `PPV: Black - white` = mean(PPV_Black - PPV_white))

# A tibble: 6 x 5
  rate  race   mean    lo    hi
  <chr> <chr> <dbl> <dbl> <dbl>
1 FNR   Black 0.285 0.263 0.307
2 FNR   white 0.496 0.462 0.531
3 FPR   Black 0.424 0.399 0.449
4 FPR   white 0.221 0.198 0.244
5 PPV   Black 0.649 0.627 0.671
6 PPV   white 0.594 0.558 0.631
# A tibble: 1 x 3
  `FPR: Black - white` `FNR: Black - white` `PPV: Black - white`
                 <dbl>                <dbl>                <dbl>
1                0.203               -0.211               0.0548

Both sides were right. Among people not re-arrested, 42% of Black defendants and 22% of white defendants had been flagged. Among people re-arrested, 28% of Black and 50% of white defendants had been missed. But among people flagged, 65% of Black and 59% of white defendants were re-arrested (the positive predictive value), and the calibration curves by group nearly coincide. The score means roughly the same thing for both groups, and it makes very different kinds of mistakes for each.

cal_g <- two |> group_by(race, decile) |>
  summarise(k = sum(recid), n = n(), .groups = "drop") |>
  mutate(obs = k / n, lo = qbeta(0.025, k + 1, n - k + 1), hi = qbeta(0.975, k + 1, n - k + 1))
p1 <- ggplot(cal_g, aes(decile, obs, colour = race)) +
  geom_line() + geom_pointrange(aes(ymin = lo, ymax = hi), size = 0.25,
                                position = position_dodge(width = 0.3)) +
  scale_x_continuous(breaks = 1:10) + scale_y_continuous(labels = scales::percent) +
  labs(x = "COMPAS decile score", y = "re-arrested", colour = NULL)
p2 <- draws |> pivot_longer(c(FPR, FNR), names_to = "rate") |>
  ggplot(aes(value, fill = race)) + geom_density(colour = NA, alpha = 0.7) +
  facet_wrap(~ rate, ncol = 1, scales = "free_y") + scale_x_continuous(labels = scales::percent) +
  labs(x = "error rate at cut-off 5", y = NULL, fill = NULL) + theme(axis.text.y = element_blank()) + guides(fill = "none")
patchwork::wrap_plots(p1, p2, widths = c(1.2, 1)) + patchwork::plot_layout(guides = "collect") &
  scale_colour_manual(values = c(Black = harbour[["rose"]], white = harbour[["sky"]])) &
  scale_fill_manual(values = c(Black = harbour[["rose"]], white = harbour[["sky"]])) &
  theme(legend.position = "top")

Why can't a score have both properties? The reason is arithmetic, and it was proved independently in 2016–17 by Jon Kleinberg, Sendhil Mullainathan and Manish Raghavan, and by Alexandra Chouldechova (Kleinberg, Mullainathan and Raghavan, 2017; Chouldechova, 2017). If two groups have different base rates of the outcome (here 52% against 39% re-arrested), then no imperfect score can simultaneously be calibrated within both groups and have equal false positive and false negative rates in both. A calibrated score has to give higher scores, on average, to the group with the higher base rate, and with any fixed cut-off that group will then contain more flagged people who turn out not to be re-arrested. The only escapes are a perfect predictor or equal base rates. The widget lets you try to beat the theorem.

With a single cut-off, PPV is similar for the two groups but the error rates differ. Press "equalise false positive rates": the cut-off for Black defendants rises until the FPRs match, but now the PPVs and FNRs differ, and a higher cut-off is being applied to one group, which is itself a form of unequal treatment. There is no setting that equalises everything. Which kind of fairness should a court prefer? That is not a statistical question; it is a question about which errors are worse, and for whom. What statistics can do is make the trade-off explicit, with uncertainty.

What about the race coefficient?

A natural next step is to add race to the logistic regression and look at its coefficient. Here is the result, and then the reasons to be very careful with it.

two <- two |> mutate(black = as.integer(race == "Black"))
m_race <- brm(recid ~ 1 + black, data = two, family = bernoulli(),
              prior = c(prior(normal(0, 1.5), class = Intercept), prior(normal(0, 1), class = b)),
              seed = 4, refresh = 0, file = "fits/ch4_race")
m_race_adj <- update(m_race, formula = recid ~ 1 + black + age10 + lpriors + male, newdata = two,
                     seed = 4, refresh = 0, file = "fits/ch4_race_adj")
rbind(`race only` = fixef(m_race)["black", ], `adjusting for age, priors, sex` = fixef(m_race_adj)["black", ]) |>
  round(2)

                               Estimate Est.Error  Q2.5 Q97.5
race only                          0.54      0.06  0.42  0.65
adjusting for age, priors, sex     0.07      0.06 -0.06  0.18

Alone, being Black is associated with log-odds of re-arrest higher by about 0.54 (52% against 39%). After adjusting for age, prior arrests and sex, the coefficient is about 0.07 (95% interval −0.06 to 0.18), which includes zero. Among people of the same age, sex and number of recorded prior arrests, Black and white defendants in these data were re-arrested at about the same rate.

What does that mean? Draw the graph, using what is known about how these data are made.

The coefficient's meaning depends on which arrows are real. If recorded prior arrests reflect both behaviour and the intensity of policing, then "adjusting for prior arrests" partly adjusts away a pathway through policing, and a near-zero coefficient cannot be read as "no racial disparity in the system". It says only that, given the recorded history, recorded re-arrest rates are similar: exactly the calibration property again. And because the outcome is itself an arrest, the path $R\to P\to Y$ affects what counts as "re-arrest" in the first place. Whether the coefficient is a direct effect (Section 3.7), a descriptive conditional comparison, or nothing interpretable at all depends on an estimand you must state and on arrows you must defend. Race is also not a treatment that can be assigned, which is why careful studies frame questions around manipulable things: a policing practice, a policy, the use of a tool.

4.8 Counts and exposure: airline safety

The second great family of GLMs is for counts: numbers of events such as accidents, arrests, emergency calls, citations or infections. Counts are whole numbers, never negative, often zero, and their spread grows with their mean. The natural likelihood is the Poisson distribution, the limit of the binomial when there are many opportunities for an event, each with a small probability, so that it inherits the binomial's maximum-entropy justification (Section 2.8). Its single parameter $\lambda$ is both its mean and its variance. Because $\lambda$ must be positive, the link is the logarithm: $\log {\lambda }_i={\eta }_i$.

Our data are from a 2014 FiveThirtyEight analysis by Nate Silver: for 56 major airlines, the number of safety incidents (events in which an aircraft is damaged or someone is injured) in 1985–1999 and in 2000–2014, together with each airline's capacity in available seat-kilometres per week. The question is one every traveller has asked: is an airline's safety record a lasting trait? If an airline had many incidents in the past, should you expect more from it in the future?

url2 <- "https://raw.githubusercontent.com/fivethirtyeight/data/master/airline-safety/airline-safety.csv"
air <- read_csv(url2, show_col_types = FALSE) |>
  mutate(exposure = avail_seat_km_per_week * 52 * 15 / 1e9,     # billions of seat-km per 15 years
         past = log((incidents_85_99 + 0.5) / exposure),         # log past incident rate
         past_z = (past - mean(past)) / sd(past))
nrow(air)
air |> summarise(across(c(incidents_85_99, incidents_00_14), list(mean = mean, var = var)))
air |> arrange(-incidents_85_99 / exposure) |>
  select(airline, exposure, incidents_85_99, incidents_00_14) |> head(5)

[1] 56
# A tibble: 1 x 4
  incidents_85_99_mean incidents_85_99_var incidents_00_14_mean
                 <dbl>               <dbl>                <dbl>
1                 7.18                122.                 4.12
# i 1 more variable: incidents_00_14_var <dbl>
# A tibble: 5 x 4
  airline                exposure incidents_85_99 incidents_00_14
  <chr>                     <dbl>           <dbl>           <dbl>
1 Aeroflot*                  934.              76               6
2 Ethiopian Airlines         381.              25               5
3 Pakistan International     272.               8              10
4 Xiamen Airlines            336.               9               2
5 Philippine Airlines        322.               7               2

Airlines differ enormously in size: the biggest flies about 28 times as many seat-kilometres as the smallest. Twenty incidents in fifteen years mean something very different for a giant carrier than for a small one. The model must compare rates, incidents per unit of flying, and the GLM does this with an offset. If airline $i$ has exposure $E_i$ (billions of seat-kilometres) and an incident rate $r_i$ per billion, then its expected count is ${\lambda }_i=E_i{r}_i$, so

The term $\log E_i$ is part of the linear predictor with its coefficient fixed at 1: an offset. Everything else models the rate. Our predictor $x_i$ is the airline's standardised log incident rate in 1985–99 (adding 0.5 to each count so that airlines with zero incidents get a finite rate). One caveat: FiveThirtyEight recorded capacity only for 2014, so we use the same exposure for both periods. That is a measurement assumption worth stating in a paper.

m_pois <- brm(incidents_00_14 ~ 1 + past_z + offset(log(exposure)), data = air,
              family = poisson(),
              prior = c(prior(normal(-5, 2), class = Intercept),   # about 7 incidents per 1,000 bn seat-km
                        prior(normal(0, 1), class = b)),
              seed = 4, refresh = 0, file = "fits/ch4_pois")
fixef(m_pois) |> round(3)
exp(fixef(m_pois)["past_z", -2]) |> round(2)    # rate ratio per standard deviation of past rate

          Estimate Est.Error   Q2.5  Q97.5
Intercept   -5.479     0.067 -5.614 -5.351
past_z       0.474     0.070  0.336  0.615
Estimate     Q2.5    Q97.5 
    1.61     1.40     1.85 

The prior for the intercept, Normal(−5, 2), says the typical airline has $e^{-5}\approx 0.007$ incidents per billion seat-kilometres, about 7 per thousand billion (roughly the 15-year exposure of a mid-sized airline), with a wide range of plausible values. The posterior says that past record matters: an airline whose 1985–99 incident rate was one standard deviation above average had, in 2000–14, an incident rate about 1.6 times higher. Past safety records do predict future ones. But how sure should we be? That depends on whether the Poisson model is right.

4.9 Overdispersion and the negative binomial

The Poisson distribution has a rigid property: its variance equals its mean. Real counts are almost always more variable than that, because the units differ in ways the model does not capture (airlines differ in fleet age, regulators, training, weather and routes, beyond what their past record reveals). Extra variation of this kind is called overdispersion, and ignoring it makes a model too confident. A posterior predictive check shows whether we have it.

m_nb <- brm(incidents_00_14 ~ 1 + past_z + offset(log(exposure)), data = air,
            family = negbinomial(),
            prior = c(prior(normal(-5, 2), class = Intercept), prior(normal(0, 1), class = b),
                      prior(exponential(0.2), class = shape)),
            seed = 4, refresh = 0, file = "fits/ch4_nb")
set.seed(8)
y <- air$incidents_00_14
chisq <- function(fit, model) {             # Pearson discrepancy, for observed and replicated data
  mu   <- posterior_epred(fit, ndraws = 1000)
  yrep <- posterior_predict(fit, ndraws = 1000)
  tibble(model, observed = rowSums(sweep(mu, 2, y)^2 / mu), replicated = rowSums((yrep - mu)^2 / mu))
}
ppc <- bind_rows(chisq(m_pois, "Poisson"), chisq(m_nb, "negative binomial"))
ppc |> group_by(model) |> summarise(`Pr(replicated >= observed)` = mean(replicated >= observed))
rbind(Poisson = fixef(m_pois)["past_z", ], `negative binomial` = fixef(m_nb)["past_z", ]) |> round(2)

# A tibble: 2 x 2
  model             `Pr(replicated >= observed)`
  <chr>                                    <dbl>
1 Poisson                                  0    
2 negative binomial                        0.589
                  Estimate Est.Error Q2.5 Q97.5
Poisson               0.47      0.07 0.34  0.62
negative binomial     0.52      0.11 0.32  0.75

ppc |> mutate(model = factor(model, levels = c("Poisson", "negative binomial"))) |>
  ggplot(aes(observed, replicated)) +
  geom_abline(linetype = "dashed", colour = harbour[["muted"]]) +
  geom_point(data = \(d) slice_head(d, n = 400, by = model), colour = harbour[["sky"]], alpha = 0.3, size = 0.9) +
  facet_wrap(~ model, scales = "free") +
  labs(x = "discrepancy for the real counts", y = "discrepancy for simulated counts")

Which statistic should the check use? Simple summaries such as the spread of the raw counts or the number of airlines with zero incidents turn out to be reproduced well by the Poisson model, because differences in exposure and past record already create a lot of variation between airlines. The statistic that targets overdispersion directly compares each airline's count with its own expected value: the ${\chi }^2$ discrepancy ${\sum }_i(y_i-{\mu }_i{)}^2/{\mu }_i$, computed under each posterior draw for the real data and for data simulated from that draw (Gelman, Meng and Stern, 1996). Under the Poisson model the real counts are far more dispersed than the model's own simulations in every one of the 1,000 draws. The negative binomial (also called gamma-Poisson) model fixes this by letting each airline's true rate vary around the regression line, following a gamma distribution, before the counts are drawn. It has one extra parameter, a shape $\varphi$: the variance of the counts is $\mu +{\mu }^2/\varphi$, so small $\varphi$ means a lot of extra variation and very large $\varphi$ gives back the Poisson.

The substantive conclusion survives, but the uncertainty changes. Under the negative binomial, the coefficient on past record is about the same (0.52 against 0.47), but its 95% interval is about 1.5 times as wide. The Poisson model was not wrong about the direction; it was overconfident about the precision, because it treated every airline's count as if its rate were known exactly given the past record. In a paper, the Poisson intervals would have overstated the evidence, and with a weaker effect they could have turned noise into a "finding".

Finally, look at what the model says about individual airlines. The figure plots each airline's incident rate in the second period against its rate in the first.

grid_a <- tibble(past_z = seq(min(air$past_z), max(air$past_z), length.out = 50), exposure = 1000)
fut <- posterior_epred(m_nb, newdata = grid_a)              # expected count per 1,000 bn seat-km
grid_a <- grid_a |>
  mutate(past_rate = exp(past_z * sd(air$past) + mean(air$past)) * 1000,
         rate = colMeans(fut), lo = apply(fut, 2, quantile, 0.025), hi = apply(fut, 2, quantile, 0.975))
air <- air |> mutate(past_rate = exp(past) * 1000,
                     future_rate = (incidents_00_14 + 0.5) / exposure * 1000)
ggplot(air, aes(past_rate, future_rate)) +
  geom_abline(linetype = "dashed", colour = harbour[["muted"]]) +
  geom_ribbon(data = grid_a, aes(past_rate, ymin = lo, ymax = hi), inherit.aes = FALSE,
              fill = harbour[["sky"]], alpha = 0.25) +
  geom_line(data = grid_a, aes(past_rate, rate), colour = harbour[["sky"]], linewidth = 1) +
  geom_point(colour = harbour[["ink"]], alpha = 0.7) +
  geom_text(data = filter(air, past_rate > 40 | future_rate > 25),
            aes(label = str_remove(airline, "\\*")), nudge_y = 0.12, size = 3, colour = harbour[["ink"]]) +
  scale_x_log10() + scale_y_log10() +
  labs(x = "incidents per 1,000 bn seat-km, 1985–99", y = "incidents per 1,000 bn seat-km, 2000–14")

The model's line is much flatter than the diagonal. Airlines with very bad records in 1985–99 did much better in 2000–14, and airlines with spotless records did somewhat worse: this is Galton's regression to the mean (Section 3.1) again. Aeroflot had 76 incidents in the first period and 6 in the second. Part of an extreme record is persistent (the slope is clearly positive) and part is luck, and the model separates the two. The same logic will reappear in Chapter 6 as shrinkage, where a multilevel model estimates each airline's rate by pulling its own noisy record toward the average of all airlines by just the right amount.

4.10 Exercises

Each exercise has a hint and an answer button.

Further reading and references

1. Angwin, J., Larson, J., Mattu, S. and Kirchner, L. (2016). Machine bias. ProPublica, 23 May 2016. Data and code: github.com/propublica/compas-analysis.
2. Berkson, J. (1944). Application of the logistic function to bio-assay. Journal of the American Statistical Association, 39(227), 357–365.
3. Bortkiewicz, L. von (1898). Das Gesetz der kleinen Zahlen. Leipzig: Teubner.
4. Burgess, E. W. (1928). Factors determining success or failure on parole. In A. A. Bruce et al., The Workings of the Indeterminate-Sentence Law and the Parole System in Illinois, 205–249. Springfield: Illinois State Parole Board.
5. Chouldechova, A. (2017). Fair prediction with disparate impact: a study of bias in recidivism prediction instruments. Big Data, 5(2), 153–163.
6. Collins, G. S., Reitsma, J. B., Altman, D. G. and Moons, K. G. M. (2015). Transparent reporting of a multivariable prediction model for individual prognosis or diagnosis (TRIPOD). Annals of Internal Medicine, 162(1), 55–63.
7. Dieterich, W., Mendoza, C. and Brennan, T. (2016). COMPAS Risk Scales: Demonstrating Accuracy Equity and Predictive Parity. Northpointe Inc.
8. Dressel, J. and Farid, H. (2018). The accuracy, fairness, and limits of predicting recidivism. Science Advances, 4(1), eaao5580.
9. Gelman, A., Meng, X.-L. and Stern, H. (1996). Posterior predictive assessment of model fitness via realized discrepancies. Statistica Sinica, 6(4), 733–760.
10. Gelman, A., Hill, J. and Vehtari, A. (2020). Regression and Other Stories. Cambridge University Press.
11. Kleinberg, J., Mullainathan, S. and Raghavan, M. (2017). Inherent trade-offs in the fair determination of risk scores. Proceedings of Innovations in Theoretical Computer Science, 43:1–43:23.
12. Liang, W., Yuksekgonul, M., Mao, Y., Wu, E. and Zou, J. (2023). GPT detectors are biased against non-native English writers. Patterns, 4(7), 100779.
13. Mitchell, S., Potash, E., Barocas, S., D'Amour, A. and Lum, K. (2021). Algorithmic fairness: choices, assumptions, and definitions. Annual Review of Statistics and Its Application, 8, 141–163.
14. Poisson, S. D. (1837). Recherches sur la probabilité des jugements en matière criminelle et en matière civile. Paris: Bachelier.
15. Silver, N. (2014). Should travelers avoid flying airlines that have had crashes in the past? FiveThirtyEight, 18 July 2014. Data: github.com/fivethirtyeight/data/tree/master/airline-safety.
16. Verhulst, P.-F. (1838). Notice sur la loi que la population suit dans son accroissement. Correspondance mathématique et physique, 10, 113–121.