Negative Binomial Models for Overdispersed Treatment Counts

Why Model Overdispersed Treatment Counts?

A count-valued treatment might record the number of doses, contacts, or service units assigned to a person. A Poisson model is a natural starting point, but it requires the conditional variance of the count to equal its conditional mean. When treatment counts vary more than that restriction allows, fitted Poisson probabilities can misrepresent the chance of the observed count. That matters for inverse probability weighting: the denominator probability for a person’s observed treatment must reflect the treatment mechanism well enough to avoid systematically extreme weights.

The negative binomial model relaxes the Poisson variance restriction while retaining a count distribution on the nonnegative integers. The focus here is its regression specification, fitting by likelihood, and calculation of the probability assigned to an observed count. The earlier material on Poisson models for count-valued treatments gives the corresponding equidispersed model.

Definition and Parameterization

Definition: Under the negative binomial mean-dispersion parameterization, a count $A$ has mean $\mu>0$ and dispersion $\alpha>0$ if, with $r=1/\alpha$, $$ P(A=a)=\frac{\Gamma(a+r)}{\Gamma(r)\,\Gamma(a+1)} \left(\frac{r}{r+\mu}\right)^r \left(\frac{\mu}{r+\mu}\right)^a,\qquad a=0,1,2,\ldots. $$ Its mean is $\mu$ and its variance is $\mu+\alpha\mu^2$. The parameter $r$ is sometimes called the size parameter; $\alpha$ is the overdispersion parameter. This tutorial uses $\alpha$.

For treatment regression, let $A_i$ be the observed count and $L_i$ the covariate vector for individual $i$. A common specification is $A_i\mid L_i \sim \operatorname{NB}(\mu_i,\alpha)$, with $\log(\mu_i)=L_i^\top\beta$. The log link ensures $\mu_i>0$; $\alpha$ is usually shared across individuals. More elaborate models can allow dispersion to depend on covariates, but that adds parameters and needs adequate data.

This model is appropriate when the outcome is a nonnegative integer and conditional variability exceeds the Poisson level. It does not make every count distribution plausible: the fitted mean structure, dispersion pattern, and tail probabilities still need to be checked. In a treatment-weight application, also require that the model assign positive probability to the observed count at the covariate values in question. A probability that is positive mathematically but extremely small can still produce a highly influential weight.

Why the Variance Exceeds the Mean

A useful derivation represents negative binomial variation as a Poisson rate that varies across individuals. Suppose, conditionally on a latent rate $\Lambda$, $A\mid\Lambda\sim\operatorname{Poisson}(\Lambda)$, and suppose $\Lambda\sim\operatorname{Gamma}(r,\text{rate}=r/\mu)$. Thus $E[\Lambda]=\mu$ and $\operatorname{Var}(\Lambda)=\mu^2/r$.

Result: Integrating out this Gamma-distributed rate gives the negative binomial probability mass function above, with $\alpha=1/r$, $E[A]=\mu$, and $\operatorname{Var}(A)=\mu+\alpha\mu^2$.

To verify the probability mass function, multiply the Poisson probability by the Gamma density and integrate over $\lambda>0$:

\begin{aligned} P(A=a) &=\int_0^\infty \frac{e^{-\lambda}\lambda^a}{\Gamma(a+1)} \frac{(r/\mu)^r}{\Gamma(r)}\lambda^{r-1}e^{-r\lambda/\mu}\,d\lambda\\ &=\frac{(r/\mu)^r}{\Gamma(a+1)\Gamma(r)} \frac{\Gamma(a+r)}{(1+r/\mu)^{a+r}}\\ &=\frac{\Gamma(a+r)}{\Gamma(a+1)\Gamma(r)} \left(\frac{r}{r+\mu}\right)^r \left(\frac{\mu}{r+\mu}\right)^a. \end{aligned}

The last integral uses the Gamma integral $\int_0^\infty \lambda^{c-1}e^{-b\lambda}\,d\lambda=\Gamma(c)/b^c$, for $b,c>0$. The mean and variance follow from conditional expectation and variance:

\begin{aligned} E[A]&=E\{E[A\mid\Lambda]\}=E[\Lambda]=\mu,\\ \operatorname{Var}(A) &=E\{\operatorname{Var}(A\mid\Lambda)\} +\operatorname{Var}\{E[A\mid\Lambda]\}\\ &=E[\Lambda]+\operatorname{Var}(\Lambda) =\mu+\frac{\mu^2}{r} =\mu+\alpha\mu^2. \end{aligned}

The mixture interpretation explains the extra variance: individuals with the same measured covariates can still have different underlying rates. It is a probability model, not proof that a literal latent Gamma rate exists in the treatment process. When $\alpha$ approaches zero, the extra term vanishes and the distribution approaches a Poisson distribution with mean $\mu$.

Worked Example: Estimating Dispersion From Count Moments

Suppose six observed treatment counts are $0,0,1,1,2,4$. Compute the empirical mean and variance using divisor $n$, then use the method of moments to estimate $\alpha$.

The count sum is $8$, so $\bar A=8/6=4/3$. The squared-count sum is $0+0+1+1+4+16=22$. The empirical second moment is $22/6=11/3$, giving empirical variance $11/3-(4/3)^2=11/9$ (about $1.222$). This particular sample does not show empirical overdispersion: $11/9<4/3$.

The method-of-moments equation for the model is $s_n^2=\bar A+\alpha\bar A^2$. Solving gives $\hat\alpha=(s_n^2-\bar A)/\bar A^2$, which here is $(11/9-4/3)/(16/9)=(-1/9)(9/16)=-1/16$. A negative value is not allowed by the negative binomial model. Thus these six counts do not yield a positive method-of-moments dispersion estimate. This does not establish that the population is underdispersed; it illustrates that small samples can give noisy moment estimates. Likelihood fitting constrains $\alpha$ to be positive, and a boundary estimate near zero corresponds to the Poisson limit.

Likelihood Fitting and the Mean-Parameter Score

For independent individuals conditional on their covariates, the conditional likelihood is the product of their negative binomial probabilities. Write $r=1/\alpha$ and $\mu_i=\exp(L_i^\top\beta)$. The log-likelihood is

\ell(\beta,\alpha)= \sum_{i=1}^n\left[ \log\Gamma(A_i+r)-\log\Gamma(r)-\log\Gamma(A_i+1) +r\log\left(\frac{r}{r+\mu_i}\right) +A_i\log\left(\frac{\mu_i}{r+\mu_i}\right) \right].

One generally maximizes this likelihood numerically over $\beta$ and $\alpha>0$. For fixed $\alpha$, differentiating with respect to the linear predictor $\eta_i=L_i^\top\beta=\log\mu_i$ yields the score contribution $(A_i-\mu_i)/(1+\alpha\mu_i)$. Therefore the score for $\beta$, holding $\alpha$ fixed, is as follows.

Result: For the log-mean negative binomial regression with fixed $\alpha>0$, the score is $$ \frac{\partial\ell}{\partial\beta} =\sum_{i=1}^n L_i\frac{A_i-\mu_i}{1+\alpha\mu_i}. $$

Indeed, the terms depending on $\eta_i$ are $-r\log(r+\mu_i)+A_i\log\mu_i-A_i\log(r+\mu_i)$, up to constants. Since $d\mu_i/d\eta_i=\mu_i$, their derivative is $A_i-\mu_i(r+A_i)/(r+\mu_i)=r(A_i-\mu_i)/(r+\mu_i)$, equal to $(A_i-\mu_i)/(1+\alpha\mu_i)$. Multiplying by $\partial\eta_i/\partial\beta=L_i$ and summing proves the result. $\square$

This score shows how the model uses covariate patterns to explain conditional means; estimating $\alpha$ additionally adjusts the likelihood for residual variation. In practice, inspect convergence, fitted means, the dispersion estimate, and observed versus predicted count frequencies. A fitted model can converge while fitting the tail poorly, so convergence alone is not evidence of an adequate treatment model.

Worked Example: Calculating Conditional Observed-Count Probabilities

Suppose a fitted model has $\log(\mu_i)=\beta_0+\beta_1x_i$, with $\beta_0=\log(1.2)$, $\beta_1=\log(2)$, and $\alpha=0.5$. Calculate the probability of observing $A_i=3$ for $x_i=0$ and $x_i=1$.

