STATISTICS GRADUATE LEVEL • LINEAR MODELS AND GLMS

Poisson Regression

Modeling count data through the generalized linear model framework with a log link function.

Historical Context & Motivation

The need to model count data—discrete, non-negative integers representing the number of occurrences of an event in a fixed interval—has driven statistical innovation for over two centuries. Classical ordinary least squares regression assumes a continuous, normally distributed response variable, which makes it fundamentally unsuitable for counts that are bounded below by zero and often exhibit skewed distributions with variance that grows alongside the mean. Poisson regression arose as the principled solution to this modeling challenge, embedding the Poisson distribution within a regression framework that respects the inherent properties of count data.

The intellectual lineage of Poisson regression traces back to Siméon Denis Poisson's foundational 1837 treatise, but the full regression machinery required the development of generalized linear models (GLMs) in the 1970s. Before GLMs, analysts often resorted to transforming count data—applying square root or logarithmic transformations before fitting a linear model—an approach that introduced bias, distorted standard errors, and struggled with zero counts. The GLM framework provided a unified, elegant alternative that handled counts, binary outcomes, and continuous responses under one theoretical umbrella.

1837
Poisson Distribution Introduced
Siméon Denis Poisson publishes Recherches sur la probabilité des jugements, introducing the Poisson distribution to model rare events such as wrongful convictions in French courts.
1898
Ladislaus Bortkiewicz & Rare Events
Bortkiewicz famously applies the Poisson distribution to model deaths from horse kicks in the Prussian army, demonstrating its practical value for count data in applied settings.
1972
Nelder & Wedderburn Formalize GLMs
John Nelder and Robert Wedderburn publish their seminal paper unifying regression models for different exponential family distributions under the generalized linear model framework, providing the theoretical backbone for Poisson regression.
1983
McCullagh & Nelder's Canonical Text
The publication of Generalized Linear Models provides comprehensive treatment of Poisson regression including diagnostics, deviance analysis, and overdispersion, becoming the standard graduate reference.
1990s–Present
Extensions & Modern Applications
Zero-inflated, negative binomial, and hurdle models extend the Poisson framework. Software packages (R's glm(), Python's statsmodels) make Poisson regression accessible for epidemiology, ecology, insurance, and text mining.

The central question that Poisson regression addresses is deceptively straightforward: how does the expected count of events change as a function of one or more predictor variables? Whether the outcome is hospital readmissions per month, the number of species per plot, or insurance claims per policyholder, Poisson regression provides a model that naturally accommodates the discrete, non-negative character of counts while linking predictors to the response through a logarithmic function that ensures predicted means remain positive.

Core Principles & Definitions

Poisson regression is a member of the generalized linear model family, which means it is defined by three structural components: a random component specifying the distribution of the response, a systematic component assembling the linear predictor from covariates, and a link function connecting the two. Understanding these components and the assumptions they encode is essential before fitting any model to count data.

1

Random Component

The response variable Yi follows a Poisson distribution with mean μi. The Poisson is a member of the exponential family, with probability mass function P(Y = y) = e−μμy / y!.
2

Systematic Component

The linear predictor ηi = β₀ + β₁xi1 + ⋯ + βpxip combines covariates in a linear function of the parameters, just as in ordinary linear regression.
3

Log Link Function

The canonical link for the Poisson family is the natural logarithm: g(μi) = ln(μi) = ηi. This ensures that the fitted mean μi = exp(ηi) is always strictly positive.
4

Equidispersion Assumption

Under the Poisson model, the variance equals the mean: Var(Yi) = μi. This mean–variance relationship is the most restrictive and frequently violated assumption, motivating extensions such as quasi-Poisson and negative binomial models.
5

Independence of Observations

