Zero-Inflated and Hurdle Count Models

In plain language

Zero-inflated and hurdle count models are used when a dataset has far more zero-event patients than a standard negative binomial model can explain — for example, when 70% of patients had zero hospitalizations but the model predicts only 50%. These models handle zeros by positing two separate groups: patients who structurally cannot have the event (like a non-smoker in a smoking-cessation study) and patients who could but did not during the study window. The price of this richer model is two sets of coefficients that must be interpreted separately and combined to get a population-average estimate, which is why a simpler negative binomial model is usually the right starting point and should only be replaced when there is a genuine subject-matter reason to believe in a separate zero-generating process.

Excess zeros versus overdispersion: not the same problem

Two problems can create elevated zero counts in healthcare data, and they require fundamentally different solutions. Overdispersion — the condition where Var(Y) > E(Y) — arises when patients have heterogeneous underlying event rates. The negative binomial distribution handles overdispersion correctly: patients with very low rates will produce zero-count observations frequently by chance, and the NB dispersion parameter absorbs that extra mass at zero through its variance function. This entry does not re-teach NB regression; see negative-binomial-distribution and poisson-negative-binomial-count-models for that foundation.

Zero-inflation addresses a categorically different situation: a subgroup of patients that structurally cannot generate the event of interest, regardless of follow-up length or model specification. These patients' zeros are not probabilistic outcomes of a count process — they are certainties from a separate data-generating mechanism entirely. A plain NB model treats all zeros as draws from a single continuous count distribution. When a structural-zero subgroup genuinely exists, the NB consistently underpredicts the observed zero frequency even after alpha is estimated, and its coefficients average over two scientifically distinct groups in a way that obscures both.

The key diagnostic: after fitting an NB model, compare its predicted P(Y=0) to the observed zero proportion. If the NB predicts 50% zeros but you observe 70%, and there is a credible subject-matter reason to believe a true never-user subgroup exists, zero-inflated or hurdle models deserve consideration. If the NB tracks the zero rate well, the additional complexity is unjustified. Most of the time, plain NB suffices — this is the most important sentence in this entry.

Structural zeros versus sampling zeros

A structural zero is a zero that cannot be anything else. A patient who has undergone a hysterectomy contributes a structural zero to any study of uterine procedures. A patient who is physiologically incapable of addiction, or who is enrolled in a claims database but has never sought care for the relevant condition, may be a structural zero for drug fills or specialist visits. The defining feature: no amount of time in the at-risk window can change the count from zero, because the generating mechanism is absent.

A sampling zero is probabilistic. A patient who fills opioid prescriptions occasionally might have zero fills in any given year simply because no refill fell within the observation window. A low-frequency hospitalizer will frequently have zero hospitalizations in a one-year window. These zeros are fully explained by the count distribution at low event rates; the negative binomial handles them without special treatment.

The clinical distinction determines the entire modeling strategy. In administrative claims, structural zeros are rarely directly observable: a claims record for a patient with zero fills looks identical whether the patient is a structural never-user or a sampling zero from a sporadic user. Subject-matter knowledge must therefore drive the decision. If every patient in your cohort definition can plausibly generate the event — for example, elderly patients who can be hospitalized, diabetic patients who can fill glucose-lowering drugs — the structural-zero hypothesis has no support and zero-inflation adds a scientifically meaningless latent class.

ZIP and ZINB: mixture model mechanics and two coefficient sets

A zero-inflated model is a finite mixture of two components. Let pi (π) denote the probability that a given patient belongs to the structural-zero class. With probability π the patient contributes a guaranteed zero. With probability (1 − π) the patient is in the active-count class and their count follows a Poisson distribution (for ZIP) or a negative binomial distribution (for ZINB).

The observed probability of a zero count is therefore:

P(Y = 0) = π + (1 − π) × P(Y_count = 0)

and the probability of a positive count k > 0 is:

P(Y = k) = (1 − π) × P(Y_count = k)

The inflation term π is itself modeled as a function of patient covariates via a logistic regression: log(π/(1−π)) = γ₀ + γ₁Z₁ + ... where Z variables may differ from the count-model predictors. This produces two distinct coefficient sets that must be interpreted separately:

