Poisson and Negative Binomial Count Models for HCRU and Utilization

In plain language

Poisson and negative binomial models count how many times an event happens to a patient over a period of time — for example, how many times they visited the ER or were admitted to the hospital — and compare those rates between two groups. Instead of asking "did the event happen?", they ask "how often did it happen?" and report the answer as a rate ratio (called an incidence rate ratio, or IRR) that accounts for the fact that some patients were followed longer than others. Poisson is the simpler version, but real healthcare data almost always violates its key assumption: when the variance in a count dataset is larger than its mean — a condition called overdispersion — Poisson produces confidence intervals that are misleadingly narrow, making weak findings look statistically significant. Negative binomial adds one extra parameter to absorb that extra spread, giving you honest uncertainty estimates, which is why it is the safer default for healthcare utilization data.

Poisson and negative binomial (NB) regression are the workhorse models for count outcomes in real-world HEOR: the number of hospitalizations, emergency department visits, outpatient encounters, infusions, distinct drug classes, or inpatient days a patient accrues over an observation window. Both use a log link, `log(μ) = β0 + β'X + log(person_time)`, so the exponentiated coefficient `exp(βj)` is a rate ratio — an incidence rate ratio (IRR) once person-time enters as an `offset`. The IRR is the quantity HCRU studies report ("the new therapy was associated with a 28% lower rate of all-cause hospitalization, IRR 0.72").

Core conceptual distinction

The two models share the same mean structure and the same IRR estimand; they differ only in the assumed mean–variance relationship, and that difference is doing all the work. (1) Poisson fixes Var(Y) = μ (equidispersion). This is almost never true for HCRU: a minority of frail high-utilizers create a long right tail, so the variance exceeds the mean (overdispersion). Under overdispersion the Poisson point estimate of the IRR stays roughly unbiased, but its standard errors are too small, confidence intervals too narrow, and p-values anticonservative — you declare effects "significant" that are not. (2) Negative binomial introduces a dispersion parameter α so that Var(Y) = μ + αμ² (the canonical NB-2 parameterization), absorbing the extra heterogeneity and widening the CI to its honest width; as α → 0 NB collapses back to Poisson. A formal likelihood-ratio test of H0: α = 0 (a boundary test, halve the p-value) decides between them. The offset is the second non-negotiable: it is `log(person_time)`, not a covariate — its coefficient is fixed at 1 — and omitting it silently turns an IRR back into a raw count ratio that confounds the effect with differential follow-up.

Pros, cons, and trade-offs

- vs cox-ph-regression (and recurrent-event survival models, Andersen–Gill / PWP / WLW): Count models target the total volume or rate of events and use the full count (0, 1, 2, …) in a single number per patient; they are simpler and more powerful when the scientific question is "how much utilization" rather than "how soon." Cost: they discard the timing and ordering of events and the at-risk structure, so they cannot represent time-varying exposure, the depletion of the at-risk set after death, or the gap-time dependence of recurrent events. Prefer count models when the endpoint is a rate or total over a fixed window; prefer recurrent-event survival when timing matters or competing death heavily truncates person-time (see recurrent-events-analysis-rwe). - vs logistic-regression-for-binary-outcomes: Counting events (number of admissions) uses information that "any admission, yes/no" throws away — more efficient and usually more clinically meaningful for HCRU. Cost: when only the first event matters clinically, or when the outcome is genuinely dichotomous, logistic is simpler and avoids modeling a count distribution. Prefer count models when the magnitude of utilization is the point. - vs healthcare-costs-pppm-pppy-pmpm (cost GLMs): Counts are discrete events; a one-night and a thirty-night admission are both "one admission." Dollars are continuous, heavy-tailed, and zero-inflated, and belong in a gamma/log-link or two-part model. Prefer count models for utilization volume and pair them with cost models for the full economic picture — never substitute one for the other. - Poisson vs NB within the family, and vs zero-inflated/hurdle: Plain Poisson is defensible only for rare, equidispersed counts; NB is the safe default for HCRU. When zeros are far more frequent than even NB predicts (a structural subgroup that cannot generate the event — e.g., patients who never re-engage with the system after a procedure), a zero-inflated NB (ZINB) or hurdle model separates a "who has any utilization" logistic process from a "how much among utilizers" count process. The cost of ZINB/hurdle is two coefficient sets to interpret and report.

