Negative Binomial Distribution for Overdispersed Counts

What the negative binomial distribution is

The negative binomial (NB) distribution is a discrete probability distribution on the non-negative integers {0, 1, 2, …} characterized by two parameters: a mean μ > 0 and a dispersion parameter α ≥ 0. Under the canonical NB-2 (quadratic) parameterization — the default in R's MASS::glm.nb, Stata's nbreg, and Python's statsmodels — the variance function is:

Var(Y) = μ + α μ²

When α = 0 the distribution collapses exactly to a Poisson with variance equal to the mean. When α > 0 the variance exceeds the mean by α μ², and the distribution develops a heavier right tail than Poisson: simultaneously more probability mass on zero counts and on large counts. This dual excess — more zeros and more large values than Poisson — is precisely the signature of healthcare utilization data. The probability mass function for the NB-2 family is:

P(Y = k) = Γ(k + 1/α) / [Γ(1/α) k!] × (αμ)^k / (1 + αμ)^(k + 1/α)

where Γ(·) is the gamma function and k = 0, 1, 2, …. For applied work, the two properties that matter are E(Y) = μ and Var(Y) = μ + αμ². Everything else — the PMF, likelihood, moment generating function — is machinery for understanding why software does what it does.

The Poisson-gamma mixture: why the NB arises in heterogeneous populations

The single most important conceptual result for RWE analysts is that the negative binomial distribution is not a convenient mathematical choice — it is the natural statistical consequence of patient heterogeneity. Suppose each patient has their own underlying event rate λᵢ, and that their counts conditional on that rate follow Poisson(λᵢ).

If the latent rates λᵢ vary across the population according to a gamma distribution with mean μ and shape parameter r = 1/α, then the marginal distribution of counts — the distribution you observe when you pool across all heterogeneous patients — is exactly the negative binomial with mean μ and dispersion α = 1/r. This result is sometimes called the Poisson-gamma mixture.

In pharmacoepidemiology terms: even if every individual patient's count process were exactly Poisson (memoryless, stationary events over time), the observed counts across a claims cohort will follow a negative binomial because patients differ in unmeasured ways. This unmeasured heterogeneity — called frailty in the survival analysis literature — integrates out to yield the NB marginal distribution. The parameter α directly measures the degree of population heterogeneity: larger α means greater between-patient variation in event rates.

When α is estimated from data, it is not arbitrary: it quantifies how much of the count variation is due to patient-level differences in underlying propensity that the measured covariates do not capture.

Why overdispersion is the norm in healthcare count data

In any claims or EHR dataset, observed counts of hospitalizations, emergency department (ED) visits, outpatient encounters, or prescription fills are virtually always overdispersed relative to Poisson for interconnected reasons:

• Case-mix heterogeneity: the most severely ill patients in any condition-defined cohort account for a disproportionate share of utilization, creating a long right tail. A patient with decompensated heart failure accrues far more admissions than a stable patient on the same therapy — this between-patient spread drives variance well above the mean.
• Coding and billing variation: administrative data capture differs by health system, payer, region, and care setting. Practices that unbundle care into separate claims inflate counts for the same underlying service intensity, adding non-clinical variance.
• Clustering of care: patients who have one hospitalization are more likely to have another (rehospitalization, exacerbation, index event followed by complication). This positive within-patient autocorrelation inflates the between-patient variance in total counts beyond what a Poisson with a constant rate predicts.
• Structural access barriers: a subgroup with near-zero utilization (healthy, uninsured intermittently, or disengaged from care) coexists with high-utilizers, widening the count distribution beyond the Poisson envelope in both tails simultaneously.

How to recognize and diagnose overdispersion

Before fitting an NB model it is good practice to confirm overdispersion. The standard diagnostic tools are:

