Overdispersed Poisson regression

Question

In Gelman and Hill 2007 (http://www.stat.columbia.edu/~gelman/arm/), they mention that adding an $\epsilon \sim N(0,\sigma^2)$ term to a Poisson regression can be used to account for overdispersion.

On a number of different sites I've seen this mentioned, but all refer you back to Gelman and Hill for an explanation of this, and the text doesn't address conceptually why this makes sense (although I believe them).

I understand that a two-parameter distribution (eg negative-binomial) may be used in place of a Poisson, and this makes sense to me: it has an extra parameter to capture variance. However, it is not at all obvious to me why simply adding an $\epsilon$ term should capture overdispersion, as the resulting Poisson will still have mean equal to its variance (as always). Can anyone clarify what's going on here and why this should account for overdispersion, as an alternative to using a negative binomial? And is there any simple way of interpreting how much overdispersion this $\epsilon$ should capture?

Answer 1

The "resulting Poisson" is only the conditional distribution of the response given an (unobserved) realization of the normal error term: the unconditional response (or the response conditional only on the realized values of the predictors, if these are random variables) does not have a Poisson distribution.

Suppose the response $Y$ has a Poisson distribution

$$\newcommand{\E}{\operatorname{E}}\newcommand{\var}{\operatorname{Var}}Y|M \sim \mathrm{Pois}(M)$$

where $M$ is a log-normal variate with log-location given by the sum of products of each predictor $x_i$ with its coefficient $\beta_i$, & with log-scale $\sigma$:

$$M = \exp(\beta_0 + \beta_1 x_1 + \ldots + \sigma Z)$$

Conditionally on the realized value of $M$, $m$, $Y$ indeed has variance equal to the mean

$$\E (Y|M=m) = \var (Y | M=m) = m$$

but unconditionally it doesn't: the unconditional expectation is

$$\E Y = \E \E Y|M = \E M \\ = \exp(\beta_0 + \beta_1 x_1 + \ldots + \tfrac{\sigma^2}{2}) \\ = \exp(\beta_0 + \beta_1 x_1 + \ldots) \exp(\tfrac{\sigma^2}{2})$$

& the unconditional variance is

$$\var Y = \E \var Y|M + \var \E Y|M = \E M + \var M \\ = \exp(\beta_0 + \beta_1 x_1 + \ldots + \tfrac{\sigma^2}{2}) + [\exp(\sigma^2)-1]\exp(2(\beta_0 + \beta_1 x_1 + \ldots + \tfrac{\sigma^2}{2}))\\ = \exp(\beta_0 + \beta_1 x_1 + \ldots) \exp(\tfrac{\sigma^2}{2}) + \exp(\beta_0 + \beta_1 x_1 + \ldots)^2 [\exp(\sigma^2)-1] \exp(\sigma^2)$$

increasingly larger than the mean as $\sigma$ increases.

Over-dispersion is measured by $\sigma$: as it tends to zero the model tends to a Poisson model without over-dispersion. Note that over-dispersion parameters in different model families do not bear precisely the same interpretation. There is more than one way to specify a negative binomial model: typically from

$$\var Y = \E Y + \alpha (\E Y)^p$$

where $p$=1 for the NB1 model, $p=2$ for the NB2 model; note the over-dispersion parameter $\alpha$ relates variance to mean differently for different values of $p$.