Poisson Regression for Count Data: Rates, Overdispersion, and Offsets

Emergency admissions per hospital-month, mutations per gene, insects per trap, claims per policyholder. These outcomes are counts: whole numbers, never negative, and often with a variance that grows alongside the mean. Fitting an ordinary linear regression to them is not merely inelegant — it produces predictions of minus two admissions and confidence intervals that cross zero.

Poisson regression is the standard answer, and it is easier than its reputation suggests. The model is still a linear model; the only change is that it predicts the logarithm of the expected count rather than the count itself, which keeps every fitted value positive and gives coefficients a rate-ratio interpretation.

Why a straight line fails on counts

Three properties of counts break the linear model’s assumptions, and it is worth naming them because each one has a visible symptom.

Property of counts	What the linear model assumes	Visible symptom
Cannot be negative	Errors are symmetric, unbounded	Negative fitted counts and interval bounds
Variance grows with the mean	Constant variance (homoscedasticity)	Funnel-shaped residuals
Many small values, few large ones	Approximately normal residuals	Skewed, non-normal residual distribution

You can sometimes rescue a linear model by transforming the outcome, but a log of a zero count is undefined, and adding a constant to every value makes the coefficient depend on the constant you chose. A model built for counts avoids the improvisation.

The Poisson model in one sentence

Poisson regression assumes each observed count comes from a Poisson distribution whose mean depends on the predictors through a log link: log(E[y]) = b0 + b1x1 + …. Exponentiating both sides gives E[y] = exp(b0) × exp(b1x1) × …, an expected count that is positive for any value of the predictors.

A defining feature of the Poisson distribution is that its variance equals its mean. That single constraint is the source of both the model’s efficiency and its most common failure, which is covered below.

# R
fit <- glm(admissions ~ age_group + season + beds,
           family = poisson(link = "log"), data = ward)
summary(fit)
exp(coef(fit))            # rate ratios

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

fit = smf.glm("admissions ~ age_group + season + beds",
              data=ward, family=sm.families.Poisson()).fit()
print(fit.summary())
print(fit.params.apply(np.exp))      # rate ratios

Reading the coefficients

Coefficients are on the log scale, so they are not directly interpretable as differences in counts. Exponentiate them and they become rate ratios: the multiplicative change in the expected count for a one-unit increase in the predictor, holding the others constant.

Coefficient (log scale)	Rate ratio exp(b)	Reading
0.00	1.00	No difference in rate
0.18	1.20	20% higher expected count per unit
-0.22	0.80	20% lower expected count per unit
0.69	2.00	Double the expected count per unit

Report the rate ratio with its confidence interval, not just the p-value. A rate ratio of 1.20 with a wide interval and one of 1.20 with a narrow one are different findings, and a reader cannot tell them apart from a significance star.

exp(cbind(RR = coef(fit), confint(fit)))

For a categorical predictor, the reference level carries the whole interpretation. A coefficient of 0.4 for treatment B is meaningless until the reader knows what treatment A did, so state the reference explicitly in the table.

Offsets: modelling rates rather than counts

Hospitals with more beds record more admissions, wards observed for longer record more events, and larger regions report more cases. When exposure differs across units, comparing raw counts is misleading; you want a rate per unit of exposure.

An offset handles this by entering the exposure on the log scale with a fixed coefficient of one, so the model predicts log(events) = log(exposure) + Xb, which is equivalent to modelling events per unit of exposure.

# R: admissions per 1000 bed-days
fit <- glm(admissions ~ age_group + season + offset(log(bed_days)),
           family = poisson, data = ward)

# Python
fit = smf.glm("admissions ~ age_group + season",
              data=ward, family=sm.families.Poisson(),
              offset=np.log(ward["bed_days"])).fit()

The distinction between an offset and a covariate matters. If the effect of exposure is not exactly proportional — if the rate genuinely rises with hospital size — then exposure belongs as a predictor with an estimated coefficient, not as an offset.

Overdispersion: the failure that hides

Real counts are usually more variable than a Poisson distribution allows. Patients cluster within hospitals, insects cluster within traps, and unobserved heterogeneity inflates the variance. This is overdispersion, and ignoring it does not bias the coefficients — it makes the standard errors too small, so p-values are too optimistic and intervals too narrow.