First, the inflation logistic coefficients (γ): these predict belonging to the structural-zero class. A positive γ coefficient means higher covariate values increase the probability of being a structural never-user. Exponentiated, they are odds ratios for the latent structural-zero state.

Second, the count-model coefficients (β): these predict the expected count for the active-count subgroup only. Exponentiated, they are incidence rate ratios (IRRs) conditional on the patient not being a structural zero. A treatment-arm IRR from the count component does not apply to all patients in the study — it applies only to the estimated (1 − π) fraction who are active-count patients.

The marginal expected count for any patient, combining both components, is:

E(Y) = (1 − π) × mu_count

where mu_count is the count-component mean. This marginal quantity — not the count-model IRR alone — is the right input to a budget-impact analysis, because it averages over the probability of being in each class.

Hurdle count models: a distinct data-generating story

A hurdle model is architecturally similar to a zero-inflated model — two components, two coefficient sets — but embodies a different causal story. In a hurdle model, there is no latent structural-zero class. All zeros have a single origin: the patient did not cross the utilization hurdle. A logistic model determines whether the count is zero (no event) or positive (any event). If positive, the count follows a zero-truncated count distribution, meaning the Poisson or NB probability mass function is renormalized to exclude zero.

Under this story:

P(Y = 0) = 1 − p_hurdle P(Y = k | k > 0) = p_hurdle × P(Y_truncated = k)

The hurdle is more natural for many utilization questions. Consider a study of specialty visits: some patients never seek the specialist at all (did not cross the hurdle), and among those who do, the number of visits follows a right-skewed count distribution. Every zero is explained by the logistic component alone; there is no latent-class ambiguity. The logistic hurdle coefficient is directly interpretable as a predictor of any engagement with care, which is itself a policy-relevant outcome separable from the volume question.

The practical difference between ZIP/ZINB and hurdle: for a given observed dataset, both models may produce similar predicted distributions and pass goodness-of-fit checks similarly. The difference is in interpretation and in the scientific commitment being made. Under ZIP, some patients are asserted to be structural zeros who cannot respond to treatment through the count channel; under hurdle, those same patients simply did not engage, and an intervention could move them across the hurdle. The choice should be grounded in the scientific question, not in model-selection statistics.

Model choice: theory first, Vuong test controversies, rootogram checks

The Vuong test has been widely applied to compare ZIP and plain Poisson (or ZINB and plain NB) in practice, but it is a poor gating criterion. The standard Vuong statistic does not follow the assumed standard-normal null distribution for non-nested model comparison in finite samples; simulation evidence shows that it generates inflated false-positive rates, recommending zero-inflation even when none exists. A correction term (Vuong-Clarke) restores proper calibration but is rarely implemented in standard software outputs. Because of these properties, automated Vuong-based model selection produces unreliable conclusions and encourages treating a scientific modeling decision as a mechanical test.

The recommended workflow:

1. Before looking at the data, ask whether structural zeros are scientifically plausible in this cohort. If the cohort definition requires the patient to be capable of the event, there are no structural zeros by design and the ZIP story is inappropriate.

2. Fit a plain negative binomial model and compute the predicted zero probability from the fitted model: P_NB(Y=0) = (1/(1+alpha*mu))^(1/alpha). Compare to the observed zero fraction. A rootogram (hanging rootogram from the countreg R package) displays this comparison visually for all count values simultaneously, making underfitting at zero obvious.

3. Only if a meaningful gap remains and subject matter supports a two-process story, fit ZINB or hurdle. Examine whether the inflation (or hurdle) coefficient estimates are interpretable and stable. A very large inflation probability with huge standard errors suggests the model cannot identify the latent class reliably.

4. Report the marginal expected count from the chosen model alongside the component-specific coefficients to give decision-makers the population-average quantity they actually need.

When plain negative binomial suffices — and that is usually the right answer

