Negative Binomial Regression

• Negative binomial regression is a flexible framework for modeling overdispersed count data, extending Poisson models with an explicit dispersion parameter.
• The method improves inference accuracy by addressing variance (overdispersion) through both classical and Bayesian estimation techniques.
• Advanced extensions include penalized likelihood, hierarchical mixtures, and zero-inflated models to effectively handle high-dimensional and structured data.

Negative binomial regression is a flexible modeling framework for count data exhibiting overdispersion, where the variance of the outcome exceeds its mean. It generalizes Poisson regression by introducing an explicit dispersion parameter, enabling accurate inference and prediction in diverse empirical contexts. The methodology is grounded in generalized linear modeling, but supports both frequentist and Bayesian inference, penalized high-dimensional extensions, model averaging, and mixture/hierarchical generalizations. This article provides a comprehensive exposition of negative binomial regression, from model specification to advanced inferential and computational techniques.

1. Model Specification and Likelihood

Negative binomial regression models the conditional distribution of non-negative count outcomes $Y_i$ using the negative binomial (NB) family, parameterized by a mean $\mu_i>0$ and a positive dispersion (or shape) parameter $\alpha$, $\theta=1/\alpha$, or $r$ in alternative notations. The canonical log-link is standard. For each observation $i$:

$$n_i \mid \mu_i, \alpha \sim \mathrm{NegBin}(\mu_i, \alpha)$$

with probability mass function

$$p(n_i \mid \mu_i, \alpha) = \frac{\Gamma(n_i+\alpha)}{n_i!\,\Gamma(\alpha)} \left(\frac{\alpha}{\alpha+\mu_i}\right)^{\alpha} \left(\frac{\mu_i}{\alpha+\mu_i}\right)^{n_i}$$

and

$$\mathbb{E}[n_i] = \mu_i, \qquad \mathrm{Var}(n_i) = \mu_i + \frac{\mu_i^2}{\alpha}$$

(Cabansag et al., 4 Jan 2026, Settanni, 2024). Overdispersion, distinct from the Poisson ($\mathrm{Var}=\mathbb{E}$) case, is controlled by $\mu_i>0$0.

Predictors $\mu_i>0$1 are linked to the mean via

$\mu_i>0$2

With a vector of observations, maximum-likelihood or Bayesian estimation proceeds via the joint log-likelihood or posterior.

2. Inference: Classical, Bayesian, and Bias-Reduction

Classical Inference

MLE for $\mu_i>0$3 is performed by maximizing the likelihood or solving the corresponding score equations: $\mu_i>0$4

$\mu_i>0$5

(Settanni, 2024). Standard Newton-Raphson or iteratively reweighted least squares is used.

Bias in the MLE of $\mu_i>0$6 (downward, $\mu_i>0$7) is acute in small samples and propagates to $\mu_i>0$8 estimation (Pagui et al., 2020). Adjusted score methods, such as mean- and median-bias reduction (Firth and Kenne Pagui et al.), incorporate $\mu_i>0$9 corrections: $\alpha$0 yielding estimators with substantially improved finite-sample properties, especially for coverage of interval estimates on dispersion (Pagui et al., 2020).

Bayesian Inference

A fully Bayesian NB regression is specified by priors: $\alpha$1 (Cabansag et al., 4 Jan 2026). Posterior inference is performed via Markov chain Monte Carlo, e.g., with NUTS or JAGS (Souza et al., 2015). Outputs include posterior means, medians, credible intervals, as well as functions such as incidence rate ratios $\alpha$2. Posterior predictive checks are essential to assess model adequacy, especially for the heavy tails and skewness characteristic of overdispersed counts.

In high-dimensional and sparse contexts, Bayesian hierarchical regularization (e.g. via gamma priors on dispersion, lognormal random effects, or in the zero-inflated regime) provides adaptive shrinkage, stability, and explicit uncertainty quantification (Zhou et al., 2012, Jiang et al., 2018).

3. Penalization, Model Averaging, and High-Dimensional NB Regression

Penalized Likelihood

For $\alpha$3, $\alpha$4-penalized NB regression (lasso) and elastic-net regularization are provably consistent under standard restricted eigenvalue and compatibility conditions (Xie et al., 2020, Zhang et al., 2017, Zilberman et al., 2024): $\alpha$5 Support recovery and oracle inequalities hold under sparsity and design conditions.

Complexity-penalized MLE (with $\alpha$6) is adaptively minimax across all sparsity regimes, and practical convex surrogates (lasso, SLOPE) achieve optimal rates for high-dimensional count GLMs (Zilberman et al., 2024).

Weighted-average least squares (WALS) extends model averaging to negative binomial GLMs by combining submodel estimates with Bayesian-style shrinkage (semi-orthogonal transformation, no cross-validation), yielding improved prediction and stability compared to MLE or lasso in situations with auxiliary covariate uncertainty (Huynh, 2024).

Liu-Type Estimators

When regressors are highly collinear, the Liu-type negative binomial estimator introduces two parameters $\alpha$7 to shrink the IRLS updates: $\alpha$8 This estimator strictly dominates MLE, ridge, or the classical Liu estimator under explicit matrix mean squared error spectral conditions, with substantial improvements in MSE under high multicollinearity (Asar, 2016).