The quickest check is the dispersion statistic: the Pearson chi-square divided by the residual degrees of freedom. Values near 1 are consistent with Poisson; values well above 1.5 indicate overdispersion worth addressing.

# R
sum(residuals(fit, type = "pearson")^2) / df.residual(fit)

library(performance)
check_overdispersion(fit)

# Python
pearson_chi2 = np.sum(fit.resid_pearson ** 2)
print(pearson_chi2 / fit.df_resid)

Two established fixes exist, and the choice between them depends on what you believe about the data generating process.

Approach	What it changes	When to use
Quasi-Poisson	Standard errors only; coefficients stay Poisson	You want the same mean model with honest uncertainty
Negative binomial	A variance function with an extra parameter (NB2)	Variance grows faster than the mean; you want a likelihood
Zero-inflated / hurdle	Two processes: structural zeros and counts	More zeros than any count model predicts
Poisson with clustered SEs	Accounts for grouping, not extra variance	Repeated measures within units

# R
quasi <- glm(admissions ~ age_group + season, family = quasipoisson, data = ward)
nb    <- MASS::glm.nb(admissions ~ age_group + season, data = ward)
AIC(nb)

# Python
nb = smf.negativebinomial("admissions ~ age_group + season", data=ward).fit()
print(nb.summary())

Negative binomial is the more common choice in published work because it is a genuine likelihood, so information criteria and likelihood ratio tests are available. Quasi-Poisson is the conservative option when you only want to widen the standard errors without committing to a distributional story.

Too many zeros

A count distribution with a spike at zero is often a sign of two different processes: some units cannot have the event at all, and among those that can, the count follows a Poisson or negative binomial. Zero-inflated models mix a structural-zero component with a count component; hurdle models instead fit a binary model for any event and a truncated count model for the positive counts.

The two are not interchangeable. Zero-inflated models allow zeros to arise from either component, which suits processes where a susceptible unit happens to observe nothing. Hurdle models assume every zero comes from the same barrier, which suits processes where the event is either possible or not.

library(pscl)
zinb <- zeroinfl(visits ~ age + insurance | income, data = clinic, dist = "negbin")
summary(zinb)

hurdle_fit <- hurdle(visits ~ age + insurance | income, data = clinic, dist = "poisson")

Before reaching for either, check whether the excess zeros are real. A Poisson model with an offset for a small exposure can produce apparent zero inflation simply because many units had little chance of an event.

Checking the fit

Residual diagnostics for count models are less familiar than for linear models, and the standard residual plots are hard to read because residuals are discrete. Simulation-based tools are more informative than eyeballing.

library(DHARMa)
sim <- simulateResiduals(fit)
plot(sim)
testDispersion(sim)

# rootogram: observed vs fitted counts
library(countreg)
rootogram(fit)

A rootogram compares observed and fitted frequencies across the range of counts, hanging the differences from a reference line. Systematic departure at zero shows the model under-predicts zeros; a poorly fitting upper tail shows it cannot represent the large counts.

What to report

State the model family and link, the outcome definition and its observation window, the offset if one was used, the rate ratios with confidence intervals, the dispersion statistic before and after any correction, and how overdispersion or excess zeros were handled. Naming the family matters more than it sounds: a reader who assumes OLS will misread every coefficient.

Count outcomes also appear in less obvious places — variant calls per sample, reads per gene, adverse events per patient — and the same reasoning applies there. If you want the wider modelling context, Ampersand Academy covers generalised linear models from the ground up.

Practice

1. Fit a Poisson model to a count outcome and report the dispersion statistic before interpreting any coefficient.
2. Convert three coefficients into rate ratios and write one sentence interpreting each.
3. Refit with an offset for exposure and explain how the interpretation of the coefficients changed.
4. Compare Poisson, quasi-Poisson and negative binomial fits and state which you would defend and why.

Related tutorials

• Logistic Regression Explained: Odds Ratios, Probabilities, and What Coefficients Mean — the same generalised-linear-model machinery for a binary outcome
• Distributions in Statistics: Types with Examples — the Poisson distribution these models are built on
• Checking Statistical Assumptions in Practice: What to Plot, What to Trust, What to Do Next — diagnostics beyond overdispersion, applied to the same model