1. Mean-variance plot: for the raw data, plot the sample mean against the sample variance within subgroups or arms. Under Poisson, all points should fall on the line Var = mean. Points above this line indicate overdispersion; how far above gives a sense of α.
2. Pearson dispersion statistic: after fitting a Poisson model, compute Pearson χ² divided by residual degrees of freedom. A ratio substantially greater than 1 (rule of thumb: > 1.5 to 2.0) signals overdispersion. This quantity estimates the dispersion factor φ = 1 + αμ̄, where μ̄ is the average fitted mean.
3. Likelihood-ratio boundary test: formally test H₀: α = 0 (Poisson) against the NB alternative. Because α = 0 is on the boundary of the parameter space, the standard likelihood-ratio chi-square statistic has a non-standard null distribution; its p-value must be halved relative to a chi-square(1) distribution. In any typical HCRU dataset with n > 200 patients, this test almost always rejects the Poisson — overdispersion is the rule, not the exception.

The practical implication of ignoring overdispersion is severe: fitting a Poisson model to overdispersed data produces standard errors that are too small (anticonservative), confidence intervals that are too narrow, and p-values that are too small. Effects that are not statistically significant are declared significant. This is not a minor issue — it is a fundamental validity failure.

The Pearson χ²/df from the Poisson fit provides an empirical correction factor (quasi-Poisson scaling), but NB regression is preferred because it estimates α from the data via maximum likelihood and provides a proper likelihood for AIC/BIC model comparison.

NB1 versus NB2 variance functions

Two parameterizations of the negative binomial are in common use, distinguished by how the variance relates to the mean:

• NB-2 (quadratic): Var(Y) = μ + αμ². The variance grows quadratically with the mean, meaning that at high count levels the excess variance grows rapidly. This is the canonical parameterization in R (MASS::glm.nb, with θ = 1/α), Stata (nbreg), and statsmodels in Python. NB-2 arises directly from the Poisson-gamma mixture and is the theoretically motivated choice for heterogeneous populations.
• NB-1 (linear): Var(Y) = μ(1 + δ) = μ + δμ. The variance is proportional to the mean with a constant proportionality factor (1 + δ). Equivalent to a quasi-Poisson with overdispersion factor φ = 1 + δ. NB-1 is less common in HEOR practice; it implies that overdispersion is uniform across count levels, which is less realistic for HCRU where high-utilizers generate disproportionate variance.

For most HCRU applications — hospitalizations, ED visits, infusion counts — where a minority of frail high-utilizers dominates the right tail, NB-2 is the better fit and the safe default. NB-1 (or quasi-Poisson) is occasionally preferred for mild, uniform overdispersion when the NB-2 likelihood fails to converge.

Offsets for unequal observation time

When patients have different lengths of follow-up — because of disenrollment, death, or administrative censoring — raw counts cannot be compared directly between patients or arms. A patient followed for six months with two hospitalizations is not equivalent to a patient followed for twelve months with two hospitalizations. The solution is an offset: including log(person-time) as a term in the linear predictor with its coefficient fixed at 1, converting the count model into a rate model:

log(μᵢ) = β₀ + β₁X₁ + … + βₚXₚ + log(person-timeᵢ)

The exponentiated coefficient exp(βⱼ) is then an incidence rate ratio (IRR) — events per unit of person-time — rather than a raw count ratio. This is non-negotiable in claims-based HCRU studies where enrollment gaps, disenrollment, and death create differential follow-up across treatment arms. Omitting the offset silently converts the IRR into a raw count ratio that confounds the rate effect with the length of observation, biasing the estimate in whichever arm has shorter mean follow-up.

The regression application of this distributional primitive is described in detail in the companion model concept (poisson-negative-binomial-count-models).

Zero-inflation versus overdispersion: not the same problem

Overdispersion and zero-inflation are distinct phenomena that require different solutions. The negative binomial distribution naturally accommodates some excess zeros relative to Poisson, because patients with very low event rates will have many zero-count observations by chance. This is ordinary overdispersion, and NB handles it correctly. Zero-inflation, by contrast, refers to a structural subgroup that cannot possibly generate the event — patients who have permanently left the health system, who have a physiological barrier to the event, or whose zero counts reflect data absence rather than true non-occurrence.

