4.6 Counts and the GLM family

Linear and logistic regression are two members of one family. Each member fits a straight line on a transformed scale (the link) and pairs it with a distribution that suits the outcome. For counts, the scale is the log and the distribution is the Poisson, or the negative binomial when units differ.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from matplotlib.axes import Axes
from scipy import stats
from statsmodels.base.model import LikelihoodModelResults

from statly import style

style.use()
rng = np.random.default_rng(1)

health = sm.datasets.randhie.load_pandas().data  # 20,190 people in the RAND health experiment
cancer = (
    sm.datasets.cancer.load_pandas().data
)  # breast cancer cases and population, 301 counties
engel = (
    sm.datasets.engel.load_pandas().data
)  # income and food spending, Belgian households (1857)

4.6.1 Poisson regression: a line on the log scale

A count, such as doctor visits per year, cannot be negative, and its spread grows with its mean (Chapter 2.2). A straight line would predict negative visits for healthy people. Poisson regression fits the straight line to the log of the expected count, μ:

$$\log \mu = \beta_0 + \beta_1 x \qquad \Longleftrightarrow \qquad \mu = e^{\beta_0 + \beta_1 x}$$

and assumes the counts are Poisson around μ. The log is the link. Every step in x adds $\beta_1$ to the log, which multiplies the expected count by $e^{\beta_1}$: the rate ratio.

In the RAND study, each person has a chronic disease index: higher means more or worse chronic conditions. Do visits rise with it?

poisson_one = smf.poisson("mdvis ~ disea", health).fit(disp=0)
health["bin"] = pd.qcut(health.disea, 10, labels=False, duplicates="drop")
by_bin = health.groupby("bin").agg(
    index=("disea", "mean"), visits=("mdvis", "mean"), n=("mdvis", "size")
)
grid = np.linspace(0, 40, 200)
curve = np.exp(poisson_one.params.Intercept + poisson_one.params.disea * grid)


def plot_count_scale(ax: Axes) -> None:
    """Average visits per bin with the fitted curve on the count scale."""
    ax.scatter(by_bin["index"], by_bin.visits, s=40, color=style.BLUE)
    ax.plot(grid, curve, color=style.CORAL)
    ax.set(
        title="Count scale: a curve",
        xlabel="chronic disease index",
        ylabel="average visits per year",
    )


def plot_log_scale(ax: Axes) -> None:
    """The same on the log scale: a line."""
    ax.scatter(by_bin["index"], np.log(by_bin.visits), s=40, color=style.GREEN)
    ax.plot(grid, np.log(curve), color=style.CORAL)
    style.note(
        ax,
        f"rate ratio per point: {np.exp(poisson_one.params.disea):.3f}",
        x=0.03,
        y=0.98,
        color=style.INK,
    )
    ax.set(
        title="Log scale: a line",
        xlabel="chronic disease index",
        ylabel="log of average visits",
    )


fig, axes = style.subplots(2, height=3.1, width=9.5)
plot_count_scale(axes[0])
plot_log_scale(axes[1])
fig.tight_layout()
plt.show()

Figure 4.6.1: Doctor visits against the chronic disease index for 20,190 people, grouped into ten bins by the index. Left: the average visits per bin on the count scale, with the Poisson regression curve (red): it bends upward. Right: the same averages on the log scale, where the model is a straight line.

Each point of the index multiplies the expected visits by about 1.04: 4% more. Ten points mean $1.04^{10} \approx 1.5$ times as many visits. As for logistic regression (Chapter 4.4), the model is fitted by maximum likelihood, here with the Poisson likelihood of Chapter 2.2.

4.6.2 Exposure: counts per unit

Counts often come from units of different size: deaths in counties of different population, accidents in years of different traffic, defects in batches of different volume. More population means more cases, and that alone says nothing interesting. What we want is the rate: cases per person. Poisson regression handles this with an offset: the log of the exposure enters the model with its coefficient fixed at 1:

$$\log \mu = \log(\text{population}) + \beta_0 + \beta_1 x \qquad \Longleftrightarrow \qquad \frac{\mu}{\text{population}} = e^{\beta_0 + \beta_1 x}$$

