Negative Binomial Regression

Negative Binomial Regression for Health Care Use and Event Rates

Negative binomial regression is a regression model for counts of health care events, such as GP consultations, hospital admissions or exacerbations of a chronic disease, in which some patients have far more events than others. It keeps the log link and rate ratios of Poisson regression but adds one parameter that lets the variance of the count grow faster than its mean. In health economics it supplies the event rates and rate ratios that drive resource use and costs in trial-based and model-based cost-effectiveness analysis. The extra parameter matters because a Poisson fit to the same data usually reports standard errors that are too small, so treatment effects look more certain than they are. It is a generalised linear model fitted by maximum likelihood estimation.

Why health care counts outgrow the Poisson model

A Poisson model assumes that the variance of a count equals its mean, which holds when every patient with the same covariates shares one underlying event rate. Real patients differ in severity, frailty and habits that no dataset records in full, so the observed counts spread out more than that. Gardner, Mulvey and Shaw describe the Poisson model as likely to mislead unless restrictive assumptions are met, because individual counts are usually more variable, or overdispersed, than the model implies.

The Stata reference manual explains that overdispersion comes about when some of the causes of the Poisson process are unknown. It describes underdispersion, a variance below the mean, as uncommon, and the negative binomial model does not help with it.

Cameron and Trivedi, in a 1986 paper on count data models, applied their estimators and tests to the number of doctor consultations. Jones, reviewing models for health care, notes that the Poisson model and extensions such as the negative binomial are often used in health economics to model the number of visits to a doctor.

The model as a Poisson process with a gamma-distributed frailty

The usual derivation gives each patient a personal Poisson rate. The mean rate depends on covariates through a log link, and an unobserved multiplier with mean one scatters individual rates around it.

$$y_i \mid \nu_i \sim \text{Poisson}(\nu_i \mu_i), \qquad \ln \mu_i = \beta_0 + \beta_1 x_{1i} + \dots + \beta_k x_{ki} + \ln t_i, \qquad \nu_i \sim \text{Gamma}(1/\alpha,\ \alpha)$$

where $y_i$ is the observed count for patient $i$, $\mu_i$ is the expected count, $x_{1i}$ to $x_{ki}$ are the covariates, $\beta_0$ to $\beta_k$ are the regression coefficients, $t_i$ is the follow-up time entering as an offset, $\nu_i$ is the unobserved multiplier, gamma-distributed with shape $\tfrac{1}{\alpha}$ and scale $\alpha$ so that its mean is one, and $\alpha$ is the overdispersion parameter, which is the variance of $\nu_i$.

The Stata manual shows that this gamma mixture of Poisson variables gives a negative binomial distribution for the count.

$$\Pr(y_i) = \frac{\Gamma(m + y_i)}{\Gamma(y_i + 1)\,\Gamma(m)}\, p_i^{\,m} (1 - p_i)^{y_i}, \qquad p_i = \frac{1}{1 + \alpha \mu_i}, \qquad m = \frac{1}{\alpha}$$

where $\Gamma$ is the gamma function, $p_i$ is a probability that falls as the expected count $\mu_i$ rises, and $m$ is the inverse of the overdispersion parameter $\alpha$.

Two ways of letting the variance grow

The mean and variance follow directly from the mixture. In the form above, known as NB2 in the terminology of Cameron and Trivedi and as the mean-dispersion model in Stata, the variance rises with the square of the mean.

$$\mathrm{E}(y_i) = \mu_i, \qquad \mathrm{Var}(y_i) = \mu_i (1 + \alpha \mu_i)$$

where $\mu_i$ is the expected count and $\alpha$ is the overdispersion parameter. The Poisson model is the special case $\alpha = 0$, and the larger $\alpha$ is, the greater the overdispersion.

The alternative NB1, or constant-dispersion, model sets the variance to $\mu_i (1 + \delta)$, a fixed multiple of the mean for every patient, where $\delta$ is its dispersion parameter. Some texts write NB2 with $r = 1/\alpha$, so the variance appears as $\mu_i + \mu_i^2 / r$. A reported dispersion value therefore cannot be read without knowing which notation was used. Hilbe's book on the model gives alternative variance parameterisations a chapter of their own.

Estimation, exposure and the test for overdispersion

The coefficients and the dispersion parameter are estimated together by maximising the log-likelihood. Stata searches over $\ln \alpha$ because $\alpha$ must be greater than zero.

Follow-up often differs between patients. The log of each patient's time at risk enters the linear predictor with its coefficient fixed at one, which the Stata manual calls an exposure variable. The coefficients then describe event rates per unit of time, not counts per patient.

