Count data is everywhere in the real world:

• Number of hospital admissions per day
• Customer complaints per month
• Goals scored per soccer match
• Traffic accidents per intersection
• Website clicks per user session
• Species observed per quadrat
• Patents filed per company per year

These outcomes share common characteristics: they are non-negative integers, often right-skewed, and their variance often depends on their mean. Linear regression is inappropriate for count data because it can predict negative counts, assumes constant variance, and ignores the discrete nature of the response.

Poisson regression is the natural GLM for count data. It assumes the response follows a Poisson distribution, uses a log link to ensure positive predictions, and automatically accounts for mean-dependent variance. Mastering Poisson regression opens the door to analyzing a vast array of real-world counting phenomena.

Before diving into regression, let's thoroughly understand the Poisson distribution itself.

Definition and Properties

A random variable $Y$ follows a Poisson distribution with parameter $\lambda > 0$ if:

$$P(Y = k) = \frac{\lambda^k e^{-\lambda}}{k!}, \quad k = 0, 1, 2, \ldots$$

Key properties:

• $E[Y] = \lambda$ (the parameter is the mean)
• $\text{Var}(Y) = \lambda$ (variance equals mean—equidispersion)
• $\text{Skewness} = 1/\sqrt{\lambda}$ (right-skewed, especially for small $\lambda$)
• As $\lambda \to \infty$, the Poisson approaches a Normal distribution

When Does the Poisson Arise?

The Poisson distribution emerges from the Poisson process, which models random events occurring in time or space under these conditions:

1. Independence: Events in non-overlapping intervals are independent
2. Homogeneity: The rate of events is constant over time/space
3. No simultaneous events: Events don't occur at exactly the same instant

If events occur at rate $r$ per unit time, then the count in an interval of length $t$ is $Y \sim \text{Poisson}(rt)$.

Poisson regression is a GLM that models count responses using the Poisson distribution with a log link.

Model Specification

Random Component: $$Y_i \sim \text{Poisson}(\lambda_i)$$

Each observation has its own mean $\lambda_i$ that depends on the predictors.

Systematic Component: $$\eta_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \cdots + \beta_p X_{ip}$$

Link Function (Log): $$\log(\lambda_i) = \eta_i = \mathbf{x}_i^\top \boldsymbol{\beta}$$

Equivalently: $$\lambda_i = \exp(\mathbf{x}i^\top \boldsymbol{\beta}) = e^{\beta_0} \cdot e^{\beta_1 X{i1}} \cdot \ldots \cdot e^{\beta_p X_{ip}}$$

Why the Log Link?

1. Positivity: $\lambda = e^\eta > 0$ for any $\eta$, ensuring valid predicted counts
2. Multiplicative effects: Predictors act multiplicatively on the mean count
3. Canonical link: Log is the canonical link for Poisson, giving nice mathematical properties
4. Rate interpretation: Log-linear models naturally express proportional changes

The Log-Likelihood

For independent Poisson observations $Y_1, \ldots, Y_n$:

$$\ell(\boldsymbol{\beta}) = \sum_{i=1}^n \left[ Y_i \log(\lambda_i) - \lambda_i - \log(Y_i!) \right]$$

Substituting $\log(\lambda_i) = \mathbf{x}_i^\top \boldsymbol{\beta}$:

$$\ell(\boldsymbol{\beta}) = \sum_{i=1}^n \left[ Y_i (\mathbf{x}_i^\top \boldsymbol{\beta}) - e^{\mathbf{x}_i^\top \boldsymbol{\beta}} \right] + \text{const}$$

Score Equations

Taking the derivative with respect to $\beta_j$:

$$\frac{\partial \ell}{\partial \beta_j} = \sum_{i=1}^n (Y_i - \lambda_i) X_{ij} = 0$$

In matrix form: $$\mathbf{X}^\top (\mathbf{Y} - \boldsymbol{\lambda}) = \mathbf{0}$$

This says: the residuals $(Y_i - \hat{\lambda}_i)$ must be orthogonal to each predictor—the same condition as in OLS, but with different residuals!

No Closed-Form Solution