so the coefficients describe the rate. For breast cancer cases in 301 counties, the data show both why the offset is right and why small units need care:

free_slope = smf.glm("cancer ~ np.log(population)", cancer, family=sm.families.Poisson()).fit()
with_offset = smf.glm(
    "cancer ~ 1", cancer, family=sm.families.Poisson(), offset=np.log(cancer.population)
).fit()
overall = np.exp(with_offset.params.Intercept)  # cases per person


def plot_log_log(ax: Axes) -> None:
    """Cases against population on log scales, with the fitted and the proportional line."""
    shown = cancer[cancer.cancer > 0]
    ax.scatter(shown.population, shown.cancer, s=10, color=style.BLUE, alpha=0.6)
    pop = np.geomspace(400, 100_000, 50)
    slope = free_slope.params["np.log(population)"]
    ax.plot(
        pop,
        np.exp(free_slope.params.Intercept) * pop**slope,
        color=style.CORAL,
        label=f"fitted slope {slope:.2f}",
    )
    ax.plot(pop, overall * pop, color=style.INK, ls="--", label="slope 1 (offset)")
    ax.set_xscale("log")
    ax.set_yscale("log")
    ax.legend(fontsize=8)
    ax.set(
        title="Cases grow with population",
        xlabel="population (log scale)",
        ylabel="cases (log scale)",
    )


def plot_funnel(ax: Axes) -> None:
    """Rates per 1,000 against population, with Poisson 95% limits."""
    rate = cancer.cancer / cancer.population * 1000
    pop = np.geomspace(cancer.population.min(), cancer.population.max(), 200)
    low = stats.poisson.ppf(0.025, overall * pop) / pop * 1000
    high = stats.poisson.ppf(0.975, overall * pop) / pop * 1000
    ax.fill_between(pop, low, high, color=style.GREY, alpha=0.2, lw=0)
    ax.scatter(cancer.population, rate, s=10, color=style.BLUE, alpha=0.6)
    ax.axhline(overall * 1000, color=style.CORAL, lw=1.4)
    outside = (
        (rate < np.interp(cancer.population, pop, low))
        | (rate > np.interp(cancer.population, pop, high))
    ).mean()
    style.note(ax, f"outside the funnel: {outside:.0%}", x=0.55, y=0.98, color=style.INK)
    ax.set_xscale("log")
    ax.set(
        title="Rates: a funnel",
        xlabel="population (log scale)",
        ylabel="cases per 1,000 people",
    )


fig, axes = style.subplots(2, height=3.3, width=9.5)
plot_log_log(axes[0])
plot_funnel(axes[1])
fig.tight_layout()
plt.show()

Figure 4.6.2: Breast cancer cases in 301 counties. Left: cases against population, both on log scales; the fitted slope (red) is 0.99, almost exactly the 1 an offset assumes (dashed): cases grow in proportion to population. Right: a funnel plot of each county’s rate against its population, with the range a Poisson count around the overall rate would cover 95% of the time. Small counties scatter widely by chance alone.

The funnel is worth remembering. Ranking units by their rate puts the smallest ones at both the top and the bottom of the list, not because they are special but because small counts are noisy. The same is true of hospitals, schools and towns. More counties fall outside the funnel than the 5% chance alone would put there: the counties really do differ somewhat, which is the subject of the next section.

4.6.3 Overdispersion: the negative binomial

The Poisson assumes the variance equals the mean. Real counts often vary more, because the units differ in ways the predictors do not capture (Chapter 2.2). A quick check is the dispersion: the sum of squared Pearson residuals, $\sum (y_i - \hat\mu_i)^2 / \hat\mu_i$, divided by the residual degrees of freedom. Near 1 is fine; much above 1 is overdispersion. For the cancer counts it is 2.6; for the doctor visits, about 6.

With overdispersion the Poisson estimates are still sensible, but its standard errors are too small. The negative binomial model adds one parameter for the extra variation, $\text{Var}(Y) = \mu + \alpha\mu^2$, and gives honest intervals:

FORMULA = "mdvis ~ lncoins + physlm + disea + hlthg + hlthf + hlthp"
poisson = smf.poisson(FORMULA, health).fit(disp=0)
negbin = smf.negativebinomial(FORMULA, health).fit(disp=0, maxiter=200)
COUNT_LABELS = {
    "lncoins": "paying more of the bill (per log unit)",
    "physlm": "has a physical limitation",
    "disea": "chronic disease index (per point)",
    "hlthg": "health good (vs excellent)",
    "hlthf": "health fair (vs excellent)",
    "hlthp": "health poor (vs excellent)",
}


def offset_ratios(ax: Axes, model: LikelihoodModelResults, shift: float, color: str) -> None:
    """Rate ratios with intervals, shifted up or down a little to sit next to another model."""
    intervals = np.exp(model.conf_int())
    for i, name in enumerate(COUNT_LABELS):
        low, high = intervals.loc[name]
        ax.plot([low, high], [i + shift, i + shift], color=color, lw=2.5)
        ax.plot(np.exp(model.params[name]), i + shift, "o", ms=5, color=color)


fig, ax = plt.subplots(figsize=(8, 3.4))
offset_ratios(ax, poisson, -0.14, style.BLUE)
offset_ratios(ax, negbin, 0.14, style.CORAL)
ax.axvline(1, color=style.INK, lw=1)
ax.set_xscale("log")
ticks = [0.8, 0.9, 1, 1.25, 1.5, 2]
ax.set_xticks(ticks, [str(t) for t in ticks])
ax.minorticks_off()
ax.set_yticks(range(len(COUNT_LABELS)), list(COUNT_LABELS.values()))
ax.invert_yaxis()
ax.plot([], [], color=style.BLUE, lw=2.5, label="Poisson")
ax.plot([], [], color=style.CORAL, lw=2.5, label="negative binomial")
ax.legend(fontsize=8, loc="center right")
ax.set_xlabel("rate ratio: times as many visits (log scale)")
plt.show()

Figure 4.6.3: Rate ratios for doctor visits per year, with 95% intervals, from a Poisson model (blue) and a negative binomial model (red) with the same predictors. The estimates agree, but the Poisson intervals are about half as wide: too narrow, because the counts vary more than a Poisson allows.

The RAND experiment randomly gave people insurance plans that differed in how much of the bill they paid themselves. The rate ratio for cost sharing is about 0.94 per log unit of the share they paid. Going from a free plan to paying 95% of the bill is about 4.6 log units, so 0.94 multiplied 4.6 times: about a quarter fewer visits. Because the plans were assigned at random, this is a causal effect, one of the best-known results in health economics.

4.6.4 Too many zeros

Even the negative binomial can miss one feature of counts: an excess of zeros, from units that never have the event at all (people who never go to a doctor, fields where a plant cannot grow). Compare the share of zeros each model expects with the share observed:

k = np.arange(0, 16)
observed = np.bincount(health.mdvis.astype(int), minlength=16)[:16] / len(health)
mu_poisson = poisson.predict(health).to_numpy()
mu_negbin = negbin.predict(health).to_numpy()
size = 1 / negbin.params.alpha
expected_poisson = np.array([stats.poisson.pmf(j, mu_poisson).mean() for j in k])
expected_negbin = np.array(
    [stats.nbinom.pmf(j, size, size / (size + mu_negbin)).mean() for j in k]
)


def plot_zeros(ax: Axes) -> None:
    """Observed and expected shares of each count."""
    ax.bar(k, observed, color=style.GOLD, alpha=0.6, label="observed")
    ax.plot(k, expected_poisson, "o-", ms=4, color=style.BLUE, label="Poisson model")
    ax.plot(k, expected_negbin, "o-", ms=4, color=style.CORAL, label="negative binomial model")
    ax.yaxis.set_major_formatter(style.PERCENT)
    ax.legend(fontsize=8)
    ax.set(xlabel="doctor visits per year", ylabel="share of people")


fig, ax = plt.subplots(figsize=(8, 3.2))
plot_zeros(ax)
plt.show()

