Negative Binomial Regression for Overdispersed Count Data

Introduction

Clinical trials frequently generate count outcomes: number of exacerbations, hospitalizations, infections, rescue-medication uses, seizures, falls, bleeding events, or other recurrent events observed during a defined period. These outcomes are nonnegative integers, so ordinary linear regression is usually not the natural starting point.

A common first model is Poisson regression. The Poisson distribution has an important mathematical property: the conditional variance equals the conditional mean. In notation,

Real clinical count data often violate this assumption. Patients can differ substantially in their underlying susceptibility to events, follow-up can vary, and unmeasured risk factors can create more heterogeneity than the Poisson model allows. The result is often overdispersion: the observed variance is larger than the mean.

What Is Overdispersion?

Suppose a clinical trial records the number of COPD exacerbations experienced by each patient during follow-up. If the average patient has approximately 1.8 exacerbations but the observed variance is 5.9, then the variance is much larger than the mean.

The ratio of the variance to the mean is a useful descriptive diagnostic, although it is not by itself a formal test of the Poisson assumption. In regression, overdispersion should be assessed using the fitted model and the structure of the data rather than relying only on an overall sample mean and variance.

Why Overdispersion Matters

If a Poisson model is fitted to substantially overdispersed data, the estimated regression coefficients can still sometimes be useful as descriptions of the mean structure, but the usual Poisson standard errors may be too small. That can make confidence intervals artificially narrow and hypothesis tests overly optimistic.

The practical consequence is important: a treatment effect may appear more statistically convincing than the data actually support.

The Poisson Regression Model

Let \(Y_i\) denote the count outcome for patient \(i\), and let \(\mu_i=E(Y_i\mid X_i)\) be the conditional mean. Poisson regression uses a log link:

Equivalently,

For a binary treatment variable \(T_i\), the model may be written:

The exponentiated treatment coefficient is a multiplicative effect:

When the model represents an event rate through an exposure offset, this quantity is usually called an incidence rate ratio (IRR).

The Negative Binomial Model

The negative binomial model keeps the same log-linear mean structure but permits additional variation around the mean. A common parameterization is:

with variance:

Here, \(\alpha\) is a dispersion parameter. When \(\alpha=0\), the variance reduces to \(\mu_i\), recovering the Poisson variance relationship. When \(\alpha>0\), the variance is larger than the mean.

The Mean Model Is Still Log-Linear

The negative binomial regression mean is modeled as:

This is why negative binomial regression often feels familiar to anyone who has used Poisson regression. The key difference is not the link function; it is the variance model.

Where Does the Extra Variation Come From?

One useful way to understand the negative binomial distribution is as a mixture model. Suppose patients have different underlying event rates. Even if a Poisson distribution were appropriate for each patient's latent rate, heterogeneity across patients can produce extra-Poisson variation in the combined population.

This interpretation is especially intuitive in clinical trials. Two patients receiving the same treatment may have very different baseline susceptibility to recurrent events. A simple Poisson model treats the conditional mean as completely determining the variance, whereas the negative binomial model allows additional heterogeneity.

Adding an Exposure Offset

Count outcomes are often observed over different amounts of time. One patient may contribute 12 months of follow-up while another contributes only 6 months. Comparing raw counts without accounting for exposure time can be misleading.

For event rates, include the logarithm of exposure time as a fixed coefficient of one:

where \(t_i\) is the exposure or follow-up time. Equivalently,

The term \(\log(t_i)\) is called an offset. It is not estimated as an ordinary regression coefficient; its coefficient is fixed at 1.

Why the Offset Matters

Imagine two patients who experience two hospitalizations each. If Patient A is followed for one year and Patient B is followed for six months, the underlying hospitalization rates are not the same.

The offset allows the regression model to compare rates rather than treating equal counts as equal exposure.

Interpreting the Regression Coefficients

Suppose treatment is coded as \(T=1\) for the experimental arm and \(T=0\) for control:

The treatment incidence rate ratio is:

Interpretation depends on the direction of the outcome. For an undesirable event such as hospitalization, an IRR below 1 indicates a lower event rate in the treatment group.

For a coefficient \(\beta_1\), the approximate percent change in the expected rate is:

Confidence Intervals for the IRR

If the estimated treatment coefficient is \(\widehat{\beta}_1\) with standard error \(SE(\widehat{\beta}_1)\), an approximate 95% confidence interval on the log scale is:

Exponentiating the limits gives the IRR confidence interval:

If the 95% confidence interval excludes 1, the corresponding two-sided Wald test rejects the null hypothesis that the rate ratio equals 1 at the 5% level, subject to the usual model and large-sample approximations.

Worked Clinical Example: Recurrent Hospitalizations

Suppose a randomized Phase III trial compares an experimental treatment with placebo. The primary recurrent-event analysis considers the number of disease-related hospitalizations during follow-up. The investigators anticipate substantial heterogeneity between patients and therefore plan to evaluate overdispersion.

