Negative Binomial Regression From Gamma Mixing

Why Mix Poisson Rates?

A Poisson regression model assigns a conditional mean to each observation and, under its full distributional specification, sets the conditional variance equal to that mean. Count data can be more variable than this equality permits. One way to represent the extra variation is to allow the Poisson rate to vary across otherwise comparable observations. If that rate variation follows a gamma distribution, averaging over the latent rates produces a negative-binomial distribution.

The construction gives the negative-binomial model a clear interpretation: conditional on observed covariates, observations may differ in their underlying event rates, even when they share the same modelled mean. The derivation below makes both the resulting count probabilities and the source of the additional variance explicit.

The Gamma-Mixed Poisson Model

Let \(Y_i\) be a nonnegative integer count and \(x_i\) the observed covariates for observation \(i\). Let \(\mu_i>0\) be the modelled conditional mean, commonly specified in regression as \(\mu_i=\exp(x_i^\top\beta)\). Introduce a positive latent multiplier \(U_i\). Conditional on both the covariates and that multiplier, assume \(Y_i\) is Poisson with rate \(\mu_iU_i\). The multiplier has mean 1, so it changes an observation's rate without changing the mean rate across observations.

Use the following gamma parameterization. For an overdispersion parameter \(\alpha>0\), let \(U_i\) have a gamma distribution with shape \(1/\alpha\) and scale \(\alpha\). Its density, mean, and variance are:

f_U(u)=\frac{u^{1/\alpha-1}e^{-u/\alpha}}{\Gamma(1/\alpha)\alpha^{1/\alpha}},\quad u>0,\qquad E(U_i)=1,\qquad \operatorname{Var}(U_i)=\alpha.

The scale convention matters: a gamma variable with shape \(r\) and scale \(\theta\) has mean \(r\theta\) and variance \(r\theta^2\). Here those quantities are \(1\) and \(\alpha\), respectively. Equivalently, the gamma rate is \(1/\alpha\).

Definition: A gamma-mixed Poisson regression model assumes \(Y_i\mid x_i,U_i\sim\operatorname{Poisson}(\mu_iU_i)\), with \(U_i\) gamma-distributed with shape \(1/\alpha\) and scale \(\alpha\), independent of \(x_i\), and with a common \(\alpha>0\). After integrating out \(U_i\), \(Y_i\mid x_i\) has a negative-binomial distribution with mean \(\mu_i\) and variance \(\mu_i+\alpha\mu_i^2\).

The independence assumption says that, conditional on included covariates, the distribution of the latent multiplier does not systematically change with those covariates. In a regression application, observations are typically treated as having independent multipliers. This model is a distributional representation of residual rate heterogeneity; by itself, it does not establish why such heterogeneity exists.

Deriving the Negative-Binomial Probabilities

For a fixed covariate value \(x_i\), write \(\mu=\mu_i\), and let \(y\) be a possible count. The conditional Poisson probability given \(U_i=u\) is \(e^{-\mu u}(\mu u)^y/y!\). To obtain the probability given the observed covariates alone, integrate that probability over the gamma density of \(U_i\):

\begin{aligned} P(Y_i=y\mid x_i) &=\int_0^\infty \frac{e^{-\mu u}(\mu u)^y}{y!} \frac{u^{1/\alpha-1}e^{-u/\alpha}}{\Gamma(1/\alpha)\alpha^{1/\alpha}}\,du\\ &=\frac{\mu^y}{y!\,\Gamma(1/\alpha)\alpha^{1/\alpha}} \int_0^\infty u^{y+1/\alpha-1} e^{-(\mu+1/\alpha)u}\,du. \end{aligned}

For \(a>0\) and \(b>0\), the gamma integral is \(\int_0^\infty u^{a-1}e^{-bu}\,du=\Gamma(a)/b^a\). Applying it with \(a=y+1/\alpha\) and \(b=\mu+1/\alpha\), then collecting the terms, gives:

P(Y_i=y\mid x_i)= \frac{\Gamma(y+1/\alpha)}{\Gamma(1/\alpha)y!} \left(\frac{1}{1+\alpha\mu_i}\right)^{1/\alpha} \left(\frac{\alpha\mu_i}{1+\alpha\mu_i}\right)^y, \qquad y=0,1,2,\ldots.

