AdvancedIndependently reviewed, not yet spot-checked by a human

Poisson and negative binomial regression

Why counts cannot go into a linear model, why Poisson's equidispersion assumption almost never holds, how to detect overdispersion, what quasi-Poisson and the negative binomial each repair, and the offset — the difference between a count and a rate is the thing clinical research gets wrong most often.

This kind of outcome

The third large family of clinical outcomes is the count: seizures in a year, falls in a quarter, infections during an admission, emergency visits in a month. What they have in common is that they are non-negative integers, they have no upper bound, and they are usually right-skewed — most people sit at small numbers and a few sit very high.

Linear regression handles continuous outcomes and logistic regression handles binary ones; a count needs a third link function.

Counts have a second feature, and it matters more than the first: a count means almost nothing until a window of observation is attached to it. “Five seizures” over two weeks and “five seizures” over two years are entirely different clinical situations. Section six is devoted to that, and it is the most useful section on the page.

Why you cannot simply use linear regression

MASS::epil is Thall and Vail’s trial of the antiepileptic progabide: 59 patients (28 on placebo, 31 on progabide), with a 8-week baseline seizure count recorded before randomisation and four follow-up records afterwards, one every 2 weeks.

The main analysis on this page is at the patient level — one row per person, the outcome being the total count over the whole 8 weeks of follow-up — so that the correlation between one person’s four observations is not quietly ignored.

The shape of that outcome: median 16, maximum 302, mean 33.05, variance 2,078. Hold on to those last two numbers; section four uses them.

Feed this to least squares and it breaks in three places, none of them a theoretical worry:

Fitted values go negative. This model gives 5 patients a predicted seizure count below 0, the lowest being -5.75.
Prediction intervals are worse still. For a patient with a baseline count of 6 and an age of 25, OLS returns a point estimate of -7.19 seizures with a 95% prediction interval of -59.68 to 45.29. A negative count does not exist.
The variance is not constant. The residual standard deviation over the lower half of the fitted values is 10.71 and over the upper half 33.19, a three-fold difference. Every OLS standard error, p-value and confidence interval is built on the assumption that this does not happen.

Logging the outcome first does not rescue it either: the dataset contains patients whose count is zero across the entire follow-up — 1 of them — and $\log 0$ is undefined. A patch such as adding 0.5 makes the coefficients depend on how much you added.

Poisson regression

Poisson regression models the mean count on the log scale:

\log \mathbb{E}[Y \mid x] = \beta_0 + \beta_1 x_1 + \cdots + \beta_k x_k

Three immediate consequences: exponentiating a fitted value always returns something positive; $\exp(\beta)$ is a rate ratio, also called the incidence rate ratio (IRR), a multiplicative effect just like an OR or an HR — and note that the RR column in every table below is this rate ratio, not the risk ratio that the same abbreviation stands for on the logistic-regression page; and the model carries a built-in property that the variance grows with the mean.

Term	RR	95% CI	p
Progabide vs placebo	0.983	0.894–1.080	0.715
log(baseline count), per unit	3.405	3.195–3.629	< 0.001
log(age), per unit	1.800	1.451–2.233	< 0.001

Both continuous covariates were logged before entering the model, so “per unit” in the table means one unit on the log scale, not one more seizure or one more year: multiply the baseline count by e (roughly 2.72) and the follow-up rate is multiplied by 3.405; multiply age by e and the rate is multiplied by 1.800. Reading it as “each additional baseline seizure multiplies the rate by 3.405” is wrong.

It looks clean. But every confidence interval in that table is wrong, for the reason in the next section.

Equidispersion: Poisson’s only assumption, and almost never a true one

The Poisson distribution has a single parameter, so its variance equals its mean — this is called equidispersion. It is not an option you can switch off; it is a property of the distribution itself, and the model’s standard errors rest entirely on it.

Real clinical counts almost never behave that way, for an easy reason: Poisson describes events that are independent of one another occurring at a rate that is the same for everybody, and patients are not homogeneous. Some are inherently frequent seizers and some are not, and individual variation the model never accounted for pushes the overall variance above the mean. That is overdispersion.

In this dataset the overall variance is 2,078 against a mean of 33.05, a ratio of 62.9. That is a crude ratio with no covariates adjusted for; the formal check is the model’s Pearson chi-square divided by its degrees of freedom:

\hat{\phi} = \frac{1}{n - p}\sum_i \frac{(y_i - \hat{\mu}_i)^2}{\hat{\mu}_i}

Here that is 608.0 / 55 = 11.05 (p < 0.001). Under equidispersion the value should sit near 1.