Unlike linear regression, these equations are nonlinear in $\boldsymbol{\beta}$ (because $\lambda_i = e^{\mathbf{x}_i^\top \boldsymbol{\beta}}$). We must use iterative methods—typically Newton-Raphson or equivalently Iteratively Reweighted Least Squares (IRLS).

The log link means that Poisson coefficients have multiplicative interpretations—they describe proportional changes in the expected count.

The Incidence Rate Ratio (IRR)

Consider a single predictor: $$\log(\lambda) = \beta_0 + \beta_1 X$$

When $X$ increases by 1 unit:

• Log-mean increases by $\beta_1$
• Mean is multiplied by $e^{\beta_1}$

The quantity $e^{\beta_1}$ is called the Incidence Rate Ratio (IRR) or Rate Ratio (RR):

$$\text{IRR} = \frac{\lambda(X+1)}{\lambda(X)} = \frac{e^{\beta_0 + \beta_1(X+1)}}{e^{\beta_0 + \beta_1 X}} = e^{\beta_1}$$

Interpretation Examples

Suppose we model doctor visits per year: $$\log(\text{visits}) = 0.5 + 0.3 \cdot \text{chronic_conditions} - 0.1 \cdot \text{income_10k}$$

Confidence Intervals for IRRs

To get a confidence interval for the IRR:

1. Compute CI for $\beta$: $(\hat{\beta}j \pm z{\alpha/2} \cdot \text{SE}(\hat{\beta}_j))$
2. Exponentiate the endpoints: $(e^{\hat{\beta}j - z{\alpha/2} \cdot \text{SE}}, e^{\hat{\beta}j + z{\alpha/2} \cdot \text{SE}})$

The CI for IRR is asymmetric because of the exponential transformation. An IRR of 1.0 means no effect; CIs not containing 1.0 indicate statistical significance.

Percentage Change Interpretation

For small effects, the percentage change is approximately $100 \cdot \beta_j$%:

• $\beta = 0.10 \Rightarrow$ IRR = 1.105, or about 10% increase
• $\beta = -0.05 \Rightarrow$ IRR = 0.951, or about 5% decrease

This approximation breaks down for large $|\beta|$.

Often we want to model rates rather than raw counts. For example:

• Accidents per vehicle-mile driven
• Infections per person-year at risk
• Defects per 1000 units produced
• Deaths per 100,000 population

The key is that we have different exposures (time at risk, population size, etc.) for different observations.

The Offset Formulation

If $Y_i$ is the count and $t_i$ is the exposure (time, population, etc.), we model the rate $r_i = Y_i / t_i$ by:

$$Y_i \sim \text{Poisson}(\lambda_i = t_i \cdot r_i)$$

where $r_i$ is the rate per unit exposure. Taking logs:

$$\log(\lambda_i) = \log(t_i) + \log(r_i) = \log(t_i) + \mathbf{x}_i^\top \boldsymbol{\beta}$$

The term $\log(t_i)$ is called the offset—a predictor with coefficient fixed at 1.

$$\log(\lambda_i) = \text{offset}i + \beta_0 + \beta_1 X{i1} + \cdots$$

Interpretation with Offset

With an offset, coefficients represent effects on the rate (count per unit exposure):

• $e^{\beta_j}$ is the rate ratio, the multiplicative change in rate per unit change in $X_j$
• The baseline rate (when all predictors = 0) is $e^{\beta_0}$ per unit exposure

Assessing whether a Poisson model fits your data is crucial. Poor fit can lead to invalid inference. Here are key diagnostic tools.

Deviance and Pearson Residuals

Deviance Residuals: $$d_i = \text{sign}(Y_i - \hat{\lambda}_i) \sqrt{2\left[Y_i \log\left(\frac{Y_i}{\hat{\lambda}_i}\right) - (Y_i - \hat{\lambda}_i)\right]}$$

When $Y_i = 0$, $d_i = -\sqrt{2\hat{\lambda}_i}$.

Pearson Residuals: $$r_i^P = \frac{Y_i - \hat{\lambda}_i}{\sqrt{\hat{\lambda}_i}}$$

Under correct model specification, both should be approximately standard normal for large samples.

Residual Plots

