From Poisson to the negative binomial

You fit a Poisson GLM to count data, check the diagnostics and find that the variance of the data is three, five, ten times the mean. The textbook says: “there is overdispersion, use the negative binomial”. This article explains why — not as a recipe but as the consequence of an argument. The central idea is to write the count in two stages: first nature draws a rate \(Z\) for each unit; then, given the rate, the count is Poisson. With two tools — the iterated expectation \(\mathbb E[\mathbb E[Y \mid Z]]\) and the law of total variance — and one integral, the negative binomial appears on its own. At the end we fit both models to the same dataset and measure, by simulation, what ignoring overdispersion costs. Every number comes from a reproducible script.

§1The starting point: the Poisson GLM

Count data — doctor visits per year, claims per policy, defects per part, goals per match — are non-negative integers. The simplest model for them is the Poisson:

In a Poisson GLM each observation \(i\) has its own mean, linked to the covariates through the logarithm:

The Poisson has one property that drives everything that follows: the variance equals the mean,

This is called equidispersion. It is a strong restriction: the model has no free parameter for the variance. Once the mean is chosen, the variance is decided. In practice, real data almost always have a variance larger than the mean. That is overdispersion.

§2How overdispersion shows up

If model (2) is right, the Pearson residuals \(r_i = (y_i - \hat\mu_i)/\sqrt{\hat\mu_i}\) have variance close to 1, and the sum of their squares is close to the degrees of freedom. Hence the dispersion statistic:

where \(p\) is the number of coefficients. Under the Poisson model, \(\hat\phi \approx 1\). A \(\hat\phi\) well above 1 signals overdispersion (the deviance-to-df ratio tells the same story). In the example we analyse in §11 — 500 counts with one covariate — the Poisson fit gives \(\hat\phi = 3.20\): the variance of the data is about three times what the Poisson allows.

A plot helps more than a number. Sort the observations by fitted mean, split them into groups and, in each group, compare the mean with the sample variance. Under Poisson the points would sit around the line \(\operatorname{Var} = \mu\):

Every point is above the line, and further above as the mean grows. The solid curve, \(\mu + \mu^2/\theta\), is the variance of the negative binomial — we will get there in §9. First the important question: why do data behave like this?

§3Where overdispersion comes from: unobserved heterogeneity