Analysts routinely reach for zero-inflated models upon seeing a high zero rate, but the NB distribution already produces substantial zero probability when mean rates are low and dispersion is high. The NB mass at zero is P_NB(Y=0) = (1/(1+alpha*mu))^(1/alpha), which can be 0.6 or higher at low mu with moderate alpha. Before adding model complexity, confirm quantitatively that the NB underpredicts the zero fraction by a practically meaningful margin — not just a margin that is statistically detectable at large sample sizes. In pharmacoepidemiology studies with 50,000-patient claims cohorts, even a 2-percentage-point underprediction of zeros will produce a "significant" Vuong test result without any scientific justification for structural zeros. The default should remain plain NB unless there is a specific, articulated scientific reason to believe in a separate zero-generating process.

Pros, cons, and trade-offs

Zero-inflated models (ZIP/ZINB): - Pros: correctly models a structural-zero subgroup when one genuinely exists; separates the probability of being a never-utilizer from the count rate among active users; can substantially improve predicted zero fit when the latent-class story is scientifically supported. - Cons: two coefficient sets required for every inference question; structural-zero class is latent and unverifiable from data alone; Vuong test for selection is unreliable; marginal expected counts require combining both components; overfitting is a real risk when the zero excess is ordinary overdispersion that NB would absorb; interpretation for payer and HTA audiences is substantially more complex. - When to prefer: population-based cohorts where genuine never-users are included by design; when the NB systematically underpredicts the zero frequency by more than ten percentage points and subject matter supports a structural-zero story.

Hurdle models: - Pros: cleaner data-generating story with all zeros from one source; logistic hurdle coefficient directly interpretable as a predictor of any utilization; avoids the unverifiable latent-class interpretation of ZIP/ZINB; separable scientific questions (any engagement vs how much engagement) map naturally to the two components. - Cons: the two-part structure still requires reporting and interpreting two models; forces all zeros into the logistic part even if some are sampling zeros from a low-rate count process; zero-truncated count distribution is less familiar and less commonly implemented in standard software. - When to prefer: the scientific question naturally divides into any-event versus how-many-given- event; specialty utilization, medication initiation, or referral studies where "crossing the hurdle" of first engagement is itself a meaningful outcome.

Plain negative binomial (the usual right answer): - Pros: single interpretable coefficient set; marginal IRR applicable to the full patient population; handles overdispersion and much of the zero excess without any latent-class assumption; computationally simpler and results are immediately communicable to non-technical stakeholders; the correct model in the vast majority of HCRU count analyses. - Cons: will underfit genuine structural zero processes; predicted zero probability from the fitted model should always be compared to the observed rate as a diagnostic. - When to prefer: almost always — specifically when the NB-predicted zero probability matches the observed zero fraction and when no subject-matter case for structural zeros can be made.

When to use

Use zero-inflated or hurdle count models when all of the following hold:

• The cohort definition plausibly includes patients who structurally cannot generate the count
• The fitted NB model's predicted zero probability materially underpredicts the observed zero
• A subject-matter story for the two-process model can be articulated and is scientifically
• The study audience can interpret and use two-component results, or marginal expected counts

When NOT to use

• When the NB fits the zero distribution well: if the predicted NB zero probability closely
• When the cohort is restricted to event-experiencing patients: inclusion criteria requiring
• When every patient can plausibly have the event: for endpoints that any enrolled patient
• As a mechanical response to a high zero rate: high zero rates are common in claims data
• When a single marginal rate ratio is required for communication or downstream modeling:

Interpreting the output

Using the worked example: a ZIP model for opioid prescription fills in 100 chronic pain patients produces two sets of estimated parameters. The inflation logistic component estimates π = 0.40 (40% are structural never-users). The count-component mean is μ_count with the Poisson zero probability for the active subgroup at 0.50.

(1) Formal interpretation. The ZIP model estimates two simultaneous processes. The logistic inflation component estimates that 40 out of 100 patients (π = 0.40) belong to a latent structural-zero class — patients who would produce zero opioid fills under any length of follow-up. These patients are not draws from the Poisson count distribution; their zero contributions are certain rather than probabilistic. Among the remaining 60% of patients (the active-count class, probability 1 − π = 0.60), the Poisson count component governs the fill distribution. The count-model IRR for a treatment indicator applies only within this active subgroup, not to all 100 patients. The marginal predicted zero rate for the full population is π + (1 − π) × P_count(Y=0) = 0.4 + 0.6 × 0.5 = 0.4 + 0.3 = 0.7, matching the observed 70/100 = 0.7. Any treatment-effect comparison that uses only the count-model IRR implicitly conditions on the patient not being a structural zero — an assumption that must be stated explicitly and is typically not what a budget-impact model requires.

