Negative-Binomial Predictive Distributions for Future Counts

Prediction Has Two Sources of Variation

A fitted negative-binomial regression can answer two different predictive questions. If its estimated coefficients and dispersion are treated as fixed, the model gives a distribution for a future count at specified covariates. If uncertainty in those estimates is also represented, the predictive distribution must average over possible conditional means and, when appropriate, possible dispersion values. These are not the same distribution.

The distinction matters in clinical planning. A fitted mean might be the expected number of events for a patient group, but an actual future count can differ from that mean because counts vary even when the model parameters are known. Uncertainty in the fitted mean adds another source of variation. A prediction interval for a future count should account for the first source; an interval intended to reflect parameter uncertainty should account for the second as well.

Definition: Under the mean–size parameterization of the NB2 model, a count $Y$ with conditional mean $\mu>0$ and size $k>0$ has probability mass function $$ P(Y=y\mid \mu,k)=\frac{\Gamma(y+k)}{\Gamma(k)y!} \left(\frac{k}{k+\mu}\right)^k \left(\frac{\mu}{k+\mu}\right)^y,\qquad y=0,1,2,\ldots . $$ Its conditional mean and variance are $E[Y\mid\mu,k]=\mu$ and $\operatorname{Var}(Y\mid\mu,k)=\mu+\mu^2/k$. The size $k$ is a dispersion parameter: smaller $k$ means greater variance at a fixed mean.

For a regression with coefficient vector $\beta$, a new covariate vector $x_{\mathrm{new}}$, and any relevant offset $o_{\mathrm{new}}$, the linear predictor and conditional mean are $\eta_{\mathrm{new}}=x_{\mathrm{new}}^\top\beta+o_{\mathrm{new}}$ and $\mu_{\mathrm{new}}=\exp(\eta_{\mathrm{new}})$. The usual plug-in prediction substitutes the fitted values $\hat\beta$ and $\hat k$, giving $\hat\mu_{\mathrm{new}}=\exp(x_{\mathrm{new}}^\top\hat\beta+o_{\mathrm{new}})$. The following examples first treat those estimates as fixed, then show what changes when the mean is uncertain.

Plug-In Predictive Probabilities and Intervals

With fitted parameters held fixed, the predictive distribution for a future count is the negative-binomial distribution evaluated at $\hat\mu_{\mathrm{new}}$ and $\hat k$. For example, the probability of no events is $\{\hat k/(\hat k+\hat\mu_{\mathrm{new}})\}^{\hat k}$. More generally, a probability such as $P(Y_{\mathrm{new}}\ge c)$ is obtained by summing the mass function over all integer counts at least $c$, or by subtracting the corresponding cumulative probability from one.

A prediction interval for a count is an interval of integers. One conventional equal-tail interval uses the quantile $q_\alpha=\min\{y\in\{0,1,2,\ldots\}:F(y)\ge\alpha\}$, where $F$ is the predictive cumulative distribution function. The interval $[q_{0.025},q_{0.975}]$ is a nominal 95% equal-tail interval. Since counts are discrete, its actual coverage under the specified distribution can exceed 95%; its endpoints need not leave exactly 2.5% probability in each tail.

Worked Example: A Plug-In 95% Prediction Interval

Suppose a fitted model gives $\hat\mu_{\mathrm{new}}=4$ and $\hat k=2$. Calculate the probability of zero events and the equal-tail 95% prediction interval, treating the fitted values as fixed.

The probability mass function has $k/(k+\mu)=2/6=1/3$ and $\mu/(k+\mu)=4/6=2/3$. Thus $P(Y_{\mathrm{new}}=y)=\frac{y+1}{9}(2/3)^y$, and in particular $P(Y_{\mathrm{new}}=0)=(1/3)^2=1/9\approx0.1111$.

The cumulative probabilities at the relevant endpoints are $F(0)=1/9\approx0.1111$ and $F(1)=1/9+2/27=5/27\approx0.1852$. Therefore $q_{0.025}=0$. For the upper endpoint, using the negative-binomial tail formula for size 2, $P(Y>y)=(2/3)^{y+1}\{(y+2)-(2/3)(y+1)\}$, gives $F(12)\approx0.9726<0.975$ and $F(13)\approx0.9806\ge0.975$. Hence $q_{0.975}=13$, and the equal-tail interval is $[0,13]$.

This interval describes future-count variation conditional on the fitted mean and size being the values used in the calculation. It does not, by itself, incorporate uncertainty in estimating those values.

Averaging Over Uncertainty in the Conditional Mean

To represent uncertainty in the conditional mean, regard $M$ as an uncertain positive quantity. Conditional on $M$ and a fixed size $k$, suppose the future count follows the negative-binomial mass function above with $\mu=M$. The predictive distribution after accounting for uncertainty in $M$ is a mixture: for each count, average its conditional probability over the distribution assigned to $M$. This distribution may come from a Bayesian posterior, a specified approximation to parameter uncertainty, or a simulation procedure. The interpretation depends on how that distribution was constructed.