This is a negative-binomial probability mass function. The expression uses the mean-and-overdispersion parameterization: \(\mu_i\) is the conditional mean, and \(\alpha\) controls the additional variation. An equivalent shape parameter is \(r=1/\alpha\). Parameterizations differ across software, so identifying the parameter convention is essential when interpreting a reported negative-binomial fit.

Worked Example: Calculating Probabilities From the Mixture

Suppose a hypothetical observation has conditional Poisson mean \(\mu=2\), and suppose the gamma multiplier has \(\alpha=0.5\). The assumptions are that the count is Poisson conditional on its multiplier, and that this multiplier has gamma shape \(1/\alpha=2\) and scale \(\alpha=0.5\). Calculate the probabilities of counts 0, 1, and 2.

Here \(1/\alpha=2\), while \(\alpha\mu=1\), so the two probability factors in the negative-binomial formula are both \(1/2\). For \(y=0\):

P(Y=0\mid x)= \frac{\Gamma(2)}{\Gamma(2)0!} \left(\frac12\right)^2 \left(\frac12\right)^0 =1(0.25)(1)=0.25.

For \(y=1\), the gamma-function coefficient is \(\Gamma(3)/(\Gamma(2)1!)=2\):

P(Y=1\mid x)=2\left(\frac12\right)^2\left(\frac12\right) =2(0.25)(0.5)=0.25.

For \(y=2\), the coefficient is \(\Gamma(4)/(\Gamma(2)2!)=6/2=3\):

P(Y=2\mid x)=3\left(\frac12\right)^2\left(\frac12\right)^2 =3(0.25)(0.25)=0.1875.

Thus the probabilities are 0.250, 0.250, and 0.188, respectively. These are probabilities of the specified count values under the mixed model; they are not estimates of the probability that the observation has a particular latent multiplier.

Conditional Mean and Variance

The moments follow directly from conditioning on the multiplier. For a Poisson variable, its conditional mean and variance both equal its rate. The law of total expectation therefore gives:

E(Y_i\mid x_i)=E\{E(Y_i\mid x_i,U_i)\mid x_i\} =E(\mu_iU_i\mid x_i)=\mu_i.

For the variance, the law of total variance separates the average conditional Poisson variance from the variance in the conditional mean across latent multipliers:

\begin{aligned} \operatorname{Var}(Y_i\mid x_i) &=E\{\operatorname{Var}(Y_i\mid x_i,U_i)\mid x_i\} +\operatorname{Var}\{E(Y_i\mid x_i,U_i)\mid x_i\}\\ &=E(\mu_iU_i\mid x_i)+\operatorname{Var}(\mu_iU_i\mid x_i)\\ &=\mu_i+\mu_i^2\alpha. \end{aligned}

The first term, \(\mu_i\), is the average Poisson variation within a given latent rate. The second term, \(\alpha\mu_i^2\), is the extra variation caused by the changing rates. Since \(\alpha>0\), the marginal conditional variance exceeds the mean for every \(\mu_i>0\). When \(\alpha\) tends to zero, the gamma multiplier becomes increasingly concentrated at 1 and the distribution approaches the Poisson model.

Theorem: Under the gamma-mixed Poisson assumptions, \(E(Y_i\mid x_i)=\mu_i\) and \(\operatorname{Var}(Y_i\mid x_i)=\mu_i+\alpha\mu_i^2\). The parameter \(\alpha\) is the variance of the mean-one latent multiplier, so it represents the extent of latent rate heterogeneity in this model.

Worked Example: Reading the Overdispersion Parameter

Suppose a hypothetical count regression predicts a conditional mean of \(\mu_i=4.5\) for an observation and uses \(\alpha=0.2\). Assume the gamma-mixed Poisson specification above, including a mean-one multiplier independent of the covariates. Calculate the conditional variance and describe what \(\alpha\) says about the multiplier.

Substitute the predicted mean and overdispersion parameter into the variance formula:

\operatorname{Var}(Y_i\mid x_i) =4.5+0.2(4.5)^2 =4.5+0.2(20.25) =4.5+4.05 =8.55.

The variance is 8.55 counts squared, compared with a Poisson variance of 4.5 counts squared at the same mean. The model attributes the additional \(4.05\) to variation in the latent rate. Since \(\operatorname{Var}(U_i)=0.2\), the multiplier's standard deviation is \(\sqrt{0.2}\approx0.447\); its coefficient of variation is also about 0.447 because \(E(U_i)=1\). This describes the assumed distribution of latent multipliers, not the variation in the observed count alone.