The Poisson arises naturally when rare events occur independently, all with the same rate. That last condition is the one that usually fails. Think of doctor visits per year: two people with the same age, sex and income (the model's covariates) do not have the same propensity to see a doctor. One has a chronic condition, another lives far from the clinic, a third is more anxious. None of that is in the data. The model gives both the same \(\mu_i\), but their true rates differ.

The way to formalise this is to give each unit a random rate. Call \(Z\) the rate nature draws for a unit. The model now has two layers:

Conditionally on \(Z\), everything is still Poisson: variance equals mean. But we do not observe \(Z\); we only observe \(Y\). The distribution the data follow is the marginal of \(Y\) — the Poisson “mixed” over all possible rates. The whole question reduces to: what are the mean, the variance and the distribution of that marginal?

A convention for the Gamma: we use shape \(\alpha\) and rate \(\beta\), with density

(Some books and numpy use the scale \(1/\beta\) instead of the rate. The algebra is the same with \(\beta\) replaced by \(1/\text{scale}\).) Why Gamma and not some other distribution? For now, accept that it is flexible (positive values, adjustable skewness) — §8 gives the deeper reason.

§4Tool 1: the iterated expectation, \(\mathbb E[\mathbb E[Y \mid Z]]\)

The first tool is the law of total expectation (or iterated expectation):

Read it from the inside out. \(\mathbb E[Y \mid Z]\) is the mean of \(Y\) if we knew the rate \(Z\). Since \(Z\) is random, \(\mathbb E[Y \mid Z]\) is a random variable — a function of \(Z\). The outer expectation averages that function over all possible values of \(Z\). In words: the overall mean is the mean of the conditional means.

Proof (for discrete \(Y\) and continuous \(Z\), our case). Write \(p(y \mid z) = \Pr(Y = y \mid Z = z)\). Then \(\mathbb E[Y \mid Z = z] = \sum_y y\, p(y \mid z)\) and

The second line swaps sum and integral (allowed because every term is non-negative — Tonelli's theorem). The next step uses the law of total probability: \(\Pr(Y = y) = \int p(y \mid z) f_Z(z)\,dz\). Keep that integral in mind — it comes back in §7 to give the whole distribution of \(Y\).

§5Tool 2: the law of total variance

The mean alone does not reveal overdispersion; we need the variance. The law of total variance says:

Proof. Start from \(\operatorname{Var}(Y) = \mathbb E[Y^2] - (\mathbb E[Y])^2\) and apply (7) to both terms — the iterated expectation holds for any function of \(Y\), including \(Y^2\):

The second line uses, conditionally on \(Z\), the identity \(\mathbb E[Y^2 \mid Z] = \operatorname{Var}(Y \mid Z) + (\mathbb E[Y \mid Z])^2\). In the third, the key is to recognise the expression in braces as the variance of the random variable \(\mathbb E[Y \mid Z]\).

The interpretation is the heart of the argument. The variance of \(Y\) has two sources: the noise within each unit (even knowing the rate, the count varies — that is the Poisson noise) and the differences between units (rates vary from unit to unit). The Poisson GLM only sees the first.

§6Applying both laws: overdispersion appears on its own

In model (5), given \(Z\), \(Y\) is Poisson with mean \(Z\). So

The mean. By the iterated expectation (7):

The variance. By total variance (9):

Compare with the mean:

This is the main result of the article, and it is worth reading slowly. The excess variance is exactly the variance of the rates. If all units had the same rate (\(\operatorname{Var}(Z) = 0\)) we would be back to equidispersion. Any unobserved heterogeneity, however small, pushes the variance above the mean. Note too that (14) used nothing specific to the Gamma beyond \(\operatorname{Var}(Z)\): every mixture of Poissons is overdispersed. The Gamma only comes in to give the resulting distribution a name.

You can check this without trusting any algebra: just simulate the two-stage experiment. With \(\alpha = 2\) and \(\beta = 0.5\), the formulas predict \(\mathbb E[Y] = 4\) and \(\operatorname{Var}(Y) = 4 + 8 = 12\):

import numpy as np
rng = np.random.default_rng(20261109)
alpha, beta = 2.0, 0.5
z = rng.gamma(shape=alpha, scale=1/beta, size=1_000_000)   # stage 1: one rate per unit
y = rng.poisson(z)                                          # stage 2: Poisson given the rate
print(y.mean(), y.var())                                    # ≈ 4 and ≈ 12

\(\alpha = 2\), \(\beta = 0.5\) (rate), \(p = \beta/(1+\beta) = 1/3\). The variance of \(Y\) is the Poisson part (4) plus the variance of \(Z\) (8): heterogeneity triples the variance without touching the mean.

§7The marginal distribution: the integral that gives the negative binomial

Mean and variance already explain overdispersion, but to fit a model by maximum likelihood we need the distribution of \(Y\). It comes from the law of total probability — the same integral that appeared inside proof (8):

Step 1 — take out of the integral whatever does not depend on \(z\) and collect powers:

Step 2 — recognise a Gamma kernel. Up to a constant, the integrand is the density of a \(\text{Gamma}(y + \alpha,\; \beta + 1)\). Since every density integrates to 1,

Step 3 — substitute and rearrange:

Step 4 — name it. Let \(p = \beta/(\beta + 1)\), so that \(1 - p = 1/(\beta+1)\), and \(r = \alpha\). Then

This is the negative binomial \(\text{NB}(r, p)\). When \(r\) is an integer, \(\Gamma(y + r)/(\Gamma(r)\,y!) = \binom{y + r - 1}{y}\) and (19) is the classical distribution of the number of failures before the \(r\)-th success in Bernoulli trials with probability \(p\). Here \(r = \alpha\) can be any positive real; the Gamma-function form is the natural extension. The name “negative” comes from the binomial series with a negative exponent, \((1 - q)^{-r} = \sum_y \binom{y + r - 1}{y} q^y\), which also shows at once that (19) sums to 1.

Checking the moments with the NB formulas, \(\mathbb E[Y] = r(1-p)/p\) and \(\operatorname{Var}(Y) = r(1-p)/p^2\):

exactly (12) and (13). The two routes — moments through the iterated expectation and the distribution through the integral — arrive at the same place. And the simulation agrees, point by point:

A shortcut: the probability generating function

The iterated expectation also gives a route that avoids integrating densities. The probability generating function of \(Y\) is \(G_Y(s) = \mathbb E[s^{Y}]\). Conditioning on \(Z\) and using the Poisson one, \(\mathbb E[s^{Y} \mid Z] = e^{Z(s-1)}\):

where \(M_Z(t) = (\beta/(\beta - t))^{\alpha}\) is the Gamma moment generating function. The last expression is the generating function of \(\text{NB}(\alpha, p)\). Same result, with \(\mathbb E[\mathbb E[\,\cdot \mid Z]]\) doing all the work.

§8Why the Gamma, specifically?

§6 showed that any distribution for \(Z\) produces overdispersion. The Gamma is chosen because it is conjugate to the Poisson: the product “Poisson × Gamma” inside integral (15) is another Gamma, so the integral has a closed form. With a lognormal \(Z\), for instance, the Poisson-lognormal model is perfectly reasonable, but \(\Pr(Y = y)\) only comes out by numerical integration.

Conjugacy has a bonus. By Bayes' rule, the distribution of a unit's rate after seeing its count is

with mean

The best estimate of a unit's rate is a weighted average of what it showed (\(y\)) and the population mean (\(\alpha/\beta\)). This is the credibility principle in insurance: a policyholder with many claims pays more, but not in proportion to a short history.

§9From the mixture model to the negative binomial GLM

So far there was a single population. In a GLM each observation has its own mean \(\mu_i = \exp(\mathbf x_i^\top \boldsymbol\beta)\). To accommodate that, write unit \(i\)'s rate as the model mean times an unobserved multiplicative effect:

Shape and rate both equal to \(\theta\) make \(\mathbb E[\varepsilon_i] = 1\). That is not arbitrary: if \(\varepsilon_i\) had mean \(c \ne 1\), \(c\) would be absorbed into the intercept \(\beta_0\) and the two could not be separated. Fixing the mean at 1 keeps \(\mu_i\) as the mean of \(Y_i\). Note that \(Z_i = \mu_i \varepsilon_i \sim \text{Gamma}(\theta, \theta/\mu_i)\): this is model (5) with \(\alpha = \theta\) and \(\beta = \theta/\mu_i\). Substituting into (12), (13) and (19):

The variance, again by total variance: \(\mathbb E[\mu_i\varepsilon_i] + \operatorname{Var}(\mu_i\varepsilon_i) = \mu_i + \mu_i^2/\theta\). The term \(\mu_i^2/\theta\) is the excess; it grows with the square of the mean, which is why the solid curve in the §2 figure pulls away from the line. This version is called NB2 (variance quadratic in the mean) and is the default in almost all software.

The parameter \(\theta\) controls heterogeneity. Many texts use \(\kappa = 1/\theta\) (statsmodels' alpha; R reports theta). Small \(\theta\) = very different rates across units = a lot of overdispersion. And as \(\theta \to \infty\), \(\operatorname{Var}(\varepsilon_i) \to 0\), all units share the same rate and we should recover the Poisson. Indeed, rewriting (25):

The Poisson is a limiting case of the negative binomial. That nesting gives us a test in §11.

§10Is the negative binomial a GLM?

A GLM requires the distribution to belong to the exponential family. With \(\theta\) fixed, the NB does. Take the log of (25):

With canonical parameter \(\vartheta = \log(\mu/(\mu+\theta))\) we have \(1 - e^{\vartheta} = \theta/(\mu + \theta)\), and the form is \(y\vartheta - b(\vartheta) + c(y)\) with \(b(\vartheta) = -\theta \log(1 - e^{\vartheta})\). The derivatives of \(b\) return the mean and the variance, as they must in an exponential family:

The canonical link would be \(\log(\mu/(\mu+\theta))\), which depends on \(\theta\) and is hard to interpret; in practice one uses the log link, \(\log\mu_i = \mathbf x_i^\top\boldsymbol\beta\), the same as the Poisson. Coefficients then read the same way in both models: \(e^{\beta_j}\) is the ratio of means for a one-unit increase in \(x_j\).

Estimation

For fixed \(\theta\), \(\boldsymbol\beta\) comes from the usual GLM algorithm, iteratively reweighted least squares (IRLS). With the log link, \(d\mu/d\eta = \mu\), and the weights and score equations are

Compare with the Poisson, whose score is \(\sum_i (y_i - \mu_i)\,\mathbf x_i = \mathbf 0\) and whose weights are \(w_i = \mu_i\). The NB gives less weight to high-mean observations, precisely the noisiest ones. For \(\theta\), maximise the log-likelihood with \(\boldsymbol\beta\) fixed and alternate the two steps until convergence (this is what R's MASS::glm.nb does). Since \(\mathbb E[\partial^2\ell/\partial\boldsymbol\beta\,\partial\theta] = 0\), the estimates of \(\boldsymbol\beta\) and \(\theta\) are asymptotically independent, and the covariance of \(\hat{\boldsymbol\beta}\) is \((\mathbf X^\top \mathbf W \mathbf X)^{-1}\) with weights (29).

What exactly the Poisson gets wrong

The Poisson score equation has expectation zero whenever the mean is right, \(\mathbb E[Y_i] = \mu_i\) — whatever the variance. That is why \(\hat{\boldsymbol\beta}_{\text{Poisson}}\) remains consistent under overdispersion. The error is in the uncertainty. The Poisson GLM reports \((\mathbf X^\top \operatorname{diag}(\mu_i)\,\mathbf X)^{-1}\), which assumes \(\operatorname{Var}(Y_i) = \mu_i\). The estimator's true covariance is the “sandwich”

with \(\mathbf M = \operatorname{diag}(\mu_i)\), larger than the reported one. There are three ways to fix it:

1. Quasi-Poisson: assumes \(\operatorname{Var}(Y_i) = \phi\,\mu_i\) and multiplies the standard errors by \(\sqrt{\hat\phi}\). Keeps the Poisson estimates. Simple, but the shape of the variance (proportional to the mean) may be wrong, and there is no likelihood (no AIC, no predicted probabilities).
2. Robust (sandwich) standard errors: estimate (30) without assuming a form for the variance. Good with large \(n\).
3. Negative binomial: models the variance explicitly through (25). Full likelihood, probabilities for every count, and more efficient when the model is right.

§11A complete example

We simulate \(n = 500\) counts exactly along the path of (24): one covariate \(x_i \sim \mathcal N(0, 1)\), \(\log\mu_i = 1 + 0.5\,x_i\), an effect \(\varepsilon_i \sim \text{Gamma}(1.5,\ 1.5)\) for each unit and \(Y_i \mid \varepsilon_i \sim \text{Poisson}(\mu_i\varepsilon_i)\). Since the simulation is ours, we know the answer: \(\beta_0 = 1\), \(\beta_1 = 0.5\), \(\theta = 1.5\). Then we fit Poisson, quasi-Poisson and negative binomial:

n = 500 counts simulated with \(\log\mu_i = 1 + 0.5\,x_i\) and \(\theta = 1.5\). Standard errors in brackets. Sample mean of \(y\): 2.89; variance: 11.49; zeros: 124. Quasi-Poisson has the same estimates as Poisson and only multiplies the standard errors by \(\sqrt{\hat\phi}\).

Read the table row by row:

• Point estimates nearly coincide (\(\hat\beta_1 = 0.463\) for Poisson, \(0.478\) for NB). As expected: both score equations are unbiased when the mean is right.
• Standard errors do not. Poisson's for \(\beta_1\) is \(0.028\); the NB's is \(0.053\) — almost double. The Poisson 95% interval, \([0.407,\ 0.518]\), contains the true value 0.5, but it is little more than half as wide as the NB interval, \([0.375,\ 0.581]\): it promises a precision the data do not have. It worked in this sample; §12 shows how often it fails.
• The dispersion diagnostic goes from 3.20 under Poisson to 1.02 under NB: once heterogeneity is modelled, Pearson residuals have variance 1 again.
• \(\hat\theta = 1.29\) (standard error 0.13), compatible with the true 1.5.
• AIC drops from 2594.9 to 2146.9 — a huge difference for a single extra parameter.

Likelihood-ratio test: \(H_0\colon\) Poisson

Since the Poisson is the \(\theta \to \infty\) limit of the NB (26), we can test \(H_0\colon \kappa = 1/\theta = 0\) with a likelihood ratio:

A detail many forget: \(\kappa = 0\) is on the boundary of the parameter space (\(\kappa\) cannot be negative). The null distribution of LR is then not \(\chi^2_1\) but a 50:50 mixture of a point mass at zero and a \(\chi^2_1\). In practice: halve the \(\chi^2_1\) p-value. Here it changes nothing — the p-value is below \(10^{-99}\) — but in borderline cases it does.

The frequencies

The most concrete check is to compare how many observations of each value the model expected with how many showed up:

The Poisson fails exactly where the theory says it would. The mixture spreads mass to both sides: low-rate units produce too many zeros, high-rate units produce too many large values. The Poisson, tied to \(\operatorname{Var} = \mu\), cannot do both and piles everything in the middle. Much of what is called “excess zeros” is just this — overdispersion — and the NB handles it without a zero-inflated model.

§12What ignoring overdispersion costs

One example is one sample. To measure how the methods behave, we repeat the experiment 2,000 times for each level of heterogeneity: generate data with \(n = 200\), fit all three models and check whether the 95% Wald interval for \(\beta_1\) contains the true value 0.5. A good method succeeds 95% of the time.

2,000 datasets per row, n = 200, \(\log\mu = 1 + 0.5x\). Each cell: coverage of the 95% Wald interval for \(\beta_1\) · mean width. Monte Carlo error of the coverage: ± 0.5 pp. The last column is the standard deviation of \(\hat\beta_1\) across datasets: what a good standard error should estimate.

Four conclusions:

1. With Poisson data (\(\theta = \infty\)), all three are right. The NB costs nothing here: \(\hat\theta\) simply goes to huge values and the model becomes Poisson. Using NB “just in case” is cheap.
2. Poisson deteriorates quickly with heterogeneity. With \(\theta = 2\) (rates with a 71% standard deviation around the mean), the “95%” interval covers 73%; with \(\theta = 0.5\), less than half the time. A 5% test would reject a true \(H_0\) more than 50% of the time.
3. Quasi-Poisson helps a lot, but not fully. It fixes the average size of the variance, not its shape: it assumes \(\operatorname{Var} \propto \mu\) when the truth is \(\mu + \mu^2/\theta\). With \(\theta = 0.5\) it covers 88.6%.
4. The NB stays close to 95% in every case and is also more efficient: the actual standard deviation of \(\hat\beta_1\) is smaller with the NB (0.116 versus 0.132 for Poisson at \(\theta = 0.5\)) because the weights (29) downweight the noisiest observations.

§13In practice: Python and R

import statsmodels.api as sm
import statsmodels.formula.api as smf

pois = smf.glm("y ~ x", data=df, family=sm.families.Poisson()).fit()
phi = pois.pearson_chi2 / pois.df_resid                 # dispersion (4): close to 1?

quasi = smf.glm("y ~ x", data=df, family=sm.families.Poisson()).fit(scale="X2")   # quasi-Poisson
robust = smf.glm("y ~ x", data=df, family=sm.families.Poisson()).fit(cov_type="HC0")  # sandwich (30)

nb = smf.negativebinomial("y ~ x", data=df).fit()      # NB2; reports alpha = 1/theta
lr = 2 * (nb.llf - pois.llf)                            # p-value: chi2(1) halved
print(nb.summary())

pois  <- glm(y ~ x, family = poisson, data = d)
sum(residuals(pois, type = "pearson")^2) / df.residual(pois)   # dispersion

quasi <- glm(y ~ x, family = quasipoisson, data = d)

library(MASS)
nb <- glm.nb(y ~ x, data = d)          # reports theta (= 1/alpha)
lr <- 2 * (logLik(nb) - logLik(pois))
pchisq(lr, df = 1, lower.tail = FALSE) / 2   # boundary test

Watch the parameterisations: R reports \(\theta\); statsmodels reports \(\alpha = 1/\theta\); some texts call \(1/\theta\) \(k\) or \(\kappa\). It is always the same model (25).

§14Caveats

• Get the mean right first. A missing important covariate, a missing interaction or a non-linear relationship also produce \(\hat\phi \gt 1\). Overdispersion can be a symptom of a wrong linear predictor; the NB “absorbs” the problem, but the coefficients stay wrong.
• Genuine excess zeros. If a separate process produces structural zeros (non-smokers smoke zero cigarettes, always), even the NB will not cope; zero-inflated or hurdle models separate the two processes. Compare zero frequencies, as in §11, before concluding.
• Dependence. Clustered observations (students in the same school, repeated measures on the same patient) produce overdispersion that is not individual heterogeneity. There the route is a mixed model (a random effect per group) or GEE.
• NB1 vs NB2. If the variance grows linearly with the mean, \(\operatorname{Var} = \mu(1 + \delta)\), NB1 or quasi-Poisson describe it better. A plot like the one in §2 shows which shape the data support.
• Underdispersion (\(\hat\phi \lt 1\)) exists but is rare; a mixture never produces it, by (14). When you see it, think of regular processes (capped counts, events that inhibit each other) and of models such as Conway–Maxwell–Poisson or the binomial.

§15Summary

1. The Poisson GLM imposes \(\operatorname{Var}(Y_i) = \mu_i\). Real data usually have a larger variance; the diagnostic is \(\hat\phi = X^2/(n-p) \gt 1\).
2. The most common cause is unobserved heterogeneity: units with the same covariates have different rates. We model it in two layers: \(Z \sim \text{Gamma}(\alpha, \beta)\) and \(Y \mid Z \sim \text{Poisson}(Z)\).
3. The iterated expectation gives \(\mathbb E[Y] = \mathbb E[\mathbb E[Y \mid Z]] = \mathbb E[Z]\). The law of total variance gives \(\operatorname{Var}(Y) = \mathbb E[Z] + \operatorname{Var}(Z)\). So \(\operatorname{Var}(Y) - \mathbb E[Y] = \operatorname{Var}(Z) \gt 0\): the excess variance is the variance of the rates.
4. Integrating the Poisson against the Gamma, the marginal is the negative binomial \(\text{NB}(r = \alpha,\; p = \beta/(1+\beta))\). The Gamma is chosen because it is conjugate: the integral has a closed form.
5. For the GLM, \(Y_i \mid \varepsilon_i \sim \text{Poisson}(\mu_i\varepsilon_i)\) with \(\varepsilon_i \sim \text{Gamma}(\theta, \theta)\), giving \(\mathbb E[Y_i] = \mu_i\) and \(\operatorname{Var}(Y_i) = \mu_i + \mu_i^2/\theta\). As \(\theta \to \infty\) we are back to the Poisson.
6. With \(\theta\) fixed the NB is an exponential family, so it is a GLM. It uses the log link, its coefficients read like Poisson ones, and \(\theta\) is estimated by maximum likelihood.
7. Ignoring overdispersion does not bias the estimates but shrinks the standard errors: in the simulation the Poisson “95%” interval covered 48% of the time at \(\theta = 0.5\). The NB stayed close to 95% in every scenario and cost nothing when the data were Poisson.

§16Further reading

• M. Greenwood and G. U. Yule (1920). An inquiry into the nature of frequency distributions representative of multiple happenings… Journal of the Royal Statistical Society, 83(2). The origin of the Poisson–Gamma mixture, with accident data.
• P. McCullagh and J. A. Nelder (1989). Generalized Linear Models, 2nd ed., ch. 6. Chapman & Hall. Quasi-likelihood and overdispersion.
• J. F. Lawless (1987). Negative binomial and mixed Poisson regression. Canadian Journal of Statistics, 15(3).
• A. C. Cameron and P. K. Trivedi (2013). Regression Analysis of Count Data, 2nd ed. Cambridge. The complete reference on count models.
• J. M. Hilbe (2011). Negative Binomial Regression, 2nd ed. Cambridge. Practical, with many examples and the NB1, NB2 and other variants.
• A. Agresti (2015). Foundations of Linear and Generalized Linear Models, ch. 7. Wiley.
• W. N. Venables and B. D. Ripley (2002). Modern Applied Statistics with S, 4th ed., section 7.4. Springer. The basis of glm.nb.