Left: the two models' expected count distributions barely differ — the model contains log(baseline count), a very strong predictor, so Poisson also predicts patients with dozens of seizures. This panel is therefore not the evidence of overdispersion. Right: patients grouped into fifths by fitted value; each group's observed variance sits far above the equidispersion line (variance = mean). The evidence for overdispersion is this panel together with φ̂ = 11.05, not the left one.Plotting script figures/scripts/B2-03-poisson.R

Seizure count	Observed patients	Expected under Poisson	Expected under negative binomial
0-4	3	3.0	4.4
5-9	6	9.5	9.2
10-14	17	9.2	8.7
15-19	6	6.7	6.9
20-29	8	7.9	9.3
30-49	7	10.1	9.6
50-99	9	9.9	7.8
100+	3	2.7	3.0

Two repairs: quasi-Poisson and the negative binomial

All three coefficients from all three models side by side (same data, same formula, only the variance assumption changes):

Term	Model	Variance assumption	RR	95% CI	p	SE
Progabide vs placebo	Poisson	μ	0.983	0.894–1.080	0.715	0.048
Progabide vs placebo	quasi-Poisson	φ·μ	0.983	0.718–1.345	0.913	0.160
Progabide vs placebo	Negative binomial	μ + μ²/θ	0.769	0.574–1.031	0.079	0.149
log(baseline count), per unit	Poisson	μ	3.405	3.195–3.629	< 0.001	0.033
log(baseline count), per unit	quasi-Poisson	φ·μ	3.405	2.755–4.209	< 0.001	0.108
log(baseline count), per unit	Negative binomial	μ + μ²/θ	2.823	2.312–3.447	< 0.001	0.102
log(age), per unit	Poisson	μ	1.800	1.451–2.233	< 0.001	0.110
log(age), per unit	quasi-Poisson	φ·μ	1.800	0.879–3.684	0.114	0.365
log(age), per unit	Negative binomial	μ + μ²/θ	1.386	0.711–2.699	0.337	0.340

The conclusion about the treatment does not change; the conclusion about log(age) does. The three treatment rows reach no statistical significance in any of the three models (p 0.715 / 0.913 / 0.079). But under Poisson, log(age) comes out at RR 1.800 (1.451–2.233, p < 0.001) — apparently a strong predictor. Switch to quasi-Poisson and the point estimate does not move a digit while the interval widens until it crosses 1 (0.879–3.684, p 0.114); the negative binomial agrees (p 0.337). This is what the spurious significance of the previous callout looks like in practice: an effect the Poisson output reported at p < 0.001 falls back to no statistical significance once the variance assumption is right — the data did not change, only how large the standard error should have been. log(base), by contrast, survives all three models (p < 0.001 in every row); it is the genuinely strong predictor in this dataset.

quasi-Poisson does not change the distribution at all; it multiplies every standard error by $\sqrt{\hat{\phi}}$ . Not one point estimate moves (the first two rows for each term carry identical RRs); only the uncertainty changes. The estimated dispersion parameter is 11.05.

The negative binomial swaps in a different distribution: it layers a gamma distribution on top of the Poisson to describe how individuals differ, which makes the variance $\mu + \mu^2/\theta$ . The estimated $\theta$ here is 3.69 (standard error 0.80); as $\theta \to \infty$ it returns to Poisson. At a mean count of 33.05 it permits a variance of 329, where Poisson permits only 33.05.

The formal test: a likelihood ratio statistic of 377.8, p < 0.001. Model fit differs greatly too (AIC 848 against 472).

How to choose between them:

Whether the point estimates move. quasi-Poisson leaves them alone; the negative binomial does not (the treatment effect moves from 0.983 to 0.769 here), because it reweights the observations — high counts no longer carry as much weight. That gap is itself a signal of how severe the overdispersion is.
The negative binomial is a genuine probability model, so it has a likelihood, an AIC, likelihood ratio tests, and a predictive distribution you can simulate from. quasi-Poisson is not a full model (it has no likelihood), and its AIC is undefined.
quasi-Poisson is more robust: it assumes nothing about the shape of the overdispersion, only that the variance is proportional to the mean. When the true variance structure is not the negative binomial’s shape, quasi is less likely to mislead.
In practice, fit the negative binomial first, and switch to quasi-Poisson or robust (sandwich) standard errors if $\theta$ is estimated unstably or the residuals suggest the variance does not look like $\mu + \mu^2/\theta$ .

What not to do: seeing overdispersion and switching to linear regression, or dichotomising the count into “any versus none” and running logistic regression. The latter discards the difference between one seizure across the whole follow-up and dozens of them, which is usually the clinically important part. For residuals and influential points, see regression model diagnostics.