The counts Y₁, Y₂, …, Yn are assumed to be independent conditional on the covariates. Violation arises with clustered, longitudinal, or spatially correlated data, requiring GEE or mixed-effects Poisson models.
KEY TAKEAWAY
Think of the log link in Poisson regression as analogous to a volume knob on a stereo: the linear predictor (the knob's position) moves on an unrestricted real-number scale, but the exponentiation step converts it into a strictly positive output (the loudness). Just as you can turn the knob arbitrarily far in either direction yet never produce a negative volume, Poisson regression guarantees that predicted counts are always non-negative. The multiplicative interpretation of coefficients follows naturally: a one-unit increase in a covariate multiplies the expected count by eβ rather than adding β units.

Visual Explanation: The GLM Architecture

The following diagram illustrates the structural flow of a Poisson regression model, from raw covariates through the linear predictor, across the link function, to the predicted Poisson mean. This architecture is shared across all GLMs—only the distribution and link function change. Understanding this flow clarifies why the model's parameters live on the log scale while predictions live on the count scale.

The diagram traces the three components of a Poisson GLM: covariates enter the linear predictor (η = Xβ), which lives on the entire real line. The log link function maps η to a positive mean μ = exp(η), and the response Y is drawn from a Poisson distribution with that mean. The lower panel emphasizes the multiplicative coefficient interpretation.

Notice that the architecture enforces a separation of concerns. The linear predictor captures all covariate effects additively on the log scale, while the inverse link (exponentiation) maps these effects to the response scale. This means that the effect of each predictor is multiplicative on the count scale rather than additive: increasing xj by one unit multiplies the expected count by eβⱼ. This multiplicative structure naturally prevents negative predictions and captures the common empirical pattern where the variance of counts grows with their level.

Mathematical Framework

We now develop the mathematical foundations of Poisson regression, starting from the probability mass function, deriving the log-likelihood, and arriving at the score equations solved by iteratively reweighted least squares (IRLS). This derivation reveals why the Poisson distribution's membership in the exponential family guarantees a concave log-likelihood and thus a unique maximum likelihood estimate.

POISSON PROBABILITY MASS FUNCTION
P(Yᵢ = yᵢ | μᵢ) = exp(−μᵢ) × μᵢʸⁱ / yᵢ!
Where Yi is the observed count for observation i, and μi = E[Yi] is the mean parameter. Both E[Yi] and Var(Yi) equal μi.
LOG-LIKELIHOOD FUNCTION
ℓ(β) = Σᵢ [ yᵢ ln(μᵢ) − μᵢ − ln(yᵢ!) ] = Σᵢ [ yᵢ xᵢᵀβ − exp(xᵢᵀβ) − ln(yᵢ!) ]
Substituting μi = exp(xᵢᵀβ) from the log link. The last term −ln(yᵢ!) does not depend on β and is dropped during optimization. The log-likelihood is globally concave in β because −exp(xᵢᵀβ) is concave and yᵢxᵢᵀβ is linear.
SCORE EQUATIONS
∂ℓ/∂β = Xᵀ(y − μ) = 0
Here X is the n × (p+1) design matrix, y is the n × 1 vector of observed counts, and μ is the n × 1 vector of fitted means with μi = exp(xᵢᵀβ). The score equations are nonlinear in β and are solved iteratively via IRLS or Newton–Raphson.
IRLS UPDATE (FISHER SCORING)
β⁽ᵗ⁺¹⁾ = (XᵀWX)⁻¹ XᵀWz
Where W = diag(μi) is the weight matrix (since Var(Yi) = μi for Poisson), and z is the working response vector zi = ηi + (yi − μi)/μi. Each iteration resembles a weighted least squares regression on a linearized pseudo-response.
📐 Offset Terms for Rates
When observations have different exposure levels (e.g., person-years, area, time), include an offset term: ln(μᵢ) = xᵢᵀβ + ln(tᵢ), where tᵢ is the exposure. This is equivalent to modeling the rate μᵢ/tᵢ. In R: glm(y ~ x, family=poisson, offset=log(t)). The offset coefficient is fixed at 1 and is not estimated.

Model Diagnostics & Goodness of Fit

Fitting a Poisson regression model is only the first step; evaluating whether the model adequately describes the data is equally critical. The two primary diagnostics are the deviance and the Pearson chi-squared statistic, both of which compare the fitted model to a saturated model that perfectly reproduces the observed data. In well-specified models with moderate-to-large counts, both statistics follow an approximate χ² distribution with n − p − 1 degrees of freedom, providing a global goodness-of-fit test.

Left panel: under a well-specified Poisson model, Pearson residuals scatter uniformly in a band around zero with roughly constant spread across fitted values, yielding deviance/df ≈ 1. Right panel: overdispersion produces a fan-shaped residual pattern with increasing spread as fitted values grow, and deviance/df substantially exceeds 1.

Key Diagnostic Statistics

Common diagnostics for evaluating Poisson regression models
DiagnosticFormulaInterpretation
DevianceD = 2 Σᵢ [ yᵢ ln(yᵢ/μ̂ᵢ) − (yᵢ − μ̂ᵢ) ]D/df ≈ 1 indicates adequate fit. D/df ≫ 1 suggests overdispersion or model misspecification.
Pearson χ²χ² = Σᵢ (yᵢ − μ̂ᵢ)² / μ̂ᵢAlternative to deviance. Also approximately χ² distributed with n − p − 1 df under the null.
Dispersion Parameter (φ̂)φ̂ = χ² / (n − p − 1)φ̂ = 1 for true Poisson. φ̂ > 1 indicates overdispersion; φ̂ < 1 indicates underdispersion.
AICAIC = −2ℓ(β̂) + 2(p + 1)Penalized likelihood for model comparison. Lower AIC is preferred among competing models.
⚠️ Detecting Overdispersion in Practice
A quick rule of thumb: compute the ratio of the residual deviance to its degrees of freedom. If this ratio substantially exceeds 1 (say, greater than 1.5), consider fitting a quasi-Poisson model (which scales standard errors by √φ̂) or switching to a negative binomial regression (which adds an explicit overdispersion parameter). A formal test is the Cameron–Trivedi regression-based test for overdispersion.

Worked Example: Modeling Species Counts

An ecologist records the number of bird species observed across 200 forest plots. She hypothesizes that species count depends on plot area (in hectares) and whether the plot is adjacent to a water source (a binary indicator). We fit a Poisson regression and interpret the results.

Poisson Regression for Bird Species Counts
1
Step 1 — Specify the ModelLet Yi be the count of species in plot i. The model is: ln(μi) = β₀ + β₁(Areai) + β₂(Wateri), where Wateri = 1 if adjacent to water and 0 otherwise.
2
Step 2 — Fit Via Maximum LikelihoodUsing R: fit <- glm(species ~ area + water, family = poisson(link='log'), data = birds). The IRLS algorithm converges and returns the following estimates:
β̂₀ = 1.386, β̂₁ = 0.223, β̂₂ = 0.405 (all p < 0.01)
3
Step 3 — Interpret the Interceptβ̂₀ = 1.386, so for a plot with area = 0 hectares and no water adjacency, the expected species count is μ̂ = exp(1.386) ≈ 4.0. While area = 0 is extrapolation, this baseline anchors interpretation of the other coefficients.
Baseline expected count: exp(1.386) ≈ 4.0 species
4
Step 4 — Interpret the Area Coefficientβ̂₁ = 0.223 means that each additional hectare of area multiplies the expected species count by exp(0.223) ≈ 1.25. This is a 25% increase in the expected count per additional hectare, holding water adjacency constant. To construct a 95% confidence interval for the rate ratio: exp(0.223 ± 1.96 × SE(β̂₁)).
exp(β̂₁) = 1.25 → 25% increase per hectare
5
Step 5 — Interpret the Water Coefficientβ̂₂ = 0.405 gives exp(0.405) ≈ 1.50. Plots adjacent to water are estimated to have 50% more species on average than plots without water adjacency, controlling for area.
exp(β̂₂) = 1.50 → 50% increase for water-adjacent plots
6
Step 6 — Assess Goodness of FitThe residual deviance is 210.5 on 197 degrees of freedom, giving deviance/df = 210.5/197 ≈ 1.07. This ratio is close to 1, suggesting the Poisson model provides an adequate fit and there is no substantial evidence of overdispersion. The AIC is 842.3, which can be compared against alternative models.
Deviance/df ≈ 1.07 — no evidence of overdispersion
7
Step 7 — Predicted Count for a Specific PlotFor a plot with area = 5 hectares and water adjacency (Water = 1): ln(μ̂) = 1.386 + 0.223(5) + 0.405(1) = 1.386 + 1.115 + 0.405 = 2.906. Therefore μ̂ = exp(2.906) ≈ 18.3 species. The model predicts approximately 18 species for this particular combination of covariates.
Predicted count: exp(2.906) ≈ 18.3 species

Strengths, Limitations & Alternatives

Poisson regression occupies a central position in the analysis of count data, but like every model it rests on assumptions that may be violated in practice. Understanding its strengths relative to alternatives—and knowing when to move beyond the basic Poisson model—is a hallmark of rigorous applied statistical practice. The table below provides a structured comparison.

Comparison of common approaches for modeling count data
CriterionPoisson RegressionNegative Binomial RegressionOLS on Transformed Counts
Distributional assumptionVar(Y) = μ (equidispersion)Var(Y) = μ + αμ² (overdispersion modeled)Normal errors on log(Y + c) scale
Handles zerosYes, naturallyYes, naturallyRequires arbitrary constant c > 0
Coefficient interpretationexp(β) = rate ratio (multiplicative)Same as Poisson: exp(β) = rate ratioApproximate percentage change; back-transformation introduces bias
Robustness to overdispersionPoor — underestimates SE, inflates Type I errorGood — models extra-Poisson variationModerate — variance stabilizing, but model misspecified
EstimationMLE via IRLS; concave log-likelihoodMLE; requires iterative estimation of αOLS (closed form)
Best suited whenTrue equidispersion; or as a starting modelOverdispersed counts without excess zerosQuick exploratory analysis; counts are large
KEY TAKEAWAY
Think of Poisson regression as the default "first draft" for count data, much as simple linear regression is the first draft for continuous responses. Its equidispersion assumption is rarely perfectly satisfied in real data, but starting with Poisson regression establishes a principled baseline against which overdispersion diagnostics can be evaluated. If the deviance-to-df ratio is close to 1, you have a parsimonious, interpretable model. If it is much larger, you have clear evidence motivating a more flexible alternative like the negative binomial or zero-inflated Poisson model.

Connection to Advanced Theory & Extensions

Poisson regression serves as the gateway to a rich family of models for count and rate data. Understanding how it connects to more advanced methods—particularly within the exponential family and Bayesian frameworks—prepares you for the modeling challenges encountered in contemporary applied statistics and biostatistics.

Poisson regression and its advanced extensions
Basic Poisson RegressionAdvanced ExtensionKey Difference
Fixed dispersion (φ = 1)Quasi-PoissonEstimates dispersion φ̂ from data; adjusts SE by √φ̂ without specifying a full likelihood
Var(Y) = μNegative BinomialVar(Y) = μ + αμ²; adds a Gamma-distributed heterogeneity parameter to Poisson rates
Single source of zerosZero-Inflated Poisson (ZIP)Mixture model: zeros from structural process + Poisson counts (including random zeros)
Independent observationsPoisson GLMMAdds random effects for clustered/longitudinal data: ln(μᵢⱼ) = xᵢⱼᵀβ + uⱼ
Frequentist MLEBayesian Poisson RegressionPlaces priors on β; posterior estimated via MCMC. Natural for small samples and informative priors
Parametric link functionGAM with Poisson FamilyReplaces linear predictor with smooth functions: ln(μ) = β₀ + s(x₁) + s(x₂), capturing nonlinearity

A deep insight from exponential family theory connects Poisson regression to other GLMs: the canonical link for any exponential family distribution is the function that makes the sufficient statistic equal to the linear predictor, simplifying the score equations to Xᵀ(y − μ) = 0. For the Poisson, the canonical link is the natural log; for the binomial, it is the logit; for the Gaussian, it is the identity. This unification is the core achievement of the Nelder–Wedderburn GLM framework and provides a common template for estimation, inference, and diagnostics across all three families.

🔭 Looking Ahead
In modern data science, Poisson regression underpins models for web click counts, text mining (word frequencies), single-cell RNA-seq gene expression, and spatial point process intensities. The Poisson likelihood trick also enables fitting certain survival and logistic regression models through equivalences with Poisson regression on restructured data, a connection exploited in Whitehead (1980) and Holford (1980).

Practice Problems

PROBLEM 1CONCEPTUAL
Explain why applying ordinary linear regression to count data is problematic. Describe at least three specific violations or shortcomings that arise when the response is a non-negative integer, and explain how the Poisson regression framework addresses each one.
PROBLEM 2BASIC CALCULATION
A Poisson regression model for the number of emergency room visits per patient has the equation ln(μ) = 0.80 + 0.15 × Age_Group − 0.40 × Insured, where Age_Group is coded 1 for elderly and 0 for non-elderly, and Insured is coded 1 for insured and 0 for uninsured. Compute the expected number of ER visits for an elderly insured patient and interpret each coefficient as a rate ratio.
PROBLEM 3INTERMEDIATE
A researcher fits a Poisson regression to model the number of traffic accidents per intersection per year, using traffic volume (in thousands of cars/day) and number of lanes as predictors. The residual deviance is 385.2 on 147 degrees of freedom. Compute the dispersion estimate, assess whether overdispersion is present, and describe two approaches the researcher could take to address it.
PROBLEM 4APPLIED
An epidemiologist studies tuberculosis (TB) case counts across 50 districts over varying time periods (1 to 5 years of observation per district). She has data on population density and poverty rate as covariates. Explain why an offset term is necessary in this analysis, write out the full model specification including the offset, and interpret a hypothetical coefficient β̂₁ = 0.03 for population density (measured in people per km²).
PROBLEM 5CRITICAL THINKING
Prove that the log-likelihood function for Poisson regression is globally concave in the parameter vector β, and explain the practical consequence of this concavity for maximum likelihood estimation. Then discuss a scenario where, despite concavity, numerical issues can arise during IRLS fitting, and explain why.

Poisson Regression — Summary

Poisson regression is the canonical generalized linear model for count data. It assumes the response follows a Poisson distribution with mean μ linked to a linear predictor through the natural log link function, ensuring predicted counts are strictly positive. Parameters are estimated via maximum likelihood using the IRLS algorithm, and the globally concave log-likelihood guarantees a unique solution. Coefficients have an elegant multiplicative interpretation: exp(βⱼ) gives the rate ratio associated with a one-unit increase in the jth covariate.

The critical assumption of equidispersion (Var(Y) = μ) should be assessed via the deviance-to-degrees-of-freedom ratio and residual plots. When overdispersion is detected, alternatives include quasi-Poisson models, negative binomial regression, or zero-inflated models. For data with varying exposure, an offset term converts the model from counts to rates. Poisson regression is both a powerful stand-alone tool and the foundation for a broad family of extensions—mixed-effects, Bayesian, and semiparametric models—that address clustered data, prior information, and nonlinear covariate effects.

Varsity Tutors • Statistics Graduate Level • Poisson Regression