Whether a negative binomial model is needed can be tested with a likelihood ratio test of $\alpha = 0$ against the Poisson fit. The manual explains that this test sits on the boundary of the parameter space, so the usual normal-theory test is omitted and the null distribution of the likelihood ratio statistic is an equal mixture of a point mass at zero and a chi-squared distribution with one degree of freedom. The Akaike information criterion is often used alongside it when comparing count models that are not nested.

Reading the coefficients as rate ratios

Because of the log link, each coefficient acts multiplicatively on the expected rate. Exponentiating a coefficient gives an incidence rate ratio, the ratio of expected event rates for a one-unit change in that covariate with the others held fixed.

$$\mathrm{IRR}_j = e^{\beta_j}$$

where $\mathrm{IRR}_j$ is the incidence rate ratio for covariate $j$ and $\beta_j$ is its coefficient. A treatment coefficient of minus 0.223 gives a rate ratio of 0.80, a 20% lower event rate at any level of the other covariates. The dispersion parameter adds what the rate ratio leaves out, how unevenly events fall across patients.

Worked example with exacerbation counts in a two-arm trial

All numbers in this example are illustrative. A trial randomises 200 patients with chronic lung disease to usual care and 200 to a self-management programme, and counts moderate or severe exacerbations over follow-up. An NB2 model with a log follow-up offset is fitted.

$$\ln \mu_i = 0.405 - 0.223\, T_i + \ln t_i, \qquad \hat{\alpha} = 0.8$$

where $\mu_i$ is the expected number of exacerbations for patient $i$, $T_i$ is 1 for the programme and 0 for usual care, $t_i$ is follow-up in years and $\hat{\alpha}$ is the estimated overdispersion parameter.

1. Convert the coefficients. The usual-care rate is $e^{0.405} = 1.50$ exacerbations per patient-year. The rate ratio is $e^{-0.223} = 0.80$, so the programme rate is 1.50 × 0.80 = 1.20 per patient-year.

2. Use the offset. A usual-care patient who withdrew after six months has $\ln t_i = \ln 0.5 = -0.693$, so the expected count is 1.50 × 0.5 = 0.75.

3. Compare the variances. Over one year, the usual-care factor is 1 + 0.8 × 1.5 = 2.2, so NB2 gives a variance of 1.5 × 2.2 = 3.30, against 1.50 under a Poisson model. For the programme arm the factor is 1 + 0.8 × 1.2 = 1.96 and the variance is 1.2 × 1.96 = 2.352, against 1.20.

4. Compare the share with no exacerbations. Under NB2 the probability of a zero count is $(1 + \alpha \mu)^{-1/\alpha}$. For usual care this is $(2.2)^{-1.25} = 0.373$, against $e^{-1.5} = 0.223$ under Poisson. The negative binomial predicts that about 37% of usual-care patients have a year free of exacerbations.

5. Compare the precision of the rate ratio. With equal one-year follow-up the variance of the log rate in each arm is approximately $\frac{1}{n\mu} + \frac{\alpha}{n}$, where $n$ is the number of patients in the arm. For usual care the two terms are 1/300 = 0.00333 and 0.8/200 = 0.004, which sum to 0.00733. For the programme they are 1/240 = 0.00417 and 0.004, which sum to 0.00817. The variance of the log rate ratio is their sum, 0.0155, so the standard error is 0.1245. A Poisson analysis drops the $\alpha/n$ terms, giving 0.00333 + 0.00417 = 0.0075 and a standard error of 0.0866.

6. Compare the intervals. The negative binomial 95% interval for the log rate ratio is minus 0.223 plus or minus 1.96 × 0.1245 = 0.244, or minus 0.467 to 0.021, which is a rate ratio of 0.63 to 1.02. The Poisson interval is minus 0.223 plus or minus 1.96 × 0.0866 = 0.170, or minus 0.393 to minus 0.053, a rate ratio of 0.68 to 0.95. The same estimate looks conclusive under Poisson and uncertain once patients are allowed to differ.

7. Turn rates into costs. At an illustrative GBP 1,200 per exacerbation, expected costs are 1.50 × 1,200 = GBP 1,800 per patient-year under usual care and 1.20 × 1,200 = GBP 1,440 on the programme, a saving of GBP 360. The mean saving is the same under either model, but the standard error of the log rate ratio behind it is 0.1245 / 0.0866 = 1.44 times as large under the negative binomial.

Using negative binomial results in cost-effectiveness work

In an economic evaluation the expected rate per arm, multiplied by unit costs and by the QALY loss per event, gives the mean costs and outcomes. The rate ratio and its standard error typically feed a decision model, where the log rate ratio is drawn from a normal distribution in probabilistic sensitivity analysis. A Poisson standard error used at that point understates decision uncertainty and the value of further research.