When to use

The endpoint is a count of discrete events over a defined window (admissions, ED visits, infusions, distinct NDCs, inpatient days); follow-up varies across patients (then carry `log(person_time)` as an offset and report an IRR); multiple events per patient are common and their total is of interest; you want a single interpretable rate ratio that maps cleanly to budget-impact and value narratives. NB is the default; reach for Poisson only after a dispersion test supports equidispersion, and quasi-Poisson when you want robust dispersion-scaled SEs without committing to the NB likelihood.

Interpreting the output

Consider the worked example: Drug A generates 50 ED visits in 100 person-years (rate 0.50/year); Drug B generates 30 visits in 100 person-years (rate 0.30/year). The crude IRR = 0.50 / 0.30 = 1.67. After fitting a negative binomial model with covariate adjustment, suppose the reported adjusted IRR is 1.54 (95% CI 1.18–2.01).

Formal interpretation: Among patients starting Drug A, the adjusted incidence rate of ED visits was 1.54 times the rate in Drug B patients per person-year of follow-up, holding baseline covariates fixed and accounting for the person-time offset. The offset log(person_years) is what converts raw counts into rates — omitting it would produce a count ratio confounded with differential follow-up, not an IRR. The negative binomial model is reported here rather than Poisson because individual patient counts showed variance well above the mean (overdispersion); the Poisson CI for the same data would have been fictitiously narrow, producing anticonservative inference.

Practical interpretation: Drug A patients accrue approximately 54% more ED visits per year than Drug B patients after adjustment. Because the NB model absorbs overdispersion from a minority of high-utilizers, the 95% CI of 1.18–2.01 is the honestly calibrated uncertainty range — not the artificially narrow interval a Poisson fit would have reported. In budget-impact terms, multiply the IRR by the per-visit cost to estimate the incremental HCRU burden per person-year. Always verify the dispersion decision (LR test of alpha = 0) before reporting and note whether any zero-inflation adjustment was applied.

When NOT to use — and when it is actively misleading or dangerous

- Overdispersed data fit with plain Poisson. This is the classic error: the IRR looks fine but the CI is fictitiously narrow, so a null effect is reported as protective with "p<0.05." Always check Pearson χ²/df (≫1 signals overdispersion) and test α; if in doubt, report NB or quasi-Poisson SEs. - No offset when follow-up differs by arm. If the exposed arm disenrolls or dies earlier (shorter person-time), modeling raw counts without `log(person_time)` makes the exposed look lower-utilizing purely because they were observed less — immortal-time and informative-censoring artifacts masquerading as a treatment benefit (see immortal-time-bias-handling). - Within-person correlation ignored. Multiple events from the same patient (or panel counts across periods) violate independence; naive SEs are too small. Use cluster-robust/GEE variance or a random-effects (NB-mixed) model. - Competing death treated as "zero utilization." In elderly or oncology claims, a patient who dies early accrues few events not because therapy reduced utilization but because they left the at-risk set. Differential mortality by arm biases the IRR; censor person-time at death and consider a recurrent-event-with-terminal-event or joint model. - The real question is binary or about timing. If only the first event is clinically meaningful, or the policy question is time-to-readmission, a count model answers the wrong question.

Data-source operational depth