As for what the trial itself concluded: both models detected no statistically significant association between progabide and seizure count (negative binomial RR 0.769, 0.574–1.031, p 0.079). The lower bound of that interval still accommodates a fairly substantial reduction, so the correct way to write it is “this analysis did not reach statistical significance, and the interval is not narrow enough to rule out a moderate effect” — not “progabide does not work”.

The offset: a count and a rate are two different things

This is the most useful section on the page, and the place clinical research goes wrong most often.

The same count attached to different observation windows describes completely different clinical facts. Five seizures in a year and five seizures in a month are the same number and twelve times apart in severity. So what a model usually needs to estimate is a rate — events divided by person-time — rather than the count itself.

The way to do it is to add the observation time on the log scale as a term whose coefficient is fixed at 1. That term is the offset:

\log \frac{\mathbb{E}[Y]}{T} = \beta_0 + \beta_1 x_1 + \cdots \quad\Longleftrightarrow\quad \log \mathbb{E}[Y] = \log T + \beta_0 + \beta_1 x_1 + \cdots

offset(log(T)) is not the same as putting log(T) in as an ordinary covariate: the latter estimates a coefficient for it, which lets the data decide “how much time matters” — when by definition the relationship between time and count is one-to-one.

This dataset demonstrates the cost of ignoring the offset on its own. Each baseline record covers 8 weeks and each follow-up record 2 weeks. Stack the two kinds of record together and ask whether seizures became less frequent after randomisation:

Two bar charts. The left compares the mean seizure count per record, where the baseline period is clearly higher than the follow-up period. The right compares the seizure rate per patient per week, where the two are almost exactly equal. — Left: mean count per record, 31.2 at baseline against 8.3 during follow-up. Right: the same data converted to events per patient per week, 3.90 against 4.13. The entire difference comes from the length of the observation window.Plotting script `figures/scripts/B2-03-poisson.R`

Model	Follow-up vs baseline ratio	95% CI	Reading
No offset (comparing counts)	0.265	0.248–0.282	“Seizures fell by 74%“
With offset (comparing rates)	1.059	0.993–1.128	No change in the seizure rate was detected

The first row’s apparent 74% reduction is an artefact of window length, one hundred per cent of it, with no treatment effect in it at all. Each baseline record covers 8 weeks and each follow-up record 2 weeks, and the ratio 0.265 is almost exactly 2 divided by 8.

Its confidence interval is very narrow and its p-value tiny — “highly significant” statistically and entirely worthless clinically. Which is precisely why significance cannot be used to judge whether an analysis is right.

What this mistake looks like in real studies:

Follow-up time differs between groups (people withdraw early, people die), yet event counts are compared directly
People are enrolled in a national health insurance database for different numbers of months, yet numbers of visits are compared directly
Hospitals recruited over different periods, yet numbers of reports are compared directly
Wards have different numbers of beds, yet numbers of infections are compared directly (the denominator should be bed-days)

The rule of thumb is simple: ask what the outcome is “per”. Per person-year, per thousand bed-days, per hundred person-months — whatever answers “per what” belongs in the offset.

Run it yourself

library(MASS)
data(epil, package = "MASS")

# Patient level: one row per person, so the correlation among one person's
# four observations is not ignored
pt <- aggregate(y ~ subject + trt + base + age, data = epil, FUN = sum)

# Poisson
m_p <- glm(y ~ trt + log(base) + log(age), data = pt, family = poisson)
summary(m_p)

# ⚠️ Always compute the dispersion yourself; summary() will not give it to you
sum(residuals(m_p, type = "pearson")^2) / m_p$df.residual

# quasi-Poisson: point estimates unchanged, standard errors inflated
summary(glm(y ~ trt + log(base) + log(age), data = pt, family = quasipoisson))

# Negative binomial
m_nb <- glm.nb(y ~ trt + log(base) + log(age), data = pt)
summary(m_nb)                       # theta is printed at the end
2 * (logLik(m_nb) - logLik(m_p))    # LR statistic; halve the p-value

# offset: put the observation time in and the comparison becomes a rate, not a count
long <- rbind(
  data.frame(subject = pt$subject, phase = "Baseline",  count = pt$base, weeks = 8),
  data.frame(subject = epil$subject, phase = "Follow-up", count = epil$y, weeks = 2))
glm(count ~ phase, data = long, family = poisson)                      # ❌ compares counts
glm(count ~ phase + offset(log(weeks)), data = long, family = poisson) # ✅ compares rates

Verified with R 4.6.0 and MASS 7.3.65. glm.nb() lives in MASS, so there is nothing extra to install.

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

epil = sm.datasets.get_rdataset("epil", "MASS").data
pt = (epil.groupby(["subject", "trt", "base", "age"], as_index=False)["y"].sum())