Respiratory trials show why the choice matters. Keene and colleagues, analysing exacerbations in the TRISTAN trial of 1,465 patients with chronic obstructive pulmonary disease, found that the key assumptions of Poisson regression often appear unreasonable for exacerbation data. They favoured a negative binomial model, which corresponds to assuming a separate Poisson parameter for each patient. Suissa, using a cohort of 5,454 patients with the disease arranged to emulate a trial, showed that ignoring between-patient variability in exacerbation rates turned a rate ratio of 0.79 with a p-value of 0.12 into one with a p-value of 0.0007.

Choosing between negative binomial and neighbouring count models

The negative binomial is one of several answers to overdispersion. Each alternative changes a different part of the model.

1. Poisson with robust standard errors. Jones notes that the Poisson estimator is consistent so long as the mean is correctly specified, even if the variance is wrong, with sandwich estimator standard errors replacing the model-based ones. This gives valid rate ratios but no credible predicted distribution of counts per patient.
2. Overdispersed or quasi-Poisson models. Gardner and colleagues describe a scale factor that corrects the inferential statistics instead of adding a random term for each patient.
3. Zero-inflated and hurdle models. Mihaylova and colleagues describe two-part models with a truncated negative binomial for users, and report that zero-inflated Poisson and negative binomial models can fail to converge if the same covariates are used in both parts.
4. Finite mixtures. Deb and Trivedi fitted a negative binomial finite mixture model to six measures of medical care demand by older people and found it preferred to the standard and hurdle versions, with two latent groups described as healthy and ill.
5. Repeated counts per patient. Counts recorded over several periods for the same patient are correlated, which generalised estimating equations or a mixed-effects model can handle.
6. Time-to-event analysis. Keene and colleagues note that clinical interest in a direct comparison of rates favours a count-based analysis over multiple time-to-event methods.

Common misreadings and limits

Most errors come from treating the dispersion parameter as a nuisance or from reading coefficients on the wrong scale. A few are specific to count data.

1. Reading a coefficient as a count. A coefficient of minus 0.223 is a 20% lower rate, and the number of events avoided depends on the baseline rate and follow-up.
2. Ignoring follow-up time. Without an exposure offset, patients who leave early look like low users and the rate ratio is biased whenever follow-up differs between the arms.
3. Reporting a dispersion value without its parameterisation. An $\alpha$ of 0.8 under NB2 is $r = 1.25$ in the inverse notation.
4. Testing overdispersion against the ordinary chi-squared. The boundary problem makes that test conservative.
5. Treating the model as a cure for every zero or underdispersion. Structural zeros may still need a two-part model, and a variance below the mean is outside its scope.

Sources

• Cameron AC, Trivedi PK. Econometric models based on count data: comparisons and applications of some estimators and tests. Journal of Applied Econometrics. 1986;1(1):29-53. https://ideas.repec.org/a/jae/japmet/v1y1986i1p29-53.html
• Deb P, Trivedi PK. Demand for medical care by the elderly: a finite mixture approach. Journal of Applied Econometrics. 1997;12(3):313-336. https://ideas.repec.org/a/jae/japmet/v12y1997i3p313-36.html
• Gardner W, Mulvey EP, Shaw EC. Regression analyses of counts and rates: Poisson, overdispersed Poisson, and negative binomial models. Psychological Bulletin. 1995;118(3):392-404. https://doi.org/10.1037/0033-2909.118.3.392
• Hilbe JM. Negative Binomial Regression. 2nd ed. Cambridge University Press; 2011. https://www.cambridge.org/core/books/negative-binomial-regression/12D6281A46B9A980DC6021080C9419E7
• Jones AM. Models for health care. HEDG Working Paper 10/01. University of York; January 2010. https://www.york.ac.uk/media/economics/documents/herc/wp/10_01.pdf
• Keene ON, Jones MR, Lane PW, Anderson J. Analysis of exacerbation rates in asthma and chronic obstructive pulmonary disease: example from the TRISTAN study. Pharmaceutical Statistics. 2007;6(2):89-97. https://doi.org/10.1002/pst.250
• Mihaylova B, Briggs A, O'Hagan A, Thompson SG. Review of statistical methods for analysing healthcare resources and costs. Health Economics. 2011;20(8):897-916. https://doi.org/10.1002/hec.1653
• StataCorp. nbreg, negative binomial regression. Stata Base Reference Manual. College Station, TX: StataCorp LLC; 2025. https://www.stata.com/manuals/rnbreg.pdf
• Suissa S. Statistical treatment of exacerbations in therapeutic trials of chronic obstructive pulmonary disease. American Journal of Respiratory and Critical Care Medicine. 2006;173(8):842-846. https://doi.org/10.1164/rccm.200508-1338PP