- Claims (FFS): The substrate of most HCRU count studies. Build person-time from enrollment spans, censoring at disenrollment, death, and data end; count events from all relevant claim files (inpatient facility + professional, outpatient, carrier) using validated code lists, because partial capture undercounts. Failure mode unique to claims: Medicare Advantage (MA) encounter data are notoriously incomplete and inconsistently submitted, so MA person-time can show artificially low hospitalization counts versus fee-for-service Parts A/B — restrict to FFS A/B (or commercial with a real medical benefit) and exclude MA-only person-time, or the count comparison is biased by data completeness rather than care. Mail-order/90-day fills and stockpiling distort `days_supply`-derived counts of "drug fills." - EHR: Counts come from orders, encounters, or NLP, and suffer visit-driven (informed-presence) bias: sicker patients generate more records simply by being seen more, inflating counts in both arms and especially in the sicker arm. Adjust for visit intensity, or link to claims for complete capture. A patient who leaves the health system is differentially lost, so define the observation window explicitly and treat loss as potentially informative. - Registry: Often carries clean structured counts (lines of therapy, transfusions, cycles) and is excellent for validating claims-based count algorithms, but is typically thin on complete utilization across all settings; link to claims for full event capture and to a death index for censoring. - Linked claims–EHR–vital records: The ideal substrate (EHR severity for confounding control + claims completeness for counts + reliable mortality for censoring), but linkage selection and order/fill/service-date discrepancies must be reconciled before person-time is built.

Worked claims example

Question: rate of all-cause inpatient hospitalization on an SGLT2 inhibitor vs a DPP-4 inhibitor among adults with type 2 diabetes in a commercial + Medicare FFS database. (1) New-user, active-comparator cohort: ≥365 days continuous A/B/D (or commercial medical+pharmacy) enrollment with no fill of either class in the lookback; time zero = first qualifying `fill_date`; arm = the NDC dispensed that day. (2) Person-time: from time zero to the earliest of disenrollment, death, or 2024-12-31 end of data — for patient A this is 365 days (1.00 person-year), for patient B who dies at day 200, 0.55 person-years. (3) Outcome count: number of distinct inpatient admissions in the window from the facility claim file, collapsing same-stay transfer claims so a transfer is not double-counted; patient A has 2 admissions, patient B has 1. (4) Offset: model `hosp_count ~ arm + age + sex + baseline_comorbidity` with `offset = log(person_years)`, so `exp(β_arm)` is the IRR (rate per person-year). (5) Distribution: fit Poisson, then NB; the LR test of α gives p<0.001 (overdispersion from a few high-utilizers), so report the NB IRR — say 0.78 (95% CI 0.69–0.88) — with cluster-robust SEs because some patients contribute multiple admissions. (6) Critically, exclude any MA-only enrollment months: a patient enrolled in MA for part of follow-up has incomplete inpatient capture there, which would deflate their count and the resulting person-time, biasing the IRR. (7) Sensitivity: vary the window, censor at treatment discontinuation (last `days_supply` end + grace period) for an as-treated analysis, and run a ZINB if the proportion of patients with zero admissions exceeds NB's prediction.

Worked example

Scenario

A database study compares rates of emergency department (ED) visits between adults newly starting Drug A versus Drug B for a chronic condition. Each patient is enrolled for a fixed one-year window. You want to compute the crude ED visit rate per person-year in each arm and the incidence rate ratio (IRR), and then decide whether Poisson or negative binomial is more appropriate.

Dataset

Summary counts for each treatment arm over one year of follow-up.

Steps

• Compute the rate for Drug A: 50 ED visits divided by 100 person-years = 0.50 ED visits per person-year.

• Compute the rate for Drug B: 30 ED visits divided by 100 person-years = 0.30 ED visits per person-year.

• Compute the incidence rate ratio (IRR): 0.50 divided by 0.30 = 1.67, meaning Drug A patients had ED visits at 1.67 times the rate of Drug B patients.

• Decide between Poisson and negative binomial: inspect the spread of individual patient counts in each arm. If a few patients account for the majority of visits — common in healthcare data — the variance will exceed the mean (overdispersion). In that case, fitting Poisson will produce a confidence interval around the IRR of 1.67 that is too narrow, making the difference look more certain than the data support. Fitting negative binomial instead lets the model absorb that extra spread and gives an honest, wider confidence interval.

Result

