R - Poisson Regression

Modern Poisson Regression in R

Poisson Regression is a type of Generalized Linear Model (GLM) used when the response variable represents count data. This includes things like the number of customer complaints per day, the number of traffic accidents at an intersection, or the number of species found in a habitat. The counts must be non-negative integers (0, 1, 2, …).

The core assumption is that the response variable follows a Poisson distribution. The model predicts the logarithm of the expected count as a linear combination of predictor variables:

log(E[y]) = β₀ + β₁x₁ + β₂x₂ + … + βₙxₙ

y is the count response variable.
E[y] is the expected (mean) count.
β values are the model coefficients.
x values are the predictor variables.

The `glm()` Function

The primary function for fitting a Poisson regression in R is glm().

Syntax

glm(formula, data, family = poisson)

formula: Defines the relationship, e.g., counts ~ predictor1 + predictor2.
data: The data frame containing the variables.
family: Specifies the error distribution and link function. For Poisson regression, we use family = poisson, which implies a log link function.

Example: Analyzing Warp Breaks in Looms

We’ll use the built-in warpbreaks dataset. It contains data on the number of breaks in yarn during weaving, based on the wool type (A or B) and tension level (L, M, H). Our goal is to model how wool type and tension affect the number of breaks.

# Load modern libraries for data handling and model interpretation
# install.packages(c("tidyverse", "broom"))
library(tidyverse)
library(broom)

# Inspect the data
data(warpbreaks)
glimpse(warpbreaks)

The glimpse() output shows our variables:

Rows: 54
Columns: 3
$ breaks  <dbl> 26, 30, 54, 25, 70, 52, 51, 26, 67, 18, ...
$ wool    <fct> A, A, A, A, A, A, A, A, A, B, ...
$ tension <fct> L, L, L, L, L, L, M, M, M, L, ...

Fitting and Interpreting the Regression Model

# Fit the Poisson GLM. The formula models 'breaks' as a function of 'wool' and 'tension'.
poisson_model <- glm(breaks ~ wool + tension, data = warpbreaks, family = poisson)

# Get a comprehensive summary of the model
summary(poisson_model)

The summary output is key to understanding the model:

Call:
glm(formula = breaks ~ wool + tension, family = poisson, data = warpbreaks)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-3.6871  -1.6503  -0.4269   1.1902   4.2616

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  3.69196    0.04541  81.302  < 2e-16 ***
woolB       -0.20599    0.05157  -3.994 6.49e-05 ***
tensionM    -0.32132    0.06027  -5.332 9.73e-08 ***
tensionH    -0.51849    0.06396  -8.107 5.21e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 297.37  on 53  degrees of freedom
Residual deviance: 210.39  on 50  degrees of freedom
AIC: 493.06

Number of Fisher Scoring iterations: 4

Interpreting the Coefficients

In a Poisson model, the Estimate column shows the change in the log of the expected count. To make this intuitive, we must exponentiate the coefficient: exp(Estimate).

Baseline: The (Intercept) corresponds to the baseline group: wool ‘A’ and tension ‘L’ (the first level of each factor). The log-count of breaks for this group is 3.69.
woolB: The coefficient is -0.206. This means wool B has a log-count of breaks that is 0.206 lower than wool A, holding tension constant. The multiplicative effect is exp(-0.206) = 0.814. So, wool B is associated with about 18.6% fewer breaks than wool A.
tensionM & tensionH: Similarly, medium and high tension are associated with fewer breaks compared to low tension. For example, high tension (tensionH) is associated with exp(-0.518) = 0.596, or about 40.4% fewer breaks than low tension.
Significance: The low p-values (Pr(>|z|)) for all predictors suggest that both wool type and tension level have a statistically significant effect on the number of breaks.

Modern Interpretation with `broom`

Using broom::tidy() simplifies extracting and working with these coefficients.

# Tidy the model output
tidy_model <- tidy(poisson_model)

# Add the exponentiated coefficients for easier interpretation
tidy_model %>%
  mutate(rate_ratio = exp(estimate))

# Output:
# A tibble: 4 × 6
#  term        estimate std.error statistic      p.value rate_ratio
#  <chr>          <dbl>     <dbl>     <dbl>        <dbl>      <dbl>
#1 (Intercept)    3.69     0.0454     81.3  1.33e-164        40.1
#2 woolB         -0.206    0.0516     -3.99 6.49e-  5         0.814
#3 tensionM      -0.321    0.0603     -5.33 9.73e-  8         0.725
#4 tensionH      -0.518    0.0640     -8.11 5.21e- 16         0.595

Practical Considerations & Best Practices

Overdispersion: A key assumption of the Poisson model is that the mean of the count data is equal to its variance. In reality, the variance is often greater than the mean (a situation called overdispersion). You can check for this by seeing if the Residual Deviance is much larger than the Degrees of Freedom (in our case, 210.39 >> 50, suggesting overdispersion).
Alternative Models: If overdispersion is present, a Negative Binomial Regression is a better choice. You can fit it by using family = "negbin" from the MASS package (glm.nb() is a convenient wrapper).
Exposure: If counts are collected over different periods (e.g., accidents per month vs. per year), you should include an exposure or offset term in your model to normalize the rates. For example: glm(counts ~ predictor + offset(log(time)), family = poisson).