m_p = smf.glm("y ~ trt + np.log(base) + np.log(age)", data=pt,
              family=sm.families.Poisson()).fit()
print(m_p.summary())
print(m_p.pearson_chi2 / m_p.df_resid)        # dispersion

m_nb = smf.glm("y ~ trt + np.log(base) + np.log(age)", data=pt,
               family=sm.families.NegativeBinomial(alpha=1 / 3.69)).fit()

# For the offset use exposure= (statsmodels logs it for you) or offset=np.log(weeks)
# smf.glm("count ~ phase", data=long, family=sm.families.Poisson(),
#         exposure=long["weeks"]).fit()

statsmodels estimates the negative binomial's alpha (= 1/theta) first and feeds it back in, whereas R's glm.nb() estimates it jointly, so the two thetas may differ slightly.

Common misuses

Misuse	Why it is wrong
Running counts through linear regression	Fitted values can be negative and the variance is not constant
Logging the outcome and running OLS	Undefined when a zero is present; adding a constant makes the coefficients depend on how much you added
Fitting a Poisson model and not checking the dispersion	Overdispersion underestimates the standard errors and manufactures significance
Comparing counts directly when observation time differs between units	Use an offset and compare rates, or what you measure is window length
Entering `log(T)` as an ordinary covariate instead of an offset	Time and count are one-to-one by definition; there is nothing to estimate
Dichotomising the count into “any versus none” and running logistic regression	Throws away the difference between one seizure and dozens
Drawing conclusions from a Poisson p-value when the dispersion is large	That p-value was computed under the wrong variance assumption
Comparing models using a quasi-Poisson AIC	quasi has no likelihood, so AIC is undefined
Reaching for a zero-inflated model as soon as there are many zeros	First establish whether the zeros are structural (some people simply cannot have the event) or merely overdispersion; the negative binomial handles the latter
Treating repeated measurements as independent rows	Between-subject standard errors come out too small; use GEE or a mixed effects model
Reading an IRR as though it were an OR or an HR	The denominators differ: IRR compares person-time, OR compares odds, HR compares instantaneous hazard
Writing “no effect” when the result was not significant	Write “no difference was detected” and give the interval

Reproducing every number on this page

/opt/homebrew/bin/Rscript figures/scripts/B2-03-poisson.R

Read the figure

The answer comes from the same statistical output that produced this page's figures, not from a number typed in beside them.

Which number decides whether these seizure-count data need overdispersion handled?

Show the answer and why

Correct answer: 11.1 - the dispersion left over after fitting the Poisson model

The figure to read is 11.1, the dispersion remaining after the fit (Pearson chi-square over degrees of freedom). The raw 62.9 is only the first signal: covariates are expected to absorb some of the variance, so the pre-fit ratio always runs high and judging by it overreacts. 33.1 is the mean count, which says nothing about dispersion - overdispersion is about the variance relative to the mean, not about the mean being large. Here the ratio is still far above one after fitting, so the standard errors cannot be reported on Poisson terms.

If ordinary linear regression is used for these counts, which number in the output most clearly shows the model form is wrong?

Show the answer and why

Correct answer: The smallest fitted value is -5.7 - the model predicts counts that cannot exist

-5.7 is the lowest of all the fitted values, and a seizure count cannot be negative - linear regression has no mechanism preventing it, so this is not imprecision but the wrong model form. -7.2 is just the worked example in the text, and one odd case really would prove little; the point here is that the whole lower end of the fitted values sits below zero. 0.0 is the smallest observed count: legitimate data do not make a model legitimate, and a prediction interval whose lower bound is negative is a prompt to think again about the outcome's distribution.

This Poisson model has three terms, and only one answers the trial's question. Which is the rate ratio for Progabide versus placebo?

Show the answer and why

Correct answer: 0.983 - the treatment term, whose confidence interval crosses one and does not reach statistical significance

0.983 is the treatment term, and its interval crosses one: these data did not detect a difference between Progabide and placebo, which is not the same as showing there is none. 3.405 belongs to log(baseline count) and 1.800 to log(age). Those much larger numbers read as "a big effect", but they are covariates and answer something other than the trial's question - and the one row that does answer it sits between the two conspicuous ones.

Chapters that use this method

Watch next

ENVery Normal· 12 minPuts linear, logistic and Poisson regression inside one framework. Watch it and the three chapters turn out to be the same thing with a different link function.

繁中MeDA School（洪弘）· 41 minA graduate-level lecture on generalised linear models, in Mandarin. No verified Chinese-language video covers count regression on its own; this is the closest.

Sources and licences

This page is original writing