Class 19 · Models for Overdispersed Counts

STAT 517: Advanced Statistical Models · Fall 2026

Published

October 27, 2026

This interactive companion follows Chapter 13 and the reviewed class19_overdispersed_counts.R. The Galapagos example checks a Poisson regression for overdispersion and then fits quasi-Poisson and negative binomial models; the school-absences example is the in-class lab.

Press Run Code on each box in order. Later boxes reuse earlier objects, and you can edit a box and run it again. Both data files are included, and everything uses base R except glm.nb() from MASS, which loads with the page. Questions 13.1–13.2 match the handout; Question 13.1 is the in-class lab, whose group assignments and submission template are posted where the instructor hosts the lab threads, and Question 13.2 is for after class.

Goals

  • Define overdispersion and explain why the Poisson mean–variance restriction is often too strict.
  • Detect overdispersion with the deviance test, the Cameron and Trivedi test, and \(\hat\phi\).
  • Fit a quasi-Poisson model and explain why its estimates equal the Poisson estimates while its standard errors are larger.
  • State the negative binomial distribution, its mean and variance, and the role of \(\theta\).
  • Compare count models with AIC or quasi-AIC on the right scale and check linearity.
  • Interpret incidence rate ratios, percent changes, and their confidence intervals.

Part 1 · Detecting overdispersion

1. Galapagos plant species

Botanists counted the unique plant species (Species) on 30 Galapagos Islands, with five geographic covariates: the island’s Area (km²), the Elevation of its highest point (m), the distance to the Nearest island (km), the distance from Santa Cruz island (Scruz, km), and the area of the adjacent island (Adjacent, km²).

The marginal mean and variance below pool all 30 islands and ignore the covariates. A variance far above the mean only suggests a problem: Poisson regression assumes that each conditional distribution \(y \mid \mathbf{x}\) is Poisson, and islands with different means already make the marginal variance exceed the marginal mean. So we fit the model and check it.

2. Fit the Poisson model

According to the Poisson model, all five covariates are significantly associated with the number of species. Can we trust that?

3. Is there overdispersion?

The deviance \(D\) on \(n-(p+1)\) degrees of freedom tests \(H_0\): the Poisson distribution fits the data well. The Cameron and Trivedi test tests \(H_0: \text{var}(y \mid \mathbf{x}) = \lambda\) against \(H_1: \text{var}(y \mid \mathbf{x}) = \lambda + c\,f(\lambda)\), where \(c > 0\) indicates overdispersion. In R it is AER::dispersiontest(), which by default uses \(f(\lambda) = \lambda\), tests \(c > 0\), and prints the estimated “dispersion” \(1 + c\), the factor by which the variance exceeds the mean:

AER::dispersiontest(gala.pois)   # in RStudio, after install.packages("AER")

The box below computes the same test directly, so the page does not need to load AER. Under \(H_1\), \(E[(y - \lambda)^2 - y] = c\lambda\), so the average of \([(y_i - \hat\lambda_i)^2 - y_i]/\hat\lambda_i\) estimates \(c\); its \(z\) statistic is compared with a standard normal. The output matches AER::dispersiontest() digit for digit.

Under a correct Poisson model the Pearson residuals have variance about 1, so most would fall within \(\pm 2\):

Part 2 · Quasi-Poisson regression

4. Fit the quasi-Poisson model

The quasi-Poisson model specifies only the conditional mean and variance, \(E(y_i \mid \mathbf{x}_i) = \lambda_i\) and \(\text{var}(y_i \mid \mathbf{x}_i) = \phi\lambda_i\), with the same log-linear mean as Poisson regression. It has no pmf and no likelihood. Its estimates solve the Poisson estimating equations, so they equal the Poisson estimates; every standard error is multiplied by \(\sqrt{\hat\phi}\), where \(\hat\phi = X^2/[n-(p+1)]\) is the Pearson statistic divided by the residual degrees of freedom.

5. Interpret the quasi-Poisson rate ratios

After controlling for the other covariates, we expect a one km² increase in island area to be associated with a 0.058% decrease in the mean number of plant species; we are 95% confident that the decrease is between 0.029% and 0.087%. Similarly, a one meter increase in elevation is associated with a 0.35% increase (95% CI: 0.26% to 0.45%), and a one km² increase in the area of the adjacent island with a 0.066% decrease (95% CI: 0.036% to 0.102%).

6. Quasi-AIC

R does not give an AIC for quasi-Poisson models. The quasi-AIC divides the Poisson log-likelihood at the fitted means by \(\hat\phi\) and adds the usual penalty:

\[ \text{QAIC} = -\frac{2\,\ell(\hat{\boldsymbol\beta})}{\hat\phi} + 2K, \qquad K = p + 2, \]

counting the \(p+1\) coefficients, including the intercept, and one more for \(\hat\phi\). To compare several quasi-Poisson models, use one \(\hat\phi\) for all of them, from the largest model, and give each model its own \(K\). A quasi-AIC is on a different scale from an ordinary AIC: never compare it with the Poisson or negative binomial AIC.

In RStudio, MuMIn::QAIC(gala.pois, chat = phi) gives the same value. Pass the Poisson fit: a quasi-Poisson fit has no log-likelihood, so MuMIn returns NA for it.

Part 3 · Negative binomial regression

7. The negative binomial distribution

