10.3 Other count regressions

Real count data often show overdispersion, zero-inflation, or other features that the simple Poisson model cannot capture. Alternative models add flexibility.

10.3.1 Negative binomial (Gamma-Poisson)

The negative binomial (also called Gamma-Poisson) adds a dispersion parameter, allowing the variance to exceed the mean:

$\text{Var}(y_i) = \lambda_i + \frac{\lambda_i^2}{\phi}$

When $\phi \to \infty$ , the model reduces to Poisson. Smaller $\phi$ means more overdispersion.

# Negative binomial regression for overdispersed counts
set.seed(42)
n <- 100
x <- runif(n, 0, 10)

# Overdispersed data: lambda varies with a Gamma-distributed multiplier

lambda*base <- exp(0.5 + 0.2 * x)
y <- rpois(n, lambda*base * rgamma(n, shape = 3, rate = 3))

# Fit Poisson (wrong)

m_pois <- glm(y ~ x, family = poisson)

# Fit negative binomial (correct)

library(MASS)
m_nb <- glm.nb(y ~ x)

cat("Poisson vs Negative Binomial:\n")
cat(" Poisson: intercept =", round(coef(m_pois)[1], 3),
" slope =", round(coef(m_pois)[2], 3), "\n")
cat(" NegBin: intercept =", round(coef(m_nb)[1], 3),
" slope =", round(coef(m_nb)[2], 3), "\n")
cat(" NegBin dispersion (theta):", round(m_nb$theta, 2), "\n")

cat("\n AIC comparison:\n")
cat(" Poisson:", round(AIC(m_pois), 1), "\n")
cat(" NegBin: ", round(AIC(m_nb), 1), "\n")

if (AIC(m_nb) < AIC(m_pois)) {
cat("\n NegBin preferred (lower AIC).\n")
}

Click Run to execute this code in your browser with WebR

10.3.2 Zero-inflated Poisson

When there are more zeros than expected, a zero-inflated Poisson (ZIP) model uses a mixture:

With probability $\pi$ , $y_i = 0$ (structural zero).
With probability $1 - \pi$ , $y_i \sim \text{Poisson}(\lambda_i)$ .

$P(y_i = 0) = \pi + (1 - \pi) e^{-\lambda_i}$ $P(y_i = k) = (1 - \pi) \frac{\lambda_i^k e^{-\lambda_i}}{k!}, \quad k > 0$

# Zero-inflated Poisson: simulate and visualize
set.seed(42)
n <- 200
x <- runif(n, 0, 10)

# Structural zeros: 30% of observations always produce 0

pi_struct <- 0.3
is_structural <- rbinom(n, 1, pi_struct)

# Poisson counts for non-structural

lambda <- exp(0.3 + 0.15 * x)
y_pois <- rpois(n, lambda)

# Final: structural zeros override Poisson

y <- ifelse(is_structural, 0, y_pois)

# Count zeros

cat("Zero-inflated data:\n")
cat(" Total zeros:", sum(y == 0), "out of", n, "\n")
cat(" Proportion of zeros:", round(mean(y == 0), 3), "\n")
cat(" Expected zeros from Poisson alone:", round(sum(exp(-lambda)), 1), "\n")
cat(" Extra zeros (structural):", sum(y == 0) - round(sum(exp(-lambda)), 0), "\n")

# Plot

par(mfrow = c(1, 2))
hist(y, breaks = max(y) + 1, col = "steelblue",
main = "Zero-inflated counts", xlab = "y")
plot(x, y, pch = 16, col = adjustcolor("steelblue", 0.5),
main = "Zero-inflated: y vs x", xlab = "x", ylab = "y")
par(mfrow = c(1, 1))

Click Run to execute this code in your browser with WebR

10.3.3 Hurdle models

Hurdle models separate the process into two parts:

Hurdle: Is $y_i > 0$ or $y_i = 0$ ? (Binary: logistic regression)
Count: Given $y_i > 0$ , what is the count? (Truncated Poisson)

This is different from ZIP: in a hurdle model, all zeros come from the hurdle. In ZIP, zeros can come from either the structural process or the Poisson.

# Hurdle model concept: two-part model
set.seed(42)
n <- 200
x <- runif(n, 0, 10)

# Part 1: probability of non-zero (logistic)

p_nonzero <- 1 / (1 + exp(-(-1 + 0.3 * x)))
is_nonzero <- rbinom(n, 1, p_nonzero)

# Part 2: count given non-zero (truncated Poisson)

lambda <- exp(0.5 + 0.1 * x)
y_count <- rpois(n, lambda)

# Truncate at 1 (no zeros in the count part)

y_count <- y_count + 1

# Combine

y <- ifelse(is_nonzero, y_count, 0)

cat("Hurdle model data:\n")
cat(" Zeros:", sum(y == 0), "\n")
cat(" Non-zeros:", sum(y > 0), "\n")
cat(" Mean of non-zeros:", round(mean(y[y > 0]), 2), "\n")

# Fit two-part model

# Part 1: logistic regression for y > 0

m_hurdle <- glm(I(y > 0) ~ x, family = binomial)

# Part 2: Poisson for y given y > 0

m_count <- glm(y[y > 0] ~ x[y > 0], family = poisson)

cat("\nHurdle model:\n")
cat(" Part 1 (y > 0): intercept =", round(coef(m_hurdle)[1], 3),
" slope =", round(coef(m_hurdle)[2], 3), "\n")
cat(" Part 2 (count): intercept =", round(coef(m_count)[1], 3),
" slope =", round(coef(m_count)[2], 3), "\n")

Click Run to execute this code in your browser with WebR

Rethinking

Choose the model based on the data-generating process. The Poisson, negative binomial, ZIP, and hurdle models are all tools for count data, but they answer different questions. Overdispersion? Use negative binomial. Too many zeros from a separate process? Use ZIP. Too many zeros because some observations never start? Use a hurdle. The choice should be guided by theory, not just by which model fits best. Posterior predictive checks help verify that the chosen model reproduces the key features of the data.

Summary of Section 10.3

Negative binomial (Gamma-Poisson): adds a dispersion parameter for overdispersed counts.
Zero-inflated Poisson (ZIP): mixture of structural zeros and Poisson counts.
Hurdle model: separates the decision to count from the count itself.
Choose based on the data-generating process, not just fit.
Always use posterior predictive checks to verify the model.

Poisson regression

Summary