When excess zeros genuinely reflect a "never-utilizer" structural process beyond what the NB distribution predicts (diagnosed by comparing predicted vs. observed zero probabilities from the fitted NB), a zero-inflated NB (ZINB) or hurdle model is appropriate. The hurdle model separates a logistic process (any events at all?) from a zero-truncated NB process (how many, given at least one?). Applying ZINB to ordinary overdispersion risks overfitting; applying plain NB to genuine structural zeros understates the zero probability and overstates the event rate among the non-zero subgroup.

Interpreting the output

Consider a negative binomial regression of annual hospitalization counts on a treatment indicator, controlling for age, sex, and a comorbidity index, with log(person-years) as the offset. The fitted model returns: treatment coefficient β = -0.223, 95% CI [-0.40, -0.05], dispersion parameter α = 0.8.

Formal interpretation: exp(-0.223) = 0.80 is the adjusted incidence rate ratio for the treatment group versus the control group. The expected rate of hospitalizations per person-year among treated patients is 0.80 times the rate among control patients — a 20% lower rate — holding age, sex, and comorbidity index fixed at their observed values. The 95% confidence interval on the rate-ratio scale is [exp(-0.40), exp(-0.05)] = [0.67, 0.95], which excludes 1.0 and is consistent with a statistically significant reduction at the 5% level.

The dispersion parameter α = 0.8 characterizes the variance function: Var(Y) = μ + 0.8 × μ². At a typical mean count of μ = 1.5 hospitalizations per year, this gives Var(Y) = 1.5 + 0.8 × 2.25 = 1.5 + 1.8 = 3.3 — more than twice the mean, confirming that substantial overdispersion is present. Had Poisson been fitted instead, the confidence interval around the IRR of 0.80 would have been artificially narrow (e.g., [0.71, 0.90] instead of [0.67, 0.95]), falsely suggesting higher precision.

Practical interpretation for a decision-maker: On average, patients receiving this treatment experienced approximately 20% fewer hospitalizations per year compared with the control group (rate ratio 0.80; 95% CI 0.67 to 0.95). Because this is an observational study, this is an association — unmeasured confounding could partially explain the difference — but the adjustment for measured covariates provides a more credible estimate than an unadjusted comparison.

The 95% confidence interval is consistent with reductions ranging from 5% to 33%, which is the honest uncertainty range after accounting for the variability in counts across patients. The dispersion parameter confirms that patients vary substantially in their hospitalization rates (some have zero, others have several per year), which is why negative binomial was used and why the confidence interval is wider than a Poisson analysis would have suggested.

Pros, cons, and trade-offs

Pros of the negative binomial distribution:

• Correctly models overdispersed count data: the NB variance function Var(Y) = μ + αμ² absorbs extra-Poisson heterogeneity, producing unbiased mean estimates and honest confidence intervals in the regression context.
• Theoretically motivated as the Poisson-gamma mixture: provides a mechanistic story for why counts are overdispersed in heterogeneous populations, linking the dispersion parameter to the magnitude of unobserved between-patient frailty.
• Interpretable IRR when used with a log link and offset; the rate-ratio scale maps directly to HEOR communication ("X% fewer events per patient-year").
• Graceful degradation as α → 0: the NB collapses to Poisson with no discontinuity in inference; no model change is needed if the data happen to be equidispersed.
• Available in all major statistical packages with mature implementations (MASS::glm.nb in R, NegativeBinomial in statsmodels, DIST=NEGBIN in SAS PROC GENMOD).

Cons of the negative binomial distribution:

• One additional parameter (α) to estimate, requiring a larger sample to pin down reliably; estimation can fail to converge when overdispersion is very mild (α near 0) or extreme.
• The boundary likelihood-ratio test for α = 0 requires halving the standard chi-square p-value; analysts who forget this step overstate the evidence against Poisson.
• Does not model event timing or ordering; the entire count history collapses to a single integer per patient, discarding the recurrence structure.
• Does not accommodate structural zeros without extension to ZINB or hurdle models, which add complexity and require two sets of coefficients to interpret.
• The NB-2 quadratic variance function can overfit the distributional tails in small samples where α cannot be estimated precisely.

Trade-offs between distributional choices:

• Versus Poisson: Poisson is simpler and requires no extra parameter, but its SEs are anticonservative under overdispersion — the classic error in naive HCRU analysis. NB is always at least as appropriate as Poisson and should be the default for healthcare counts. Use Poisson only when a formal dispersion test (LR boundary test; Pearson χ²/df ≈ 1) supports equidispersion.
• Versus quasi-Poisson: quasi-Poisson inflates SEs by an estimated scale factor but does not change the likelihood or provide a proper probability model. It cannot estimate α, does not support likelihood-based AIC/BIC model selection, and cannot generate predicted count distributions for simulation or power analysis. NB is generally preferred when a full probabilistic model is needed; quasi-Poisson is acceptable for robust SE correction in simple descriptive analyses.
• Versus recurrent-event survival models (Andersen-Gill, PWP): NB collapses event timing into a single count, sacrificing the at-risk structure and the ability to model gap-time dependence or the depletion of the at-risk set after terminal events. Survival models are more appropriate when the timing and ordering of events carry scientific meaning, or when competing death is severe and differential between arms.
• Versus ZINB/hurdle: appropriate only when a structural zero-generating process exists beyond what the NB naturally models. Over-applying zero-inflation to ordinary overdispersion creates two models where one suffices, reduces statistical power, and complicates interpretation.

When NOT to use

• Equidispersed processes: genuinely equidispersed counts where Pearson χ²/df ≈ 1.0 and the LR boundary test does not reject Poisson. In homogeneous clinical-trial populations or for rare events with very low α, Poisson is simpler and adds no unnecessary parameter.
• Structural zeros requiring hurdle or zero-inflated models: when the zero-event subgroup represents patients who structurally cannot generate the event — not just patients who happened to have zero events in the observation window — plain NB will misfit the zero probability. Use ZINB or hurdle NB after comparing observed versus NB-predicted zero frequencies.
• Time-varying rates and recurrent-event survival analysis: when the scientific question is the time-to-next-event, the gap-time between events, or the at-risk process over calendar time, the NB distribution discards exactly the information that matters. Use Andersen-Gill, Prentice-Williams-Peterson, or Wei-Lin-Weissfeld recurrent-event models.
• Duration or time-to-event outcomes: if the outcome is time (days to discharge, months to death, survival time), the NB distribution is inapplicable; use survival analysis (Cox, parametric AFT, or Fine-Gray competing-risks models).
• Continuous outcomes: the NB is defined only for non-negative integers; costs, laboratory values, and length-of-stay in fractional hours are not counts and require gamma GLM, log-normal, or Tweedie models.
• Uncontrolled confounding: like any GLM, an unadjusted NB comparison across treatment arms in an observational dataset estimates a biased association, not a causal rate ratio. Pair NB regression with propensity weighting, matching, or g-methods for valid causal inference.

Decision diagram

flowchart TD
  Het[Patients differ in unmeasured event-rate propensity]
  Het --> Gamma["Each patient's rate lambda_i drawn from Gamma(r, r/mu)"]
  Gamma --> PoisMix["Conditional on lambda_i, counts ~ Poisson(lambda_i)"]
  PoisMix --> NB["Marginal distribution of counts ~ Negative Binomial(mu, alpha=1/r)"]
  NB --> VF["Variance function: Var(Y) = mu + alpha*mu^2"]
  VF --> D1{alpha near 0?}
  D1 -->|Yes| Pois[Use Poisson: equidispersed]
  D1 -->|No| D2{Structural zeros beyond NB?}
  D2 -->|No| NBreg[Negative Binomial regression with offset]
  D2 -->|Yes| ZINB[ZINB or hurdle NB]
  NBreg --> IRR["exp(beta_j) = incidence rate ratio"]