Result: If $E[M]$ and $E[M^2]$ are finite, and $Y\mid M$ has the negative-binomial distribution with fixed size $k$ and mean $M$, then $$ E[Y]=E[M],\qquad \operatorname{Var}(Y)=E[M]+\frac{E[M^2]}{k}+\operatorname{Var}(M). $$ Equivalently, $$ \operatorname{Var}(Y)=E[M]+\frac{(E[M])^2}{k} +\left(1+\frac{1}{k}\right)\operatorname{Var}(M). $$

Derivation. Conditional expectation gives $E[Y]=E\{E[Y\mid M]\}=E[M]$. By the law of total variance, $\operatorname{Var}(Y)=E\{\operatorname{Var}(Y\mid M)\}+\operatorname{Var}\{E[Y\mid M]\}$. Substituting the conditional mean and variance gives $E[M+M^2/k]+\operatorname{Var}(M)=E[M]+E[M^2]/k+\operatorname{Var}(M)$. Finally, $E[M^2]=(E[M])^2+\operatorname{Var}(M)$, which yields the equivalent expression. Thus the conditional count variance is not the whole predictive variance: variation in $M$ contributes as well. $\square$

A useful approximate representation arises when the uncertain linear predictor is normal. If $\eta\sim N(m,v)$ and $M=\exp(\eta)$, then $M$ is lognormal and $E[M]=\exp(m+v/2)$, $E[M^2]=\exp(2m+2v)$, and $\operatorname{Var}(M)=\exp(2m+v)(\exp(v)-1)$. These identities follow from the normal moment-generating function $E[\exp(t\eta)]=\exp(tm+t^2v/2)$, evaluated at $t=1$ and $t=2$. If $m=\hat\eta_{\mathrm{new}}$ and $v$ is derived from coefficient uncertainty, the resulting moments illustrate how uncertainty in the fitted predictor affects future-count variation. They do not automatically define an exact frequentist predictive distribution; the source and interpretation of the uncertainty distribution must be stated.

Worked Example: Predictive Moments Under Lognormal Mean Uncertainty

Suppose the uncertain linear predictor is $\eta\sim N(\log 4,0.04)$, and the size is fixed at $k=2$. Find the predictive mean and variance of a future count.

The conditional mean is $M=\exp(\eta)$. Its first two moments are $E[M]=\exp(\log 4+0.04/2)=4\exp(0.02)\approx4.0808$ and $E[M^2]=\exp(2\log 4+2(0.04))=16\exp(0.08)\approx17.3326$. Therefore $\operatorname{Var}(M)=16\exp(0.08)-\{4\exp(0.02)\}^2 =16\exp(0.08)-16\exp(0.04)\approx0.6796$.

The predictive mean is $E[Y]=E[M]\approx4.0808$. The variance is $E[M]+E[M^2]/2+\operatorname{Var}(M)$, so $\operatorname{Var}(Y)\approx4.0808+17.3326/2+0.6796=13.4267$, or approximately $13.43$. If uncertainty in $M$ were omitted, the plug-in variance at mean 4 would be $4+4^2/2=12$. The predictive variance is larger because it includes variation in the conditional mean as well as the negative-binomial count variation.

Mixture Probabilities Are Not Plug-In Probabilities

When the mean is uncertain, averaging the mean and then evaluating a negative-binomial probability generally does not give the mixture probability. For the event $Y=0$, the conditional probability is $\{k/(k+M)\}^k$, a nonlinear function of $M$. Its average over uncertainty in $M$ need not equal the same function evaluated at $E[M]$. The general mixture mass function makes the distinction explicit.

Result: For an uncertain positive mean $M$ and fixed $k>0$, the mixture predictive mass function is $$ P(Y=y)=E\!\left[ \frac{\Gamma(y+k)}{\Gamma(k)y!} \left(\frac{k}{k+M}\right)^k \left(\frac{M}{k+M}\right)^y \right],\qquad y=0,1,2,\ldots . $$ It is a probability mass function: each integrand is nonnegative, and summing over $y$ gives $E[\sum_{y=0}^\infty P(Y=y\mid M,k)]=E[1]=1$. Interchanging expectation and the nonnegative sum is justified by Tonelli’s theorem. Its cumulative probabilities and integer quantiles can therefore be computed by summing these mixture masses.

Worked Example: A Mixture Changes the Zero-Count Probability

For a simple illustration of mean uncertainty, let $k=2$, and suppose $M$ equals 2 or 4 with probability $1/2$ each. Calculate the mixture probability of zero events, and compare it with the plug-in probability at the average mean.