Figure 4.6.4: Share of people with 0, 1, 2, … doctor visits (bars), with what the fitted Poisson model (blue) and negative binomial model (red) expect, averaged over everyone. The Poisson expects far too few zeros; the negative binomial gets them close.

Here the negative binomial is close enough. When it is not, zero-inflated models mix a group that is always zero with an ordinary count model, and hurdle models fit “any visit at all?” (a logistic model) separately from “how many, given at least one?” (a count model that starts at 1).

4.6.5 Positive amounts: gamma regression

Amounts of money, waiting times and concentrations are positive and skewed, and their spread grows with their level. Logging them and running a linear regression (Chapter 4.3) often works. A gamma regression with a log link models the amounts directly: the log of the expected amount is a straight line, and the noise follows the gamma distribution of Chapter 2.1, whose standard deviation grows in proportion to the mean.

gamma = smf.glm(
    "foodexp ~ np.log(income)", engel, family=sm.families.Gamma(sm.families.links.Log())
).fit()
plain = smf.ols("foodexp ~ income", engel).fit()
incomes = pd.DataFrame({"income": np.linspace(engel.income.min(), engel.income.max(), 200)})


def plot_gamma(ax: Axes) -> None:
    """Food spending against income with a straight line and a gamma regression curve."""
    ax.scatter(engel.income, engel.foodexp, s=10, color=style.BLUE, alpha=0.5)
    ax.plot(
        incomes.income, plain.predict(incomes), color=style.INK, ls="--", label="straight line"
    )
    b = gamma.params["np.log(income)"]
    ax.plot(
        incomes.income,
        gamma.predict(incomes),
        color=style.CORAL,
        label=f"gamma, log link: spending ∝ income^{b:.2f}",
    )
    ax.legend(fontsize=8, loc="upper left")
    ax.set(xlabel="household income", ylabel="food spending")


fig, ax = plt.subplots(figsize=(8, 3.3))
plot_gamma(ax)
plt.show()

Figure 4.6.5: Food spending of the 235 households against income. A gamma regression with a log link on log income (red) is a power curve, spending ∝ income^b, with b about 0.85: the elasticity of Chapter 4.3, estimated on the original spending scale. Its noise grows with the level, as the data’s does. A straight line (dashed) assumes equal noise everywhere.

4.6.6 The family

Linear, logistic, Poisson and gamma regression are members of one family, the generalised linear models (GLMs). Each has three parts: a linear predictor $\eta = \beta_0 + \beta_1 x_1 + \dots$, a link g that connects it to the average outcome, $g(\mu) = \eta$, and a distribution for the noise.

Outcome	Distribution	Link $g(\mu)$	Coefficients become	Chapter
a number	normal	none: μ	differences	4.1, 4.2
yes/no	Bernoulli (binomial)	log-odds: log(μ/(1 − μ))	odds ratios	4.4
a count	Poisson	log μ	rate ratios	4.5
a count, units differ	negative binomial	log μ	rate ratios	4.5
a positive amount	gamma	log μ	ratios	4.5

All are fitted by maximum likelihood, compared with deviance (Chapter 4.4) and likelihood ratio tests, and checked with residuals. Two kinds of residual are common: Pearson residuals, $(y - \hat\mu)/\sqrt{\text{Var}(\hat\mu)}$, and deviance residuals, each case’s share of the deviance with the sign of $y - \hat\mu$. Both should look like structureless noise in the plots of Chapter 4.3. Everything from Chapters 4.1 to 4.3 carries over: holding other predictors fixed, interactions, confounding, and influential points.

4.6.7 In one picture

• Poisson regression: a straight line on the log scale, $\log \mu = \beta_0 + \beta_1 x$; coefficients become rate ratios.
• Counts from units of different size need an offset, log(exposure), so the model describes rates. Small units scatter widely: look at the funnel.
• Check the dispersion. Much above 1 means the Poisson intervals are too narrow: use the negative binomial. Check the zeros too.
• Gamma regression suits positive, skewed amounts.
• All are GLMs: a linear predictor, a link and a distribution.

Next: data that come in clusters.