Worked example

Scenario

A pharmacoepidemiologist is analyzing annual COPD exacerbation counts for 10 patients in a claims database, each observed for exactly one full year on the same treatment. Before fitting a count regression model, the analyst wants to check whether the data are overdispersed — that is, whether the sample variance substantially exceeds the mean, which would make Poisson inappropriate and negative binomial necessary. The 10 patients' observed annual exacerbation counts are recorded in the table below.

Dataset

Annual COPD exacerbation counts for 10 patients, each observed for exactly 1 year. Three patients had zero exacerbations; one patient had six, pulling the variance well above the mean.

Steps

Result

Mean = 20/10 = 2 exacerbations per patient per year. Sample variance = 36/9 = 4. Overdispersion ratio = 4/2 = 2. Because the variance (4) is double what the Poisson assumption requires (variance would equal the mean, 2), a negative binomial model is needed. The dispersion parameter alpha will be estimated from the data and inserted into the variance function Var(Y) = mean + alpha times mean squared to give honest confidence intervals. Using Poisson here would produce anticonservative standard errors and falsely narrow confidence intervals.

Trade-offs

Runnable example

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy import stats

# ── Motivating dataset: 10 COPD patients, annual exacerbation counts (1 year each) ──
df = pd.DataFrame({
    "patient_id"   : [f"P{i:02d}" for i in range(1, 11)],
    "count"        : [0, 0, 0, 1, 1, 2, 3, 3, 4, 6],
    "person_years" : [1.0] * 10,      # all observed exactly 1 year
})

# ── 1. Overdispersion diagnostic: sample mean and variance ──
mean_count  = df["count"].mean()          # 2.0
var_count   = df["count"].var(ddof=1)     # 4.0  (sample variance, n-1 denominator)
disp_ratio  = var_count / mean_count      # 2.0  (should be 1.0 under Poisson)
print(f"Mean: {mean_count:.2f}  |  Variance: {var_count:.2f}  |  Var/Mean ratio: {disp_ratio:.2f}")
print("Overdispersed (ratio > 1):", disp_ratio > 1)

# ── 2. Poisson rate model (intercept-only for illustration); offset = log(person_years)
df["log_pt"] = np.log(df["person_years"])
pois = smf.glm("count ~ 1", data=df,
                family=sm.families.Poisson(),
                offset=df["log_pt"]).fit(disp=False)

# Pearson chi-square / residual df: >> 1 flags overdispersion from the Poisson fit
pois_disp = float(pois.pearson_chi2 / pois.df_resid)
print(f"\nPoisson Pearson chi2/df = {pois_disp:.2f}  (expect ~1 if Poisson holds)")

# ── 3. Negative binomial MLE to estimate dispersion alpha ──
#    NegativeBinomial.from_formula: fits NB-2 Var(Y) = mu + alpha*mu^2
nb_mle = sm.NegativeBinomial.from_formula(
    "count ~ 1", data=df, offset=df["log_pt"]
).fit(disp=0, method="nm", maxiter=500, disp_kwds={"disp": False})
alpha = float(nb_mle.params["alpha"])
print(f"\nEstimated dispersion alpha = {alpha:.4f}")
print(f"Variance function: Var(Y) = mu + {alpha:.4f} * mu^2")

# ── 4. Boundary LR test of H0: alpha = 0 (Poisson) vs NB; halve the chi-square p-value
lr_stat = 2.0 * (nb_mle.llf - pois.llf)
lr_p    = 0.5 * stats.chi2.sf(lr_stat, df=1)   # halved because alpha=0 is on boundary
print(f"\nBoundary LR test: stat = {lr_stat:.2f}, p (halved) = {lr_p:.4f}")
print("Decision: use NB" if lr_p < 0.05 else "Decision: Poisson may suffice")