In negative binomial regression, \(y_i \mid \mathbf{x}_i \sim \text{NB}(\lambda_i, \theta)\) with

\[ P(y_i = k \mid \mathbf{x}_i) = \frac{\Gamma(k+\theta)}{k!\,\Gamma(\theta)} \left(\frac{\lambda_i}{\lambda_i+\theta}\right)^{k} \left(\frac{\theta}{\lambda_i+\theta}\right)^{\theta}, \qquad k = 0, 1, 2, \ldots \]

Here \(k\) is a possible value of the count, \(\lambda_i\) is the conditional mean, and \(\theta > 0\) controls the spread: \(\text{var}(y_i \mid \mathbf{x}_i) = \lambda_i + \lambda_i^2/\theta\). A small \(\theta\) means a lot of overdispersion, and as \(\theta \to \infty\) the distribution approaches the Poisson. R calls \(\theta\) size. Change theta and run the box again.

8. Fit the negative binomial model

glm.nb() estimates \(\boldsymbol\beta\) and \(\theta\) by maximum likelihood and prints \(\hat\theta\) as Theta. Both AICs below are ordinary likelihood-based AICs for the same 30 responses, so they can be compared.

confint() prints “Waiting for profiling to be done…”; for this model it may also warn that one profiling fit did not converge. After controlling for the other covariates, we expect a one km² increase in island area to be associated with a 0.063% decrease in the mean number of plant species, and a one km² increase in the area of the adjacent island with a 0.076% decrease. The quasi-Poisson and negative binomial results are very similar.

9. Check linearity

For either model, plot the deviance residuals against the fitted values. We want no clear pattern and a least-squares line close to horizontal near 0. A plot with no clear pattern gives no clear evidence against linearity; it cannot prove that the assumption holds.

Optional · Which variance function fits?

The two models differ in how the variance grows with the mean: \(\phi\lambda\) is a straight line and \(\lambda + \lambda^2/\theta\) curves upward. Sort the islands by fitted mean, split them into 5 groups of 6, and compare each group’s average squared residual with both curves. A clear \(\hat\lambda^2\) term in the no-intercept regression would favor the negative binomial.

Part 4 · In-class lab: school absences

10. Load the data and check the Poisson model

For 314 high-school students, daysabs is the number of days absent, a count with no fixed ceiling. The predictors are math, the math score, and prog, the program of study; Academic is the first level and so the reference.

TipQuestion 13.1 (in-class lab): Check, compare, and interpret count models for absences
  1. Check. Give the conclusion of the goodness-of-fit test and of the Cameron and Trivedi test. What does the estimated dispersion say about the conditional variance?
  2. Compare. By what factor are the quasi-Poisson standard errors larger than the Poisson ones?
  3. Interpret. Interpret the negative binomial incidence rate ratio for progVocational and its 95% interval the way the progGeneral model sentence in step 12 does.

11. Fit quasi-Poisson and negative binomial models

Both models account for the overdispersion and give very similar estimates and comparable standard errors, so for concreteness we focus on the negative binomial model. Its much smaller AIC shows that it fits better than the Poisson model; the quasi-Poisson model has no likelihood, so AIC cannot compare it with either.

12. Interpret the negative binomial rate ratios

Model sentence, for progGeneral: after accounting for math score, we expect students in the general program to have, on average, 1.554 times the mean absences of students in the academic program, that is, a 55.39% higher mean number of absences; we are 95% confident that it is between 9.71% and 124.81% higher (\(p = 0.0158\)).

13. Predicted mean absences

The predicted mean number of days absent decreases with math score in all three programs and is lowest for the vocational program at every math score. These are observational data: the model does not show that moving a student into a vocational program would reduce that student’s absences.

After class · Question 13.2

TipQuestion 13.2 (after class): Galapagos plant species, from Poisson to quasi-Poisson and negative binomial
  1. The Pearson residuals of gala.pois have sum of squares \(X^2 = 761.98\), with \(n = 30\) islands and \(p = 5\) slopes. Compute \(\hat\phi = X^2/[n-(p+1)]\) and check that it reproduces summary(gala.quasi)$dispersion. What does its size say about the conditional variance of the number of species?
  2. By what factor are the quasi-Poisson standard errors larger than the Poisson ones? Check it for Elevation with the two outputs. Which covariates are significantly associated with the number of species under the Poisson model but not under the quasi-Poisson model, and which inferences should we trust?
  3. Interpret the negative binomial incidence rate ratio for Elevation and its 95% interval. Then give the multiplicative change in the mean number of species, with its interval, for a 100 m increase in elevation.

Check your understanding

  • Overdispersion means the conditional variance exceeds the conditional mean; Poisson standard errors are then too small and its tests too optimistic.
  • Fit the Poisson model first, then check it with the deviance test, the Cameron and Trivedi test, and \(\hat\phi = X^2/[n-(p+1)]\).
  • Quasi-Poisson keeps the Poisson estimates, multiplies the standard errors by \(\sqrt{\hat\phi}\), and has no likelihood and no ordinary AIC.
  • Negative binomial regression is a full count distribution with \(\text{var}(y_i \mid \mathbf{x}_i) = \lambda_i + \lambda_i^2/\theta\); a small \(\hat\theta\) means a lot of overdispersion.
  • In both models \(e^{\beta_j}\) is an incidence rate ratio for the mean count.
  • Compare likelihood-based models with AIC on the same responses; a quasi-AIC is on a different scale.