Example Results

Assume the fitted model produces:

The treatment effect is:

Thus, after accounting for the exposure time and other model terms, the estimated hospitalization rate in the experimental arm is approximately 25% lower than the rate in the placebo arm.

The 95% confidence interval of 0.60 to 0.94 does not include 1.00, indicating evidence of a difference in hospitalization rates under the fitted model.

Step-by-Step Model Construction

Poisson vs. Negative Binomial: Model Selection

A common workflow is to fit both Poisson and negative binomial models and assess whether the additional dispersion component is supported by the data and by clinical/statistical reasoning.

One descriptive quantity is the Pearson dispersion statistic:

where \(r_{P,i}\) are Pearson residuals. A value substantially greater than 1 can indicate overdispersion, although the exact interpretation depends on the model, sample size, covariates, and estimation method.

A Note About Quasi-Poisson Models

Another response to overdispersion is a quasi-Poisson model. It retains the Poisson mean structure but relaxes the variance assumption through a dispersion multiplier. A simplified representation is:

This can provide robust inference about the mean relationship when the main concern is the variance, but it differs from negative binomial regression in its variance structure and likelihood-based interpretation.

Zero Inflation: When Negative Binomial Is Not Enough

Some clinical count datasets contain many more zeros than a standard negative binomial model can reasonably accommodate. For example, a large proportion of patients may never experience a recurrent event because they belong to a structurally event-free subgroup, while the remaining patients have positive counts.

In that setting, a zero-inflated negative binomial model or hurdle model may be considered. These models answer somewhat different scientific questions and should not be selected merely because they produce a better numerical fit.

Repeated Events and Patient-Level Dependence

Negative binomial regression models the distribution of a count, but clinical recurrent events can have additional structure. Events within the same patient are not independent in a biological sense, and the timing of events may matter.

If the scientific question concerns the overall event rate over follow-up, a negative binomial rate model can be useful. If the analysis instead concerns event timing, within-patient event dependence, or recurrent-event processes, methods such as Andersen–Gill, Prentice–Williams–Peterson, frailty models, or other recurrent-event approaches may be more appropriate depending on the estimand.

Offset vs. Covariate: An Important Distinction

An exposure variable used as an offset is not the same as an ordinary covariate. For a follow-up time \(t_i\), the offset is:

Its coefficient is fixed at 1. If instead you include follow-up time as an ordinary covariate, the model estimates its coefficient from the data and makes a different assumption about how exposure relates to the mean count.

Example With a Continuous Covariate

Suppose age is included along with treatment:

The treatment IRR is still:

The age effect is:

which is the multiplicative change in the event rate associated with a one-unit increase in age, holding the other variables constant. For clinically interpretable reporting, a larger increment such as 10 years may be more useful:

Interactions

An interaction can be used when the treatment effect is expected to differ across a prespecified subgroup. For example:

The treatment IRR in the reference subgroup is \(\exp(\beta_1)\), whereas the treatment IRR in the \(G=1\) subgroup is:

Subgroup analyses should be prespecified when possible, and interaction terms should not be introduced simply because one subgroup happens to have a favorable point estimate.

Model Diagnostics

A negative binomial model should not be accepted solely because its estimated dispersion parameter is positive. Model adequacy should be evaluated using several complementary diagnostics.

1. Dispersion

Confirm that the fitted model provides a materially better representation of the observed variability than a Poisson model when such a comparison is appropriate.

2. Predicted vs. Observed Counts

Compare the distribution of observed counts with model-predicted counts. Look for systematic underprediction of high counts or other patterns that suggest misspecification.

3. Residuals

Examine Pearson, deviance, or other appropriate residuals. Patterns can indicate omitted covariates, nonlinear effects, incorrect functional forms, or influential observations.

4. Influential Patients

A small number of patients with unusually high event counts can have substantial influence on a count model. Investigate whether such observations are clinically valid and whether the fitted effect depends disproportionately on them.

5. Covariate Functional Form

Continuous predictors should not automatically be assumed to have a linear effect on the log mean. Restricted cubic splines, transformations, or clinically meaningful categories may sometimes be appropriate.

Likelihood Ratio Comparison With Poisson

If the negative binomial model is formulated as a nested extension of the Poisson model, a likelihood-based comparison can sometimes be used to assess whether the additional dispersion parameter is supported. However, the boundary nature of the null value and the exact estimation method require careful interpretation.

In practice, do not base model choice on a single automatic test. Consider the dispersion diagnostics, clinical structure, goodness of fit, prespecified analysis strategy, and scientific estimand together.

SAS Example

A common SAS implementation uses PROC GENMOD with a negative binomial distribution and log link. The following example models recurrent hospitalizations with treatment and age and uses follow-up time as an offset.

proc genmod data=analysis;
    class treatment(ref='Placebo');
    model hospitalizations =
          treatment age
          / dist=negbin
            link=log
            offset=log_followup
            type3;
    estimate 'Treatment IRR'
             treatment 1 -1
             / exp cl;