(2) Practical interpretation. For a payer or clinical decision-maker: the model estimates that roughly 4 in 10 patients in this chronic pain population do not fill opioid prescriptions at all — a structural behavioral or clinical state the count model alone cannot explain. Among the 6 in 10 who do fill at least occasionally, the count model describes how often they fill and how a treatment or policy affects that rate. For a budget-impact calculation, the population-average expected fills per patient is E(Y) = (1 − 0.40) × μ_count, not μ_count alone — using μ_count directly would overestimate the cost burden for the full population by nearly 67%. This is the interpretive price of the richer model: two sets of coefficients, a marginal calculation to reach a population-level quantity, and a latent-class story that must be defended to a non-technical audience. When that price is not worth paying because the plain NB fits well, do not pay it.

Worked example

Scenario

A pharmacy analyst studies opioid prescription fills among 100 patients diagnosed with chronic pain over one year. Seventy patients had zero fills. After fitting a negative binomial model, the predicted proportion of zero-fill patients is approximately 50 percent — substantially below the observed 70 percent. With subject-matter justification (some patients decline opioids entirely on behavioral or cultural grounds), the analyst fits a zero-inflated Poisson model. The inflation logistic component estimates that 40 percent of patients are structural never-users (π = 0.40). The Poisson count component for the active-user group yields a zero probability of 0.50 at the estimated mean. The analyst verifies that these ZIP parameters reproduce the observed zero rate.

Dataset

Observed count distribution of opioid fills in one year for 100 chronic pain patients. Seventy patients had zero fills; the remaining 30 had one or more fills, creating a zero proportion of 0.70 that a plain NB model cannot match.

Steps

• Observed zero count = 70 patients out of 100 total. Observed zero rate: 70/100 = 0.7.

• A fitted negative binomial model predicts approximately 50 patients with zero fills (observed rate 0.70 vs NB-predicted rate approximately 0.50) — a gap suggesting a zero-generating process beyond ordinary overdispersion.

• Fit a ZIP model. The logistic inflation component estimates pi = 0.40: 40 percent of patients are classified as structural never-users who contribute a guaranteed zero.

• The Poisson count component for the remaining active-user group (probability 0.60) estimates a mean fill rate such that the count-part zero probability is 0.50 for an active-user with average characteristics.

• Compute the ZIP-predicted proportion of zeros. Structural-zero share contributes 0.40. Active-user sampling-zero share: 0.6*0.5 = 0.3. Total ZIP-predicted zero rate: 0.4 + 0.3 = 0.7.

• The ZIP prediction matches the observed rate: 0.4 + 0.3 = 0.7, which equals observed 70/100 = 0.7. The model successfully reproduces the excess zero rate by attributing 40 of the 70 zero-fill patients to the structural-zero class and 30 to sampling zeros from the active-user Poisson distribution.

Result

Observed zero rate: 70/100 = 0.7. ZIP-predicted zero rate: 0.4 + 0.3 = 0.7, matching observed. Of the 70 zero-fill patients, approximately 40 are attributed to the structural never-user class (pi = 0.40) and 30 to sampling zeros from the active-user group (0.6*0.5 = 0.3). The marginal expected fill count for the full population is (1 - 0.40) times the active-user mean — not the active-user mean alone — because 40 percent of patients are expected to contribute zero fills regardless of follow-up.

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

rng = np.random.default_rng(42)
n = 200

# ── Synthetic dataset: 200 patients, structural-zero subgroup (~40%) ──
# Structural zeros: 80 patients who never fill opioids
# Active users: 120 patients; fill count ~ NegBin(mu=2.5, alpha=0.8)
structural_zero = rng.binomial(1, 0.4, n).astype(bool)
mu_active = 2.5
alpha_nb   = 0.8
p_nb = 1.0 / (1.0 + alpha_nb * mu_active)         # NB prob param for scipy
r_nb = 1.0 / alpha_nb

