Estimating and Interpreting Negative Binomial Effects

Why Negative-Binomial Estimation Matters

A log-linked count model describes how covariates relate to a count’s conditional mean. When counts vary more than a Poisson model allows, a negative-binomial model can represent that extra variation while retaining a convenient multiplicative interpretation of coefficients. Estimation involves both the regression effects and a dispersion parameter: the effects determine the conditional mean, while dispersion affects the spread of counts around that mean.

This tutorial uses the NB2 specification, whose variance is the conditional mean plus a term proportional to its square. The distinction matters in practice: at the same fitted mean, models with different fitted dispersion can predict different probabilities for zero, one, or larger counts; observation-specific dispersion would require an extension to this specification.

Definition: In an NB2 regression, the conditional mean for observation \(i\) is \(\mu_i=E(Y_i\mid x_i)=\exp(x_i'\beta)\), where \(Y_i\) is a nonnegative integer count, \(x_i\) is a vector of covariates, and \(\beta\) is a vector of coefficients. Conditional on the covariates, \(Y_i\) has a negative-binomial distribution with dispersion \(\phi>0\) and variance \(\operatorname{Var}(Y_i\mid x_i)=\mu_i+\phi\mu_i^2\). The observations are assumed conditionally independent, and the same \(\phi\) is used for all observations in this specification.

The conditional probability mass function is

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

Here \(\Gamma(\cdot)\) is the gamma function, which extends the factorial because \(\Gamma(m)=(m-1)!\) for positive integers \(m\). As \(\phi\) approaches zero, the negative-binomial distribution approaches a Poisson distribution with mean \(\mu_i\). Positive \(\phi\) permits variance greater than the mean.

Likelihood Estimation

For a sample of \(n\) conditionally independent observations, the likelihood is the product of their probability mass functions. The log-likelihood is more convenient for computation: it turns that product into a sum. Substituting \(\mu_i=\exp(x_i'\beta)\), the log-likelihood is

\ell(\beta,\phi)= \sum_{i=1}^{n}\left[ \log\Gamma(y_i+1/\phi)-\log\Gamma(1/\phi)-\log(y_i!) +y_i\log\left(\frac{\phi\mu_i}{1+\phi\mu_i}\right) -\frac{1}{\phi}\log(1+\phi\mu_i) \right].

The maximum-likelihood estimates \((\hat\beta,\hat\phi)\) are the parameter values that maximize this log-likelihood. In general, they are found numerically. The regression coefficients and dispersion are estimated jointly, rather than first fitting a Poisson model and then treating a dispersion adjustment as if it were part of the same likelihood fit.

One part of the optimization can be understood from the score for \(\beta\). For a single observation, differentiating its log-likelihood with respect to the linear predictor \(x_i'\beta\) gives a contribution proportional to \(y_i-\mu_i\), but downweighted according to the fitted mean and dispersion:

\frac{\partial\ell_i}{\partial\beta} =x_i\frac{y_i-\mu_i}{1+\phi\mu_i}.

To see the result, first differentiate with respect to \(\mu_i\). The terms that depend on \(\mu_i\) give

\frac{\partial\ell_i}{\partial\mu_i} =\frac{y_i-\mu_i}{\mu_i(1+\phi\mu_i)}.

Since \(\partial\mu_i/\partial\beta=\mu_i x_i\), the chain rule gives the stated score. In the full estimation, the optimizer also updates \(\phi\), so the fitted values must jointly maximize the likelihood rather than merely solve this coefficient score with dispersion held fixed.

Worked Example: Comparing Probabilities at the Same Mean

Consider a count with conditional mean \(\mu=2\). Compare two NB2 specifications, one with \(\phi=0.40\) and one with \(\phi=0.80\). The mean is held fixed so that the comparison isolates the role of dispersion.

For NB2, the probability of zero is \(P(Y=0)=(1+\phi\mu)^{-1/\phi}\). With \(\phi=0.40\),

P(Y=0)=(1+0.40(2))^{-1/0.40} =1.8^{-2.5}\approx 0.2300.

The probability of one is the zero probability multiplied by \(\phi\mu/(1+\phi\mu)\). Thus,

P(Y=1)=1.8^{-2.5}\left(\frac{0.40(2)}{1.8}\right) \approx 0.2300(0.4444)\approx 0.1022.

For \(\phi=0.80\), the zero probability is

P(Y=0)=(1+0.80(2))^{-1/0.80} =2.6^{-1.25}\approx 0.3028,

and the probability of one is

P(Y=1)=2.6^{-1.25}\left(\frac{0.80(2)}{2.6}\right) \approx 0.3028(0.6154)\approx 0.1863.

The higher-dispersion model has a higher probability of zero in this comparison. It also has a higher probability of one here. That direction for \(P(Y=1)\) is not a universal rule across all means and dispersion values; the full probability mass shifts as dispersion changes. The conditional mean remains 2 in both models, but the distribution around that mean differs.

Interpreting Coefficients and Rate Ratios

Because the conditional mean is log-linked, a one-unit increase in a covariate changes the log of the expected count by its coefficient, holding the other covariates fixed. Exponentiating a coefficient expresses that change on a multiplicative scale. For a continuous covariate, \(\exp(\beta_j)\) is the ratio of conditional means for a one-unit increase in that covariate. For a binary covariate, it is the ratio comparing the group coded 1 with the group coded 0, conditional on the other covariates.

Interpretation: In a log-linked count model, \(\exp(\beta_j)\) is the incidence-rate ratio (IRR) for a one-unit increase in \(x_j\), holding the other included covariates fixed. An IRR above 1 corresponds to a higher conditional mean; an IRR below 1 corresponds to a lower conditional mean. The coefficient describes a conditional mean association in the specified model, not by itself a causal effect.

If the model includes exposure, such as person-time at risk, the log of exposure can be included as an offset with coefficient fixed at 1. In that case the exponential of the remaining linear predictor represents a conditional event rate per unit of exposure. Without such an offset, the mean is an expected count for the observation period represented by the data.

Worked Example: Estimating a Treatment Incidence-Rate Ratio

Suppose a hypothetical study models the number of events per patient-month using a binary treatment indicator \(T_i\), where 1 denotes treatment and 0 denotes the comparison group. A negative-binomial likelihood fit returns \(\hat\beta_0=0.4055\), \(\hat\beta_T=0.6931\), and \(\hat\phi=0.40\). The treatment coefficient is a fitted estimate from the hypothetical likelihood optimization.

For a comparison-group patient-month, \(T=0\), so the fitted mean is

\hat\mu_0=\exp(0.4055)\approx 1.500 \quad\text{events per patient-month}.

For a treated patient-month, \(T=1\), giving

\hat\mu_1=\exp(0.4055+0.6931) =\exp(1.0986)\approx 3.000 \quad\text{events per patient-month}.

The estimated IRR is

\widehat{\operatorname{IRR}}=\exp(\hat\beta_T) =\exp(0.6931)\approx 2.000.

In this fitted model, the treatment group’s conditional mean event rate is estimated to be twice the comparison group’s rate, at the same values of other included covariates. This is a model-based association; the interpretation does not establish that treatment caused the difference.

The estimate of \(\phi\) also affects fitted probabilities. For a treated patient-month, \(\hat\mu=3.000\) and \(\hat\phi=0.40\), so the fitted probability of zero events is

\widehat{P}(Y=0\mid T=1) =(1+0.40(3))^{-1/0.40} =2.2^{-2.5}\approx 0.1393.

The fitted mean and IRR are determined by \(\hat\beta\), but the distribution’s probabilities depend on both \(\hat\beta\) and \(\hat\phi\).

Inference for Effects

A fitted coefficient is an estimate, not an exact population value. Under regularity conditions and a correctly specified likelihood, maximum-likelihood estimates have an approximately normal sampling distribution in sufficiently large samples. Their estimated covariance matrix is typically obtained from the inverse of the observed information matrix, which reflects curvature of the log-likelihood at the joint estimate. The uncertainty calculation should account for estimating dispersion along with the regression coefficients.

A Wald confidence interval for a coefficient is \(\hat\beta_j\pm z^*SE(\hat\beta_j)\), where \(z^*\) is the standard normal critical value for the chosen confidence level. Exponentiating both endpoints gives an interval for the IRR. This interval is often asymmetric around the estimated IRR because exponentiation is nonlinear.

Worked Example: A Confidence Interval for an IRR

Suppose a hypothetical negative-binomial fit estimates a treatment coefficient of \(\hat\beta_T=0.4055\), with \(SE(\hat\beta_T)=0.18\). For a two-sided 95% Wald interval, use \(z^*=1.96\). The coefficient-scale interval is

0.4055\pm1.96(0.18) =0.4055\pm0.3528 =[0.0527,\ 0.7583].

The estimated IRR is \(\exp(0.4055)\approx1.500\). Exponentiating the interval endpoints gives

[\exp(0.0527),\exp(0.7583)]\approx[1.054,\ 2.135].

The model-based 95% confidence interval for the treatment IRR is approximately 1.054 to 2.135. In repeated samples under the model and regularity conditions, intervals constructed this way would cover the true IRR about 95% of the time. This interval does not mean there is a 95% probability that the fixed parameter lies within these particular endpoints.

For testing \(H_0:\beta_T=0\) against a two-sided alternative, the Wald statistic is

z=\frac{0.4055}{0.18}\approx2.253.

The corresponding two-sided standard-normal p-value is approximately 0.02427. At significance level \(\alpha=0.05\), this is evidence against \(H_0\) under the fitted model. It is not the probability that \(H_0\) is true.

What Dispersion Changes

In NB2, dispersion affects more than a standard error. At a fixed conditional mean, it changes the shape of the fitted count distribution and therefore changes fitted probabilities such as \(P(Y=0)\). This is why \(\hat\phi\) must be considered when predicting the probability of particular counts, even if the main substantive interest is in the regression coefficients.

Dispersion also enters the likelihood and the estimated covariance matrix. A different fitted dispersion can therefore change coefficient estimates and their standard errors, because the likelihood weights observations differently and the coefficients and dispersion are optimized together. It is not generally correct to take a Poisson regression estimate and simply alter its reported standard errors to claim that an NB2 likelihood fit has been performed.

The conditional variance formula provides a direct check of the model’s implication. For example, if \(\mu=3\) and \(\phi=0.40\), then

\operatorname{Var}(Y\mid x)=3+0.40(3^2)=6.600,

which exceeds the mean of 3. This extra variability is built into the NB2 probability distribution, not added after the mean model is fitted.

Common Mistakes

Reading a coefficient as an additive change in counts. A log-linked coefficient changes the log conditional mean. Exponentiate it to obtain the multiplicative change in the mean.
Calling an IRR a causal effect. The IRR is a conditional model-based comparison. A causal interpretation requires an appropriate design and assumptions beyond fitting this regression.
Confusing the mean with the whole distribution. The same fitted mean can yield different probabilities for zero and other counts under different values of \(\phi\).
Treating dispersion as a standard-error-only correction. In NB2, dispersion is a parameter of the conditional distribution, affects the likelihood and fitted probabilities, and is estimated jointly with regression coefficients.
Interpreting a confidence interval as a probability statement about a fixed parameter. The confidence level describes the long-run coverage of the interval procedure under its assumptions.

Key takeaway: NB2 maximum likelihood estimates the log-linked conditional mean and the dispersion parameter together. Exponentiated regression coefficients are conditional incidence-rate ratios, while dispersion governs how probability is distributed around the fitted mean and contributes to likelihood-based inference.

Check Your Understanding

Answer these questions using the NB2 model and interpretations developed above.

Write the NB2 conditional variance in terms of \(\mu_i\) and \(\phi\). What happens to the variance as \(\phi\) increases while the mean is fixed?
A fitted coefficient for a one-unit increase in a covariate is \(-0.20\). What is the corresponding IRR, and how should it be interpreted?
Why can two NB2 models with the same fitted mean predict different probabilities of zero?
If a treatment coefficient is 0.30, what is the estimated treatment IRR? Does that estimate alone establish a causal effect?
Why is it inappropriate to describe an NB2 fit as a Poisson fit with only a different standard error?