4. Extensions: Mixtures, Zero-Inflation, and Robust Variants

Mixtures and Inflation

Mixtures of NB components are used when data show both overdispersion and multi-modality or excess mass at a particular count (e.g., zero or one). The $\alpha$9-inflated negative binomial mixture model augments a finite NB mixture with an additional atom at count $\theta=1/\alpha$0: $\theta=1/\alpha$1 EM algorithms provide efficient estimation (Najafabadi et al., 2017).

Zero-inflated NB (ZINB) regression decomposes the data into a structural zero component and a conditional NB process, with hierarchical Bayesian implementations facilitating integration of covariate and group effects, feature selection, and high-dimensional settings (e.g. microbiome-omics) (Jiang et al., 2018).

Multivariate Negative Binomial Models

Hierarchical Poisson-log-gamma models marginalize to multivariate NB distributions, supporting modeling of overdispersed and correlated multivariate counts with closed-form likelihoods. Randomized quantile residuals, global and local influence diagnostics, and dedicated R implementations extend standard GLM diagnostics to this correlated count domain (Fabio et al., 2021).

Robust and Semiparametric Regression

Robust negative binomial M-quantile regression solves estimating equations with bounded influence $\theta=1/\alpha$2-functions. This approach is less sensitive to outliers than likelihood-based NB regression, can flexibly model spatial heterogeneity (via smoothed M-quantile indices), and empirically yields lower RMSE in disease mapping and spatial count prediction when compared to Empirical Bayes and Besag-York-Mollié models (Chambers et al., 2013).

5. Practical Considerations and Diagnostics

Key recommendations include:

• Always test for overdispersion; if $\theta=1/\alpha$3, NB is preferred over Poisson (Wanga et al., 2022, Cabansag et al., 4 Jan 2026).
• For small or moderate sample sizes, use bias-reduced or Bayesian estimators rather than unadjusted MLE for the dispersion parameter (Pagui et al., 2020).
• Posterior predictive checks (Bayesian) and standardized residuals (frequentist) are critical for model adequacy, especially in the upper tail and for zero-inflation.
• Model selection may be performed via AIC, BIC, deviance, or more advanced criteria (DIC, predictive risk), but penalized or averaged estimators often yield superior prediction and interpretability, especially in high-dimensional/rich covariate settings (Zilberman et al., 2024, Huynh, 2024).
• For very large-scale or constrained computational scenarios, method of moments or transformer-based direct estimators offer dramatic speedups with comparable accuracy to likelihood-based approaches in simple designs (Svensson, 6 Aug 2025).

6. Applications and Interpretative Quantities

Negative binomial regression has been deployed in diverse empirical contexts: chart persistence of popular music (Cabansag et al., 4 Jan 2026), urban trail traffic (Wanga et al., 2022), disease mapping and spatial epidemiology (Chambers et al., 2013), insurance rate setting (Najafabadi et al., 2017), galaxy cluster counts (Souza et al., 2015), and omics/metagenomics (Jiang et al., 2018).

Parametric coefficients are typically interpreted via incidence rate ratios $\theta=1/\alpha$4, corresponding to multiplicative effects on the mean count per unit change in covariate $\theta=1/\alpha$5. When using the log-link, these ratios are immediately interpretable as relative changes in expected counts (Cabansag et al., 4 Jan 2026). For mixture and zero-inflated models, interpretation proceeds via predictive marginals and Bayes risk estimates (Najafabadi et al., 2017).

Posterior predictive intervals, probability of direction (e.g., $\theta=1/\alpha$6), and model comparison diagnostics formalize inference, model fit, and predictive assessment. Applications frequently validate models via cross-validation, out-of-sample prediction error (e.g. MAPE), and, in Bayesian contexts, posterior credible intervals and DIC/PPC (Wanga et al., 2022, Souza et al., 2015).

7. Methodological and Computational Issues

A summary of methodological recommendations:

• Analytical fidelity is critical in deriving derivatives, Fisher information, and inferential statistics; omitting or mishandling log-gamma terms may introduce inconsistent or biased standard errors for $\theta=1/\alpha$7 (Settanni, 2024).
• For high throughput or fixed-design screening, method-of-moments estimators are extremely efficient and have well-calibrated $\theta=1/\alpha$8-values, outperforming both MLE and novel transformer-based methods in accuracy and speed (Svensson, 6 Aug 2025).
• For robust inference under collinearity or high-dimensionality, Liu-type, lasso, elastic-net, and SLOPE penalization are theoretically justified and computationally efficient (Asar, 2016, Zhang et al., 2017, Zilberman et al., 2024).
• Diagnostic tools including randomized quantile residuals, Cook's distance (case-deletion), and local influence analysis are well developed for both univariate and multivariate NB regression and should be employed to detect misfit, leverage, and influential points (Fabio et al., 2021).
• Careful choice and consistent use of the parameterization for the dispersion/shape parameter is necessary for correct analytic derivatives and valid inference (Settanni, 2024).

In summary, negative binomial regression provides a robust, general modeling strategy for overdispersed counts, with well-established theory, practical algorithms, and advanced generalizations for high-dimensionality, mixture structure, robust estimation, and latent-variable modeling. Its ongoing refinement continues to broaden its application in empirical science and advanced statistical methodology.