active_counts = rng.negative_binomial(r_nb, p_nb, n)
counts = np.where(structural_zero, 0, active_counts).astype(int)

age = rng.normal(55, 12, n)
df = pd.DataFrame({
    "count":           counts,
    "age":             age,
    "pain_severity":   rng.uniform(0, 10, n).round(1),
    "person_years":    np.ones(n),
})
df["log_pt"] = np.log(df["person_years"])

obs_zero_rate = (df["count"] == 0).mean()
print(f"Observed zero rate: {obs_zero_rate:.3f}")

# ── Step 1: Fit plain NB and compare predicted vs observed zeros ──
nb_mle = sm.NegativeBinomial.from_formula(
    "count ~ age + pain_severity", data=df, offset=df["log_pt"]
).fit(disp=0)
alpha_est = float(nb_mle.params["alpha"])
mu_hat    = nb_mle.predict()

# P_NB(Y=0) = (1/(1 + alpha*mu))^(1/alpha) for each patient
p_zero_nb = (1.0 / (1.0 + alpha_est * mu_hat)) ** (1.0 / alpha_est)
nb_pred_zero_rate = p_zero_nb.mean()
print(f"NB predicted zero rate: {nb_pred_zero_rate:.3f}")
print(f"Gap (observed - NB): {obs_zero_rate - nb_pred_zero_rate:.3f}")
print("NB underpredicts zeros by more than 10 pp -> investigate zero-inflation")

# ── Step 2: Fit ZIP (as a diagnostic, count component is Poisson) ──
zip_res = sm.ZeroInflatedPoisson.from_formula(
    "count ~ age + pain_severity",      # count-model predictors
    exog_infl=df[["age"]].assign(const=1.0),   # inflation predictors
    data=df,
    offset=df["log_pt"]
).fit(method="bfgs", maxiter=400, disp=False)
print("\nZIP model summary (count component):")
print(zip_res.summary2().tables[1].round(3))

# ── Step 3: Fit ZINB (NB count component + inflation logit) ──
zinb_res = sm.ZeroInflatedNegativeBinomialP.from_formula(
    "count ~ age + pain_severity",
    exog_infl=df[["age"]].assign(const=1.0),
    data=df,
    offset=df["log_pt"]
).fit(method="bfgs", maxiter=400, disp=False)

# ── Step 4: Extract the structural-zero probability (pi) ──
# Inflation logit -> sigmoid -> pi for each patient
pi_hat = zinb_res.predict(which="prob-zero")  # or which="mean-main"
pi_mean = float(pi_hat.mean()) if hasattr(pi_hat, "mean") else float(pi_hat)
print(f"\nEstimated mean structural-zero probability (pi): {pi_mean:.3f}")

# ── Step 5: Compute marginal expected count E(Y) = (1 - pi) * mu_count ──
mu_count_hat = zinb_res.predict(which="mean-main")
marginal_mean = (1.0 - pi_mean) * float(mu_count_hat.mean())
raw_count_mean = float(mu_count_hat.mean())
print(f"Active-user mean (count component): {raw_count_mean:.3f}")
print(f"Marginal expected count E(Y):        {marginal_mean:.3f}")
print(f"Using count-component mean for full population would overstate by "
      f"{100 * (raw_count_mean / marginal_mean - 1):.0f}%")

# ── Step 6: Verify ZIP arithmetic from worked example ──
pi_ex     = 0.4       # structural-zero share
p_count_0 = 0.5       # count-part zero probability for active group
zip_pred_zero = pi_ex + (1 - pi_ex) * p_count_0
obs_ex        = 70 / 100
print(f"\nWorked-example verification:")
print(f"ZIP predicted zeros: {pi_ex} + {1-pi_ex}*{p_count_0} = {zip_pred_zero:.1f}")
print(f"Observed:            {obs_ex:.1f}")
print(f"Match: {abs(zip_pred_zero - obs_ex) < 1e-9}")

library(MASS)
library(pscl)
# install.packages("countreg", repos="http://R-Forge.R-project.org")
# library(countreg)   # for rootogram(); install from R-Forge if needed

set.seed(42)
n <- 200