# ── 5. NB-2 GLM using the estimated alpha; exponentiate intercept = estimated rate ──
nb_glm = smf.glm("count ~ 1", data=df,
                  family=sm.families.NegativeBinomial(alpha=alpha),
                  offset=df["log_pt"]).fit()
rate = float(np.exp(nb_glm.params["Intercept"]))
ci   = np.exp(nb_glm.conf_int().loc["Intercept"])
print(f"\nNB estimated rate = {rate:.2f} events/year  95% CI [{ci[0]:.2f}, {ci[1]:.2f}]")

# ── 6. Extended template: NB rate regression with treatment arm and covariates ──
# Replace df with your analytic dataset (one row per patient, already cohort-built):
#   analytic: person_id, arm (0/1), count, person_years, age, sex, comorb_index
#
# alpha_est = <alpha from nb_mle above, or fit on analytic>
# nb_reg = smf.glm(
#     "count ~ C(arm) + age + sex + comorb_index",
#     data=analytic,
#     family=sm.families.NegativeBinomial(alpha=alpha_est),
#     offset=np.log(analytic["person_years"])
# ).fit(cov_type="cluster", cov_kwds={"groups": analytic["person_id"]})
#
# irr = np.exp(nb_reg.params)        # incidence rate ratios
# ci  = np.exp(nb_reg.conf_int())    # 95% CIs on rate-ratio scale
# print(pd.concat([irr.rename("IRR"), ci], axis=1))

library(MASS)
library(sandwich)
library(lmtest)

# ── Motivating dataset ──
df <- data.frame(
  patient_id   = paste0("P", sprintf("%02d", 1:10)),
  count        = c(0, 0, 0, 1, 1, 2, 3, 3, 4, 6),
  person_years = rep(1.0, 10)
)

# ── 1. Overdispersion diagnostic ──
mean_count  <- mean(df$count)          # 2.0
var_count   <- var(df$count)           # 4.0  (R uses n-1 by default)
disp_ratio  <- var_count / mean_count  # 2.0
cat(sprintf("Mean: %.2f | Variance: %.2f | Var/Mean ratio: %.2f\n",
            mean_count, var_count, disp_ratio))

# ── 2. Poisson rate model (intercept-only, with offset) ──
pois_fit <- glm(count ~ 1, family = poisson(link = "log"),
                data = df, offset = log(person_years))

# Pearson dispersion statistic: >> 1 signals overdispersion
pois_disp <- sum(residuals(pois_fit, type = "pearson")^2) / df.residual(pois_fit)
cat(sprintf("Poisson Pearson chi2/df = %.2f  (expect ~1 under Poisson)\n", pois_disp))

# ── 3. Negative binomial (NB-2) via MASS::glm.nb ──
#    glm.nb estimates theta (precision); alpha = 1/theta in Var(Y) = mu + alpha*mu^2
nb_fit <- glm.nb(count ~ 1 + offset(log(person_years)), data = df)
theta   <- nb_fit$theta
alpha   <- 1 / theta
cat(sprintf("NB theta (precision) = %.4f\n", theta))
cat(sprintf("NB alpha (dispersion) = %.4f  [Var(Y) = mu + %.4f * mu^2]\n", alpha, alpha))

# ── 4. Boundary LR test: H0 alpha = 0 (Poisson) vs NB; halve the p-value ──
lr_stat <- 2 * as.numeric(logLik(nb_fit) - logLik(pois_fit))
lr_p    <- 0.5 * pchisq(lr_stat, df = 1, lower.tail = FALSE)
cat(sprintf("Boundary LR test: stat = %.2f, p (halved) = %.4f\n", lr_stat, lr_p))
cat(if (lr_p < 0.05) "Decision: use NB\n" else "Decision: Poisson may suffice\n")