Here $r=1/\alpha=2$. At $x_i=0$, $\mu_i=1.2$, so $r/(r+\mu_i)=2/3.2=5/8$ and $\mu_i/(r+\mu_i)=3/8$. Thus $P(A_i=3\mid x_i=0)=\binom{4}{3}(5/8)^2(3/8)^3 =2700/32768\approx0.0824$.

At $x_i=1$, $\mu_i=1.2\times2=2.4$, giving $r/(r+\mu_i)=5/11$ and $\mu_i/(r+\mu_i)=6/11$. Therefore $P(A_i=3\mid x_i=1)=\binom{4}{3}(5/11)^2(6/11)^3 =21600/161051\approx0.1341$.

The same count has different probabilities because the conditional mean differs. These are probabilities for the observed count under the fitted denominator model; they are not the probability of receiving any treatment at all.

From a Count Probability to a Treatment Weight

For a count-valued treatment $A_i$ and covariates $L_i$, a fitted negative binomial model supplies the denominator probability $\hat f(A_i\mid L_i)$: evaluate the fitted mass function at the actual count $A_i$ and covariates $L_i$. A stabilized weight may use a numerator probability as well, provided that numerator is specified for the target weighting procedure. The negative binomial model only supplies probabilities for whichever conditional distribution it is used to model.

Worked Example: A Simple Probability Ratio

Continue with the fitted denominator model in the previous example. For an individual with $x_i=1$ and observed count $A_i=3$, the denominator probability is approximately $0.1341$. Suppose, just for illustration, a specified numerator distribution assigns probability $0.1028$ to this observed count. The probability ratio is $0.1028/0.1341\approx0.767$.

This calculation uses the probability of the observed count in each distribution. It does not replace the need to justify the numerator model or the causal assumptions underlying the weights. If the denominator probability were $0.001$ instead, the same numerator probability would produce a ratio of $102.8$, illustrating why tiny fitted probabilities deserve careful investigation.

Diagnostics, Alternatives, and Limitations

Compare the conditional observed variance with the fitted Poisson mean restriction and examine whether the negative binomial model improves calibration across the count range. Useful checks include tabulating observed and model-predicted frequencies for counts, comparing observed and fitted tail frequencies, and inspecting fitted probabilities for the counts actually observed. Because covariates define conditional distributions, a single overall variance comparison is only a first screen; check patterns across important covariate strata or fitted-mean ranges as well.

If fitted probabilities for observed counts are very small, the resulting inverse probabilities can be large. Look for sparse covariate patterns, influential observations, poor tail fit, and model misspecification before treating the problem as merely a numerical one. Restricting or truncating weights changes the weighting procedure and does not repair a misspecified treatment model. Likewise, adding many predictors can make conditional probabilities unstable when the data do not support that detail.

The Poisson model is simpler and is recovered as $\alpha$ approaches zero; it can be more efficient when its variance restriction is credible. The negative binomial model is useful when extra-Poisson variation is present, but it is not automatically preferable: a poor mean model or inappropriate tail shape can still give poor probabilities. Empirical cell probabilities are another possibility for a small set of discrete treatment values, but become unreliable when covariate-defined cells are sparse. Flexible prediction methods offer alternatives, while requiring their own checks. These alternatives answer related estimation problems; the key requirement here is that the fitted count distribution represent the probability of the observed assignment.

Check Your Understanding

Answer each question, then compare with the solutions below.

In the mean-dispersion parameterization, what is the variance when $\mu=3$ and $\alpha=0.4$?
For $\mu=2$ and $\alpha=0.5$, calculate $P(A=0)$ and $P(A=1)$.
Why does a negative method-of-moments estimate of $\alpha$ not define a negative binomial model?
Name two checks to make if a fitted model assigns a very small probability to an observed treatment count.

Solutions. (1) $\operatorname{Var}(A)=3+0.4(3^2)=6.6$. (2) Here $r=2$, so $P(A=0)=(2/4)^2=0.25$. Also $P(A=1)=\binom{2}{1}(1/2)^2(1/2)=0.25$. (3) The model requires $\alpha>0$; a negative estimate is an invalid moment estimate, not a permissible dispersion value. (4) Check observed-versus-predicted count frequencies, especially the tail; inspect the covariate pattern for sparse support or influential observations. Also verify model convergence and the fitted mean structure.