# ── Synthetic dataset: structural-zero subgroup (~40% never fill) ──
structural_zero <- rbinom(n, 1, 0.4) == 1
active_counts   <- rnbinom(n, mu = 2.5, size = 1 / 0.8)  # size = 1/alpha
counts          <- ifelse(structural_zero, 0L, active_counts)

df <- data.frame(
  count        = counts,
  age          = round(rnorm(n, 55, 12), 1),
  pain_severity= round(runif(n, 0, 10), 1),
  person_years = 1.0
)

obs_zero_rate <- mean(df$count == 0)
cat(sprintf("Observed zero rate: %.3f\n", obs_zero_rate))

# ── Step 1: Fit plain NB; check predicted vs observed zeros ──
nb_fit <- glm.nb(count ~ age + pain_severity + offset(log(person_years)), data = df)
mu_hat <- predict(nb_fit, type = "response")
theta  <- nb_fit$theta            # MASS parameterization: size = theta = 1/alpha
# P_NB(Y=0) = (theta / (theta + mu))^theta for each patient
p_zero_nb <- (theta / (theta + mu_hat)) ^ theta
nb_pred_zero_rate <- mean(p_zero_nb)
cat(sprintf("NB predicted zero rate: %.3f\n", nb_pred_zero_rate))
cat(sprintf("Gap: %.3f\n", obs_zero_rate - nb_pred_zero_rate))

# Rootogram for visual diagnostic (requires countreg):
# rootogram(nb_fit, style = "hanging", main = "NB fit: hanging rootogram")

# ── Step 2: Fit ZINB with the same predictors for count and inflation parts ──
# Formula syntax: count-model predictors | inflation-model predictors
zinb_fit <- zeroinfl(
  count ~ age + pain_severity + offset(log(person_years)) | age,
  data = df,
  dist = "negbin"
)
cat("\nZINB coefficient summary:\n")
print(summary(zinb_fit)$coefficients)

# ── Step 3: Extract structural-zero probability pi per patient ──
pi_hat <- predict(zinb_fit, type = "zero")   # P(structural zero) per patient
cat(sprintf("\nMean estimated pi (structural-zero probability): %.3f\n", mean(pi_hat)))

# ── Step 4: Marginal expected count E(Y) = (1 - pi) * mu_count ──
mu_count_hat <- predict(zinb_fit, type = "count")   # conditional on active-user class
marginal_mean <- mean((1 - pi_hat) * mu_count_hat)
raw_count_mean <- mean(mu_count_hat)
cat(sprintf("Active-user mean (count component):  %.3f\n", raw_count_mean))
cat(sprintf("Marginal expected count E(Y):         %.3f\n", marginal_mean))
cat(sprintf("Overstatement if count mean used for full population: %.0f%%\n",
            100 * (raw_count_mean / marginal_mean - 1)))

# ── Step 5: Hurdle NB — alternative when all zeros are non-engagers ──
hurdle_fit <- hurdle(
  count ~ age + pain_severity + offset(log(person_years)) | age,
  data = df,
  dist = "negbin"
)
cat("\nHurdle NB (logistic hurdle + zero-truncated NB count):\n")
cat("Hurdle (logistic) coefficients:\n")
print(coef(hurdle_fit, model = "zero"))
cat("Count (zero-truncated NB) coefficients:\n")
print(coef(hurdle_fit, model = "count"))

# ── Step 6: Vuong test — informational only; do NOT use as sole decision criterion ──
# vuong(nb_fit, zinb_fit)   # from pscl; note the inflated false-positive rate issue
cat("\nNOTE: Vuong test is shown here for illustration only. Its null distribution\n")
cat("has inflated false-positive rates for non-nested ZI vs plain-NB comparison.\n")
cat("Use the rootogram and subject-matter knowledge as primary decision criteria.\n")

# ── Step 7: Verify worked-example arithmetic ──
pi_ex     <- 0.4
p_count_0 <- 0.5
zip_pred_zero <- pi_ex + (1 - pi_ex) * p_count_0
cat(sprintf("\nWorked-example: ZIP predicted zeros = %.1f + %.1f*%.1f = %.1f\n",
            pi_ex, 1 - pi_ex, p_count_0, zip_pred_zero))
cat(sprintf("Observed:                                           = %.1f\n", 70/100))