# ── 5. Estimated rate with 95% CI (intercept-only model = overall event rate) ──
rate <- exp(coef(nb_fit)["(Intercept)"])
ci   <- exp(confint(nb_fit))
cat(sprintf("NB rate = %.2f/year  95%% CI [%.2f, %.2f]\n", rate, ci[1], ci[2]))

# ── 6. Extended template: NB regression with treatment arm, covariates, offset ──
# Replace df with your analytic cohort dataset:
#   analytic: person_id (cluster), arm (factor), count, person_years, age, sex, comorb_index
#
# nb_reg <- glm.nb(
#     count ~ arm + age + sex + comorb_index + offset(log(person_years)),
#     data = analytic
# )
# # Cluster-robust SEs to account for within-person correlation:
# vc <- vcovCL(nb_reg, cluster = ~person_id)
# irr_table <- exp(coeftest(nb_reg, vcov. = vc)[, 1:2])   # IRR and SE
# cat("Incidence Rate Ratios:\n"); print(irr_table)

/* ── Create motivating dataset ── */
data work.exacerbations;
  input patient_id $ count person_years;
  log_pt = log(person_years);   /* offset = log(person-time) */
  datalines;
P01 0 1.0
P02 0 1.0
P03 0 1.0
P04 1 1.0
P05 1 1.0
P06 2 1.0
P07 3 1.0
P08 3 1.0
P09 4 1.0
P10 6 1.0
;
run;

/* ── 1. Overdispersion diagnostic: mean and variance ── */
proc means data=work.exacerbations mean var;
  var count;
  /* Mean = 2, Variance = 4 -> ratio = 2 -> overdispersed */
run;

/* ── 2. Poisson rate model; SCALE=PEARSON prints Pearson phi (chi2/df) ── */
proc genmod data=work.exacerbations;
  model count = / dist=poisson link=log offset=log_pt scale=pearson;
  /* Pearson dispersion (phi) >> 1 flags overdispersion requiring NB */
  /* The printed "Scale" parameter in PROC GENMOD output is sqrt(phi) */
run;

/* ── 3. Negative binomial (NB-2): Var(Y) = mu + alpha*mu^2 ── */
proc genmod data=work.exacerbations;
  model count = / dist=negbin link=log offset=log_pt;
  /* PROC GENMOD NB output: the DISPERSION parameter is k (shape = 1/alpha).
     To obtain alpha: alpha = 1 / DISPERSION from the printed "Scale" row.
     Confirm convergence in the iteration history. */
run;

/* ── 4. Extended template: NB regression with treatment, covariates, and offset ──
   Replace work.analytic with your cohort dataset (one row per patient).
   exp(ESTIMATE) for arm is the adjusted incidence rate ratio (IRR). */
proc genmod data=work.analytic;
  class arm (ref='0') sex;
  model count = arm age sex comorb_index / dist=negbin link=log offset=log_pt;
  estimate 'IRR: arm 1 vs 0' arm 1 -1 / exp;   /* exponentiated = IRR */
run;

/* ── 5. Cluster-robust (empirical) variance for within-person correlation ── */
proc genmod data=work.analytic;
  class arm (ref='0') sex person_id;
  model count = arm age sex comorb_index / dist=negbin link=log offset=log_pt;
  repeated subject=person_id / type=ind;   /* empirical (robust) sandwich variance */
  estimate 'IRR: arm 1 vs 0' arm 1 -1 / exp;
run;

Citations

1. [1]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.
2. [2]Lawless JF. Negative binomial and mixed Poisson regression. Canadian Journal of Statistics. 1987;15(3):209-225.

1. [3]Cameron AC, Trivedi PK. Regression-based tests for overdispersion in the Poisson model. Journal of Econometrics. 1990;46(3):347-364.

1. [4]Zeileis A, Kleiber C, Jackman S. Regression models for count data in R. Journal of Statistical Software. 2008;27(8):1-25.