run;

Here, log_followup should contain the natural logarithm of the exposure time. The exact treatment contrast depends on how the treatment variable is coded.

R Example

In R, the MASS package provides a widely used glm.nb() function for fitting negative binomial generalized linear models.

library(MASS)

fit_nb <- glm.nb(
  hospitalizations ~ treatment + age + offset(log(followup)),
  data = analysis
)

summary(fit_nb)

exp(coef(fit_nb))
exp(confint(fit_nb))

The exponentiated treatment coefficient is the treatment rate ratio when the model is specified with the exposure offset.

Worked Calculation of an IRR

Suppose the fitted treatment coefficient is:

Then:

The estimated event rate is therefore approximately 30% lower in the treatment group:

If the 95% confidence interval for the IRR were 0.55 to 0.89, the appropriate interpretation would be that the treatment group has an estimated 30% lower event rate, with a 95% confidence interval corresponding to rate reductions of approximately 11% to 45%.

Clinical Interpretation: Count vs. Risk

One of the most common reporting errors is to confuse a rate ratio with a risk ratio. These are different estimands.

A patient can experience zero, one, or many events. Therefore, a recurrent-event count analysis should be reported using language consistent with the estimand.

Common Mistakes

1. Using Poisson regression automatically for every count outcome. A count outcome does not guarantee that the Poisson variance assumption is adequate.
2. Checking only the overall sample variance. Regression overdispersion is conditional on the model and covariates; assess the fitted model rather than relying on one descriptive ratio.
3. Ignoring follow-up time. When the estimand is an event rate, unequal exposure should generally be handled through an appropriate offset.
4. Interpreting an IRR as a risk ratio. An IRR concerns event rates, not necessarily the probability of at least one event.
5. Calling every excess zero “zero inflation.” Excess zeros require a stronger structural interpretation than ordinary overdispersion.
6. Ignoring recurrent-event dependence. A count model may not answer a question about event timing or within-patient event processes.
7. Comparing dispersion parameters across software without checking parameterization. Negative binomial parameter definitions differ across implementations.
8. Reporting only the p-value. For clinical interpretation, report the IRR, confidence interval, and clearly defined estimand.
9. Adding interactions after looking at the data. Subgroup effects and interactions should generally be scientifically motivated and prespecified where possible.
10. Assuming negative binomial regression solves every count-model problem. Zero inflation, truncation, dependence, structural zeros, and time-to-event questions may require other methods.

A Practical Analysis Workflow

What Should Be Reported in a Statistical Analysis Plan?

For a clinical trial using negative binomial regression, the SAP should make the analysis reproducible and clinically interpretable. At minimum, describe:

• Definition of the count endpoint
• Observation and exposure period
• Primary analysis population
• Treatment coding and reference arm
• Covariates and prespecified interactions
• Link function
• Negative binomial parameterization
• Offset definition, if applicable
• Method for estimating the dispersion parameter
• Primary effect measure, usually the IRR for a rate model
• Confidence interval and hypothesis-testing method
• Handling of incomplete follow-up
• Rules for missing or unavailable event information
• Model diagnostics and sensitivity analyses
• Criteria for considering alternative count or recurrent-event models

Sensitivity Analyses

Sensitivity analyses can be particularly useful when the count endpoint has unusual features. Potential analyses may include:

• Poisson regression with an appropriate variance adjustment
• Negative binomial regression with alternative clinically justified covariate specifications
• Alternative handling of follow-up or exposure time
• Models addressing influential observations
• Alternative recurrent-event approaches when event timing is important
• Zero-inflated or hurdle models when there is a strong scientific rationale for excess zeros
• Alternative definitions of the recurrent-event endpoint

Negative Binomial Regression in a Clinical Trial CSR

A concise CSR presentation might report the treatment effect as follows:

The final wording should reflect the actual prespecified model, estimation method, confidence interval method, and analysis population.

How to Decide Whether Negative Binomial Regression Fits the Question

The Most Important Concept

Negative binomial regression is not simply “Poisson regression with a different distribution.” The central modeling change is the variance assumption.

Poisson regression imposes:

Negative binomial regression allows:

That extra variance can be crucial when clinical count outcomes show substantial heterogeneity. The mean model remains log-linear, and when an exposure offset is included, the exponentiated treatment coefficient has a direct interpretation as an incidence rate ratio.

Summary

References

Hilbe, J.M. Negative Binomial Regression. Cambridge University Press.

McCullagh, P. & Nelder, J.A. Generalized Linear Models. Chapman & Hall.

Agresti, A. Foundations of Linear and Generalized Linear Models. Wiley.

Cameron, A.C. & Trivedi, P.K. Regression Analysis of Count Data. Cambridge University Press.

Long, J.S. & Freese, J. Regression Models for Categorical Dependent Variables Using Stata. Stata Press.