1. Residuals vs. Fitted: Plot deviance or Pearson residuals against $\hat{\lambda}_i$. Look for patterns suggesting:

   • Heteroscedasticity (fan shape)
   • Nonlinearity (curved patterns)
   • Outliers (extreme residuals)

2. Residuals vs. Predictors: Check for nonlinear relationships not captured by the model

3. Q-Q Plot: Compare residual distribution to standard normal

Checking for Overdispersion

The most common problem in Poisson regression is overdispersion—variance exceeding the mean.

Detection Methods:

1. Dispersion Statistic: $$\hat{\phi} = \frac{1}{n-p} \sum_{i=1}^n \frac{(Y_i - \hat{\lambda}_i)^2}{\hat{\lambda}_i} = \frac{\text{Pearson } \chi^2}{n-p}$$

For a well-fitting Poisson model, $\hat{\phi} \approx 1$. Values significantly > 1 indicate overdispersion.

2. Deviance/df Ratio: $$\frac{D}{n-p} \approx 1$$

if the model fits. Much greater than 1 suggests overdispersion.

3. Cameron-Trivedi Test: Regress $(Y - \hat{\lambda})^2 - Y$ on $\hat{\lambda}$. A significant positive coefficient indicates overdispersion.

Consequences of Ignoring Overdispersion:

• Standard errors are too small
• Confidence intervals are too narrow
• p-values are too small (false positives)
• Model selection criteria (AIC) are invalid

When Poisson regression shows overdispersion, several approaches can address the problem.

Approach 1: Quasi-Poisson

The simplest fix is quasi-Poisson regression, which assumes: $$\text{Var}(Y_i) = \phi \cdot \lambda_i$$

where $\phi > 1$ is the dispersion parameter estimated from the data.

Effect:

• Point estimates remain the same as Poisson
• Standard errors are multiplied by $\sqrt{\hat{\phi}}$
• Confidence intervals widen; p-values increase
• No true likelihood, so AIC/BIC not directly applicable

Approach 2: Robust Standard Errors

Sandwich (robust) standard errors don't assume a specific variance structure. They estimate the variance of $\hat{\beta}$ directly from the data:

$$\hat{V}(\hat{\beta}) = (\mathbf{X}^\top \mathbf{W} \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{W} \hat{\mathbf{\Omega}} \mathbf{W} \mathbf{X} (\mathbf{X}^\top \mathbf{W} \mathbf{X})^{-1}$$

where $\hat{\mathbf{\Omega}}$ is diagonal with $(Y_i - \hat{\lambda}_i)^2$ on the diagonal.

Effect:

• Point estimates unchanged
• Standard errors robust to variance misspecification
• Valid even if overdispersion is not constant

Approach 3: Negative Binomial Regression

The negative binomial distribution explicitly models overdispersion through an additional parameter:

$$\text{Var}(Y_i) = \lambda_i + \frac{\lambda_i^2}{\theta}$$

where $\theta > 0$ is the overdispersion parameter. As $\theta \to \infty$, this reduces to Poisson.

Advantages over quasi-Poisson:

• A true likelihood exists → AIC/BIC applicable
• Explicitly models the variance structure
• Better for prediction

We'll cover negative binomial regression in detail on the next page.

Approach 4: Zero-Inflated Models

If excess zeros drive overdispersion, zero-inflated Poisson (ZIP) or zero-inflated negative binomial (ZINB) models are appropriate. These model:

1. A point mass at zero with probability $\pi$
2. A Poisson (or NB) distribution with probability $1-\pi$

$$P(Y = k) = \begin{cases} \pi + (1-\pi) e^{-\lambda} & k = 0 \ (1-\pi) \frac{\lambda^k e^{-\lambda}}{k!} & k > 0 \end{cases}$$

Let's work through a complete Poisson regression analysis, from data exploration to model diagnostics.

We've comprehensively covered Poisson regression—the canonical GLM for count data. Let's consolidate the key concepts:

What's Next:

In the next page, we'll explore negative binomial regression—the natural extension of Poisson regression that explicitly models overdispersion. We'll see how the negative binomial distribution arises as a Poisson-gamma mixture and when to prefer it over Poisson.