At $M=2$, the zero-count probability is $\{2/(2+2)\}^2=1/4=0.25$. At $M=4$, it is $\{2/(2+4)\}^2=1/9\approx0.1111$. Averaging these conditional probabilities gives $P(Y=0)=\frac12(1/4)+\frac12(1/9)=13/72\approx0.1806$.

The average mean is $E[M]=3$. A plug-in calculation at that mean instead gives $\{2/(2+3)\}^2=4/25=0.16$, which is not the mixture probability. The difference occurs because the zero-count probability is nonlinear in the mean.

The variance decomposition provides a check on the mixture’s moments. Here $E[M]=3$, $E[M^2]=(4+16)/2=10$, and $\operatorname{Var}(M)=10-9=1$. Thus $E[Y]=3$ and $\operatorname{Var}(Y)=3+10/2+1=9$. This agrees with averaging the conditional variances, $\frac12(2+2^2/2)+\frac12(4+4^2/2)=8$, and adding the variance of the conditional mean, 1.

From Regression Uncertainty to a Predictive Distribution

In a fitted regression, a common first approximation to uncertainty in the linear predictor is obtained from the estimated covariance matrix of $\hat\beta$. If that covariance matrix is $\widehat{\operatorname{Var}}(\hat\beta)$, then the estimated variance of $x_{\mathrm{new}}^\top\hat\beta+o_{\mathrm{new}}$, treating the offset as fixed, is $v_{\mathrm{new}}=x_{\mathrm{new}}^\top\widehat{\operatorname{Var}}(\hat\beta)x_{\mathrm{new}}$. Under the usual large-sample normal approximation for the coefficient estimator, one can draw predictor values around the fitted linear predictor, transform them to positive means, and draw future counts conditional on those means. If the dispersion estimate is also treated as uncertain, its uncertainty must be included in the draws or in a suitable joint approximation.

Obtain a distribution for fitted parameters.
Use a justified approximation or simulation reflecting the coefficient and, if required, dispersion uncertainty. State whether this represents posterior uncertainty or an approximate sampling distribution.

Calculate the conditional mean.
For each parameter draw, evaluate the new observation’s linear predictor, including its offset, and exponentiate to obtain a positive mean.

Generate a future count.
For each drawn mean and size, generate a negative-binomial count. Alternatively, average the conditional mass probabilities when calculating a particular event probability.

Summarize the predictive distribution.
Estimate event probabilities or take empirical integer quantiles of the simulated counts. Do not substitute a standard error for the mean directly as if it were the count’s prediction interval.

This procedure separates two operations: drawing uncertain model quantities and drawing the outcome conditional on them. A confidence interval for the regression mean addresses uncertainty in a conditional mean, not the full variation of a future count. Conversely, the plug-in negative-binomial interval includes count variation under the fitted model but holds the estimated parameters fixed. The intended target determines which distribution is appropriate.

Common Mistake: Confusing Mean Intervals with Count Prediction Intervals

A frequent error is to report $\hat\mu_{\mathrm{new}}$ plus or minus a standard error as an interval for the future count. That calculation concerns uncertainty in an estimated mean on its chosen scale; it does not include the negative-binomial variation of individual counts and may even produce non-integer or negative endpoints. The correction is to specify the target first: use an interval for the conditional mean when inference about the mean is wanted, and use the negative-binomial predictive distribution for future counts. If parameter uncertainty is also part of the prediction target, integrate or simulate over it and then include the conditional count draw.

Another mistake is to treat the overdispersion parameter as if it were a standard error for the fitted mean. In the NB2 variance $\mu+\mu^2/k$, $k$ governs count variability conditional on the mean. Coefficient uncertainty governs uncertainty in the regression predictor and hence in the mean. They enter predictive calculations differently, even though both can increase the overall spread of predictions.

Key takeaway: A negative-binomial plug-in prediction describes count variation conditional on fitted parameters. A prediction that also represents parameter uncertainty is a mixture over possible means (and possibly sizes); its probabilities and intervals should be computed from that mixture, not from a confidence interval for the mean alone.

Check Your Understanding

Use these questions to check that you can distinguish the predictive targets and carry out the main calculations.

For a future count with fixed mean $\mu$ and size $k$, write its zero-count probability and explain what variation its predictive distribution includes.
Why can a mixture probability $E[\{k/(k+M)\}^k]$ differ from the plug-in probability $\{k/(k+E[M])\}^k$?
State the law-of-total-variance expression for a negative-binomial count with uncertain mean $M$ and fixed size $k$.
In a fitted regression, what role does $x_{\mathrm{new}}^\top\widehat{\operatorname{Var}}(\hat\beta)x_{\mathrm{new}}$ play, and what source of variation does it not by itself represent?
Explain why a confidence interval for the conditional mean is not automatically a prediction interval for a future count.