Why the Mixing Interpretation Matters

The mixing construction connects a familiar conditional count model to a marginal distribution that allows extra variation. In a regression, covariates determine \(\mu_i\), often through a log link, while \(\alpha\) describes how much additional conditional variance is associated with gamma-distributed rate variation. The quadratic term means that the size of the extra variance depends on the mean: it is \(\alpha\mu_i^2\), not a fixed amount for every observation.

The parameter also gives a compact description of the latent multiplier's spread. A larger \(\alpha\) means a more dispersed multiplier distribution and therefore more variation in Poisson rates among observations with the same covariates. This can help explain why a Poisson variance assumption may be too restrictive. It does not, however, prove that gamma-distributed unobserved heterogeneity is the true data-generating mechanism. Other mechanisms can also produce extra variation, so the mixing model is an assumption with an interpretable consequence.

A useful way to check the algebra is to compare the variance-to-mean ratio with one:

\frac{\operatorname{Var}(Y_i\mid x_i)}{E(Y_i\mid x_i)} =\frac{\mu_i+\alpha\mu_i^2}{\mu_i} =1+\alpha\mu_i.

For positive \(\alpha\) and \(\mu_i\), this ratio exceeds one. It also changes with \(\mu_i\), even when \(\alpha\) is common across observations. The overdispersion parameter is therefore not itself the variance-to-mean ratio; rather, it determines how that ratio increases with the conditional mean.

Worked Example: Separating Poisson Noise From Rate Heterogeneity

Consider a hypothetical observation with \(\mu_i=3\) and \(\alpha=1/3\). Its multiplier has gamma shape \(1/\alpha=3\) and scale \(\alpha=1/3\), so its mean is 1 and variance is \(1/3\). Verify the total variance by calculating the two components in the law of total variance.

The expected conditional Poisson variance is:

E\{\operatorname{Var}(Y_i\mid x_i,U_i)\mid x_i\} =E(3U_i\mid x_i)=3E(U_i)=3.

The conditional Poisson mean is \(3U_i\), whose variance across multipliers is:

\operatorname{Var}\{E(Y_i\mid x_i,U_i)\mid x_i\} =\operatorname{Var}(3U_i\mid x_i) =3^2\operatorname{Var}(U_i) =9\left(\frac13\right)=3.

Adding the components gives \(\operatorname{Var}(Y_i\mid x_i)=3+3=6\), which agrees with \(3+(1/3)(3)^2=6\). Half of the total variance in this example comes from the average conditional Poisson variation and half from variation in the latent rates. This decomposition is a property of the assumed model, not a direct observation of each component for an individual count.

Common Pitfall

A common mistake is to call \(\alpha\) the variance-to-mean ratio or to interpret \(\alpha\) as the amount added directly to every observation's variance. Neither is correct. The excess variance over the conditional mean is \(\alpha\mu_i^2\), and the variance-to-mean ratio is \(1+\alpha\mu_i\). Both depend on the observation's conditional mean, while \(\alpha\) is the variance of the mean-one latent multiplier under the gamma-mixing specification.

A second caution is interpretive: the mixing representation does not establish that a specific omitted variable or mechanism caused the observed excess variation. It states how a gamma distribution of latent rates would generate the negative-binomial count distribution. That interpretation is useful, but it remains tied to the model's assumptions.

Key takeaway: Mixing Poisson rates with a mean-one gamma multiplier produces a negative-binomial conditional distribution. Its mean remains \(\mu_i\), while its variance becomes \(\mu_i+\alpha\mu_i^2\); \(\alpha\) measures the assumed latent rate heterogeneity.

Check Your Understanding

Use the gamma-mixed Poisson assumptions and parameterization from this tutorial.

For \(\alpha=0.25\), what are the shape, scale, mean, and variance of the latent multiplier?
Starting from the law of total variance, identify which term represents average conditional Poisson variation and which represents variation across latent rates.
If \(\mu_i=5\) and \(\alpha=0.1\), calculate the conditional variance and variance-to-mean ratio.
Explain why \(\alpha\) is not itself the excess variance for an observation with conditional mean \(\mu_i\).
What does the gamma-mixing interpretation assume about the relationship between the latent multiplier and observed covariates?