Drug A rate = 50 / 100 = 0.50 ED visits per person-year. Drug B rate = 30 / 100 = 0.30 ED visits per person-year. IRR = 0.50 / 0.30 = 1.67. When individual patient counts show variance greater than the mean (overdispersion), report the negative binomial IRR with its wider, correctly calibrated confidence interval rather than the Poisson result.

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

def fit_count_models(analytic: pd.DataFrame, formula: str = "hosp_count ~ arm + age + sex + baseline_comorb_index"):
    df = analytic[analytic["person_years"] > 0].copy()
    offset = np.log(df["person_years"].to_numpy())

    # Poisson rate model (IRR = exp(beta)); cluster-robust SE on person_id.
    pois = smf.glm(formula, data=df, family=sm.families.Poisson(), offset=offset).fit(
        cov_type="cluster", cov_kwds={"groups": df["person_id"]}
    )

    # Overdispersion check: Pearson chi-square / residual df should be ~1 under Poisson.
    disp = float(pois.pearson_chi2 / pois.df_resid)

    # Negative binomial (NB-2): estimate dispersion alpha, then refit GLM with that alpha.
    nb_mle = sm.NegativeBinomial.from_formula(formula, data=df, offset=offset).fit(disp=0)
    alpha = float(nb_mle.params["alpha"])
    nb = smf.glm(formula, data=df, family=sm.families.NegativeBinomial(alpha=alpha), offset=offset).fit(
        cov_type="cluster", cov_kwds={"groups": df["person_id"]}
    )

    # Boundary LR test of H0: alpha = 0 (Poisson) vs NB; halve the chi-square(1) tail probability.
    lr = 2.0 * (nb_mle.llf - pois.llf)
    lr_p = 0.5 * stats.chi2.sf(lr, df=1)

    def irr_table(fit):
        out = pd.DataFrame({"IRR": np.exp(fit.params)})
        ci = np.exp(fit.conf_int())
        out["ci_low"], out["ci_high"] = ci[0], ci[1]
        return out

    return {
        "pearson_dispersion": disp,
        "nb_alpha": alpha,
        "lr_stat": lr,
        "lr_p_boundary": lr_p,
        "model_choice": "negative_binomial" if lr_p < 0.05 else "poisson",
        "poisson_irr": irr_table(pois),
        "nb_irr": irr_table(nb),
    }

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

fit_count_models <- function(analytic,
                             fml = hosp_count ~ arm + age + sex + baseline_comorb_index) {
  df <- subset(analytic, person_years > 0)
  off <- log(df$person_years)

  # Poisson rate model with person-time offset.
  pois <- glm(fml, family = poisson(link = "log"), data = df, offset = off)

  # Overdispersion: Pearson chi-square / residual df (~1 under Poisson).
  disp <- sum(residuals(pois, type = "pearson")^2) / df.residual(pois)

  # Negative binomial (NB-2).
  nb <- glm.nb(update(fml, . ~ . + offset(log(person_years))), data = df)

  # Boundary LR test of theta -> Inf (Poisson) vs NB; halve the p-value.
  lr   <- 2 * (logLik(nb) - logLik(pois))
  lr_p <- 0.5 * pchisq(as.numeric(lr), df = 1, lower.tail = FALSE)

  # Cluster-robust (by person_id) IRRs with 95% CI.
  irr <- function(fit) {
    vc  <- vcovCL(fit, cluster = df$person_id)
    est <- coeftest(fit, vcov. = vc)
    data.frame(IRR = exp(est[, 1]),
               ci_low  = exp(est[, 1] - 1.96 * est[, 2]),
               ci_high = exp(est[, 1] + 1.96 * est[, 2]))
  }

  list(pearson_dispersion = disp,
       nb_theta = nb$theta,
       lr_p_boundary = as.numeric(lr_p),
       model_choice = ifelse(lr_p < 0.05, "negative_binomial", "poisson"),
       poisson_irr = irr(pois),
       nb_irr = irr(nb))
}