Two-Part and Hurdle Models for Semicontinuous Costs

In plain language

Many healthcare datasets have a large fraction of patients with exactly zero cost — they were enrolled and observable but never used the service being measured. A two-part model handles this by fitting two separate equations: the first uses logistic regression to estimate the probability that any cost occurs at all, and the second uses a gamma regression model (fitted only on the patients who did have costs) to estimate how much those users spent. Multiplying the two predictions gives the expected cost for any patient or group, on the same dollar scale that budget and payer models need. A key advantage is that a treatment can affect who seeks care and how much users spend by different amounts — effects that are hidden in a single model.

The semicontinuous cost distribution and why it breaks one-part models

Healthcare cost data from claims and linked sources routinely exhibit a structural pattern that statisticians call semicontinuous: a large mass of exact zeros from patients who had no qualifying utilization during the measurement window, combined with a positively skewed continuous distribution among those who did incur costs. The zeros are not missing data, rounding artifacts, or deductible truncation — they represent genuine non-users: patients who were enrolled and observable but who never filled a relevant prescription, visited an outpatient specialist, or were admitted to hospital during the study window. In a typical commercially insured cohort with an uncommon specialty drug, 30–60% of observed patients may have zero total cost in any given month.

This creates a fundamental problem for one-part models. The gamma distribution, the standard for right-skewed continuous costs (see gamma-distribution entry), is defined only on strictly positive values (Y > 0, ∞). Fitting a gamma GLM to a dataset that includes zeros ignores those observations or provokes convergence failures. Log-ordinary-least-squares cannot accommodate a zero outcome without a small offset constant, and the offset distorts both the distributional shape and the estimated mean. Tobit models (censored regression) treat zeros as left-censored observations, which is the wrong assumption: a zero cost is a true zero, not a cost below a detection limit. Only a two-part architecture directly models both mechanisms — the binary any-cost decision and the conditional cost level — as separate but compatible statistical processes.

The two-part model: Part 1 (any cost?) and Part 2 (how much?)

The two-part (hurdle) model factorizes the distribution of Y into two estimating equations:

Part 1. A logistic regression models P(Y > 0 | X) — the probability that a patient incurs any cost at all — as a function of covariates X. The binary outcome is 1{Y > 0}. This part answers the clinical question: "Does the treatment or exposure change who enters the healthcare system?" The logistic coefficients are log-odds ratios for any-cost occurrence. For the mechanics of logistic regression, see the logistic-regression-for-binary-outcomes entry.

Part 2. A gamma GLM with a log link (or a log-normal GLM) is fitted only among patients with Y > 0, modelling E[Y | Y > 0, X]. This part answers: "Among those who do use care, does the treatment change how much they spend?" Because Part 2 is restricted to the positive sample, the gamma GLM's requirement that Y > 0 is automatically satisfied. For the gamma GLM mechanics, variance function V(μ) = φμ², the Modified Park test, and coefficient interpretation, see the gamma-distribution entry — this entry does not duplicate those details.

Combining the two parts into the overall expected cost

The unconditional (population-average) expected cost integrates the two estimating equations:

E[Y | X] = P(Y > 0 | X) × E[Y | Y > 0, X]

In practice this is computed by recycled (marginal) standardisation: for each patient in the dataset, predicted probabilities from Part 1 are multiplied by predicted positive costs from Part 2, and these products are averaged over the sample to obtain the arm-specific expected cost. The treatment contrast is the difference (or ratio) of these marginalised expectations. This approach is estimand-transparent: the analyst recovers E[Y], the arithmetic mean on the dollar scale, which aggregates correctly to population totals for budget-impact and payer decision models.

Standard errors for the combined effect: delta method and bootstrap

Because the overall expected cost is a nonlinear function of coefficients from two separately fitted models, standard errors cannot be read directly from either model's output table. Two approaches are standard:

Delta method (analytical). By the multivariate delta method, the variance of E[Y | X] can be approximated from the gradient of the estimand with respect to the stacked coefficient vector and the asymptotic variance-covariance matrices of the two models. This is implemented in the margins and predictNL packages in R, and in recycled-prediction commands in Stata. The delta method is fast but relies on asymptotic normality and can understate variance in small or highly heterogeneous samples.

Bootstrap (nonparametric). The nonparametric bootstrap — re-fitting both parts on each bootstrap resample and computing the combined expected cost each time — captures variance from both models and from the dependence between them without distributional assumptions. For details and implementation of the bootstrap itself, see the bootstrap-resampling-methods entry. The bootstrap is the more robust default for complex two-part estimands and is recommended whenever computing resources permit (typically 500–2000 resamples).

The interpretive payoff: covariate effects that differ by part

The most powerful feature of the two-part model — and the primary reason to prefer it over a one-part Tweedie alternative in HEOR — is that it allows covariates to have entirely different effects in each part. This decomposition is clinically interpretable and practically important:

• A novel specialty drug may shift P(any cost) upward because initiating therapy requires an
• A care-management programme may not change who enters the system (Part 1 coefficient near zero)
• A comorbidity burden score may be a strong predictor of any utilisation (Part 1) but a weak

Reporting part-specific coefficients alongside the combined expected cost is therefore not optional detail — it is the substantive contribution that motivates choosing the two-part architecture in the first place.

Tweedie GLM as the one-model alternative

The Tweedie distribution is a compound Poisson-gamma family that is defined at zero (the Poisson component generates the zero mass) and on the positive real line (the gamma component generates the positive costs). A GLM with a Tweedie variance function — V(μ) = φμ^p where the power parameter p is estimated from data, typically 1 < p < 2 for cost data — handles zeros and positive costs in a single estimating equation. This is computationally simpler than fitting two separate models and produces a single coefficient per covariate that captures the combined effect on E[Y].

The trade-off is that the Tweedie constrains the covariate effect to be the same in both the probability-of-any-cost part and the conditional-cost part. When this constraint is reasonable (which must be checked by comparing Tweedie predictions to two-part predictions), the Tweedie is more parsimonious. When covariates genuinely differ in their Part 1 and Part 2 effects — the typical situation in drug-comparative HEOR — the Tweedie blends the two effects into one coefficient, hiding the mechanism and potentially misrepresenting the effect size. The two-part model is therefore preferred when the scientific question requires understanding how treatment affects each part separately, and the Tweedie is a reasonable default when only the combined expected cost is needed and the constraint can be empirically validated.

Hurdle vs two-part: terminology

In applied HEOR literature, "hurdle model" and "two-part model" are used interchangeably for the logistic + GLM architecture described here. A technical distinction sometimes drawn is that the hurdle model uses a zero-truncated distribution (e.g., zero-truncated gamma, zero-truncated Poisson) for the positive part, while the "two-part" label is reserved for any distribution on the positive part. In practice, when the positive-part model is a gamma or log-normal GLM — both defined on Y > 0 — the two terms describe identical models. Analysts should specify the positive-part distribution explicitly rather than relying on terminology alone.

Zero-inflated count cousins

When the outcome is a count (emergency department visits, hospitalisations, prescription fills) rather than continuous cost, the analogous model is a zero-inflated Poisson or zero-inflated negative-binomial regression. These share the two-mechanism logic (structural zero mass plus a count distribution for positive values) but are fitted with count-specific methods. For count outcomes, see the zero-inflated-count-models entry. The current entry is specifically for semicontinuous cost and other positive-continuous outcomes with a structural zero mass.

When a one-part gamma is sufficient

When the proportion of zero-cost patients is small — typically below 5–10% of the analytic cohort — the bias from fitting a one-part gamma GLM is modest. In this case, a constant small offset (e.g., Y + 0.01) or quasi-Poisson as a workaround is acceptable for a sensitivity analysis, though the two-part model remains the principled primary choice. The 5–10% threshold is a rule of thumb, not a hard cutoff; the decision should also weigh whether the zeros are informative (structural non-users vs. random non-claimants) and whether the part-specific decomposition is scientifically important.

Pros, cons, and trade-offs

Pros: - Correctly handles the semicontinuous distribution that is the empirical norm for healthcare cost data in claims, EHR, and registry sources with variable follow-up. - E[Y] = P(Y > 0) × E[Y | Y > 0] produces an unbiased estimate of the arithmetic mean on the dollar scale — the quantity that payers, HTA bodies, and budget models need. - Covariate effects can differ by part, enabling mechanistic interpretation (who uses care vs. how much users spend) that is invisible in one-part or Tweedie models. - Both parts (logistic and gamma GLM) are standard methods with well-understood diagnostics, implemented natively in all major statistical software. - Marginal standardisation of combined predictions accommodates IPTW weighting, subgroup effects, and covariate adjustment in a single coherent framework.

Cons: - Two separate models must be specified, diagnosed, and reported; twice the model-checking burden relative to a one-part Tweedie. - Standard errors for the combined estimand require bootstrap or delta method — they are not printed automatically by any model summary table. - The Part 2 gamma GLM is fitted only on the positive subsample, which reduces effective sample size and raises standard errors, especially when the zero fraction is large. - Separation or near-separation in Part 1 (all zeros in a covariate stratum) causes logistic instability; data sparsity in Part 2 (very few positives in a covariate stratum) causes gamma GLM instability. - When the Tweedie constraint (common covariate effects across both parts) holds approximately, the two-part model is slightly less efficient than the Tweedie.

When to use

Use the two-part model when: - The outcome is a continuous cost or resource-use measure with 10% or more zero values in the analytic cohort. Below 5%, a one-part gamma is often adequate as a sensitivity analysis; above 10%, the two-part model is the principled primary analysis. - The scientific question requires distinguishing effects on utilisation (Part 1) from effects on cost per user (Part 2) — a drug comparison, programme evaluation, or subgroup analysis where the mechanism matters as much as the overall expected cost. - The target estimand is the arithmetic mean E[Y] on the dollar scale, required for budget-impact models, payer dossiers, cost-effectiveness submissions, and value frameworks. - IPTW-weighted or stratified analyses where marginal standardisation of the combined prediction is the appropriate final step after weighting.

When NOT to use

• Very few zeros (<5%): when the zero fraction is small, fitting a one-part gamma GLM (possibly
• Count outcomes with excess zeros: zero-inflated Poisson or zero-inflated negative-binomial
• When the estimand is not the mean: if the policy question concerns the median cost, use
• When zeros are informative censoring, not structural zeros: if patients with zero cost are
• When causal language is required but confounding is uncontrolled: the two-part model is an

Interpreting the output

In the worked example below, the control arm has P(any cost) = 0.6 and E[cost | cost > 0] = 5000, so E[Y] = 0.6 5000 = 3000. The treated arm has P(any cost) = 0.5 and E[cost | cost > 0] = 4800, so E[Y] = 0.5 4800 = 2400. The difference is 3000 - 2400 = 600.

(1) Formal statistical interpretation. The two-part model produces two component estimates: a logistic odds ratio quantifying the change in the log-odds of any cost occurrence between arms (Part 1), and a gamma GLM mean cost ratio quantifying the change in expected cost among users between arms (Part 2). The combined estimand E[Y] = 3000 (control) versus E[Y] = 2400 (treated) is the marginalised (population-average) expected annual cost per patient, averaging over both non-users and users in each arm. The difference of $600 is the estimated adjusted mean cost reduction per patient per year in the treated arm. A 95% bootstrap confidence interval (from resampling both models simultaneously) represents values of the true mean cost difference that are compatible with the data at the conventional significance threshold in repeated sampling — it is not a Bayesian probability statement about the true difference.

(2) Practical interpretation for a decision-maker. "On average, treated patients cost $600 less per year than control patients — but this composite hides two separate effects. In the treated arm, fewer patients sought care at all (50% vs 60% in controls — a 10-percentage-point reduction in the probability of any utilisation). Among those who did seek care, treated patients also spent less ($4,800 vs $5,000 on average — a 4% reduction in cost per user). Both effects work in the same direction, but a drug payer considering formulary placement would want to know which effect is driving savings, because each implies a different care-management intervention."

Worked example

Scenario

A health outcomes analyst is comparing total annual allowed costs between a control arm and a treated arm in a retrospective cohort study. The analyst notices that some patients have zero cost in the study year — these are enrolled patients who never filed a relevant claim. The analyst applies a two-part model: Part 1 estimates the probability of any cost per arm; Part 2 estimates the mean cost among users per arm. Multiplying the two parts gives the overall expected cost per arm, and the difference is the estimated treatment contrast.

Dataset

Total annual allowed cost (USD) for nine patients across two arms. Zero-cost patients are structural non-users enrolled throughout the year who filed no qualifying claims.

Steps

• Count zero-cost patients in the control arm: 2 out of 5 have zero cost, so 3 out of 5 have positive cost. Part 1 probability (control): 3/5 = 0.6.

• Count zero-cost patients in the treated arm: 2 out of 4 have zero cost, so 2 out of 4 have positive cost. Part 1 probability (treated): 2/4 = 0.5.

• Compute mean cost among users in the control arm. Sum of user costs = 5000 + 5000 + 5000 = 15000. Mean among users (control) = 15000 / 3 = 5000.

• Compute mean cost among users in the treated arm. Sum of user costs = 4800 + 4800 = 9600. Mean among users (treated) = 9600 / 2 = 4800.

• Multiply Part 1 by Part 2 to get the overall expected cost. Control arm: 0.6 5000 = 3000. Treated arm: 0.5 4800 = 2400.

• The difference in overall expected cost is 3000 - 2400 = 600. Treated patients cost on average $600 less per year than control patients. Two components drive this: fewer treated patients used care at all (50% vs 60%) and those who did use care spent slightly less ($4800 vs $5000 per user).

Result

Control P(any cost) = 3/5 = 0.6, mean among users = 15000/3 = 5000, overall mean = 0.6 5000 = 3000. Treated P(any cost) = 2/4 = 0.5, mean among users = 9600/2 = 4800, overall mean = 0.5 4800 = 2400. Difference = 3000 - 2400 = 600. A sequential Shapley decomposition shows the utilisation component (Part 1) accounts for approximately $490 of the reduction and the conditional-cost component (Part 2) accounts for approximately $110; they do not contribute equally because the larger absolute gap is in the probability of use, not the cost per user.

Runnable example

import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf

# ── Simulate semicontinuous cost data ───────────────────────────────────────────────────────
rng = np.random.default_rng(seed=42)
n = 600
treatment = np.repeat([0, 1], n // 2)       # 300 control, 300 treated
age       = rng.normal(58, 12, n)

# Part 1 truth: logit(P(cost > 0)) = 0.5 + 0.3*treatment + 0.02*age
p_pos  = 1 / (1 + np.exp(-(0.5 + 0.3 * treatment + 0.02 * age)))
has_cost = rng.binomial(1, p_pos)

# Part 2 truth (among users): log(E[cost | cost>0]) = 7.5 - 0.2*treatment + 0.01*age
mu_pos = np.exp(7.5 - 0.2 * treatment + 0.01 * age)
cost   = np.where(has_cost, rng.gamma(shape=2.0, scale=mu_pos / 2.0), 0.0)

df = {"cost": cost, "treatment": treatment, "age": age,
      "any_cost": has_cost.astype(float)}
import pandas as pd
df = pd.DataFrame(df)

# ── Part 1: logistic regression for P(any cost) ─────────────────────────────────────────────
m1 = smf.logit("any_cost ~ treatment + age", data=df).fit(disp=0)

# ── Part 2: gamma GLM with log link on positive-cost patients ────────────────────────────────
pos = df.loc[df["cost"] > 0].copy()
m2  = smf.glm("cost ~ treatment + age", data=pos,
              family=sm.families.Gamma(link=sm.families.links.Log())).fit()

print("=== Part 1: Logistic for P(any cost) ===")
print(m1.summary().tables[1])
print("\n=== Part 2: Gamma GLM for E[cost | cost>0] ===")
print(m2.summary().tables[1])

# ── Combined expected cost via recycled prediction ───────────────────────────────────────────
def expected_cost(data, t_val):
    """Marginalised E[cost] under treatment assignment t_val."""
    d = data.assign(treatment=t_val)
    p_any = m1.predict(d)                        # P(cost > 0)
    mu    = m2.predict(d)                        # E[cost | cost > 0]
    return float((p_any * mu).mean())

e_control = expected_cost(df, 0)
e_treated = expected_cost(df, 1)
contrast  = e_control - e_treated
print(f"\nE[cost | control] = ${e_control:,.0f}")
print(f"E[cost | treated] = ${e_treated:,.0f}")
print(f"Mean cost reduction (control - treated) = ${contrast:,.0f}")

# ── Bootstrap 95% CI for the combined contrast ───────────────────────────────────────────────
B = 500
boot = np.empty(B)
for b in range(B):
    idx   = rng.choice(len(df), size=len(df), replace=True)
    df_b  = df.iloc[idx].reset_index(drop=True)
    pos_b = df_b.loc[df_b["cost"] > 0]
    try:
        b1 = smf.logit("any_cost ~ treatment + age", data=df_b).fit(disp=0)
        b2 = smf.glm("cost ~ treatment + age", data=pos_b,
                     family=sm.families.Gamma(link=sm.families.links.Log())).fit(disp=0)
        ec0 = (b1.predict(df_b.assign(treatment=0)) * b2.predict(df_b.assign(treatment=0))).mean()
        ec1 = (b1.predict(df_b.assign(treatment=1)) * b2.predict(df_b.assign(treatment=1))).mean()
        boot[b] = ec0 - ec1
    except Exception:
        boot[b] = np.nan
boot = boot[~np.isnan(boot)]
lo, hi = np.percentile(boot, [2.5, 97.5])
print(f"Bootstrap 95% CI for mean cost reduction: [${lo:,.0f}, ${hi:,.0f}]")

set.seed(42)
n         <- 600
treatment <- rep(0:1, each = n / 2)
age       <- rnorm(n, mean = 58, sd = 12)

# Part 1 truth: logit(P(cost > 0)) = 0.5 + 0.3*treatment + 0.02*age
p_pos    <- plogis(0.5 + 0.3 * treatment + 0.02 * age)
has_cost <- rbinom(n, 1, p_pos)

# Part 2 truth: log(E[cost | cost>0]) = 7.5 - 0.2*treatment + 0.01*age
mu_pos <- exp(7.5 - 0.2 * treatment + 0.01 * age)
cost   <- ifelse(has_cost == 1, rgamma(n, shape = 2, scale = mu_pos / 2), 0)

df <- data.frame(cost = cost, treatment = treatment, age = age,
                 any_cost = as.numeric(has_cost))

# ── Part 1: logistic regression for P(any cost) ────────────────────────────────────────────
m1 <- glm(any_cost ~ treatment + age, family = binomial, data = df)
cat("=== Part 1: Logistic for P(any cost) ===\n")
print(summary(m1)$coefficients)

# ── Part 2: gamma GLM with log link on positive-cost subset ───────────────────────────────
pos <- df[df$cost > 0, ]
m2  <- glm(cost ~ treatment + age, family = Gamma(link = "log"), data = pos)
cat("\n=== Part 2: Gamma GLM for E[cost | cost>0] ===\n")
print(summary(m2)$coefficients)

# ── Combined expected cost via recycled prediction ─────────────────────────────────────────
expected_cost <- function(data, t_val) {
  d <- data; d$treatment <- t_val
  p_any <- predict(m1, newdata = d, type = "response")
  mu    <- predict(m2, newdata = d, type = "response")
  mean(p_any * mu)
}
e_ctrl <- expected_cost(df, 0)
e_trt  <- expected_cost(df, 1)
cat(sprintf("\nE[cost | control] = $%.0f\n", e_ctrl))
cat(sprintf("E[cost | treated] = $%.0f\n", e_trt))
cat(sprintf("Mean cost reduction = $%.0f\n", e_ctrl - e_trt))

# ── Bootstrap 95% CI ──────────────────────────────────────────────────────────────────────
B    <- 500
boot <- numeric(B)
for (b in seq_len(B)) {
  idx  <- sample(nrow(df), replace = TRUE)
  db   <- df[idx, ]; posb <- db[db$cost > 0, ]
  tryCatch({
    b1   <- glm(any_cost ~ treatment + age, family = binomial, data = db)
    b2   <- glm(cost ~ treatment + age, family = Gamma(link = "log"), data = posb)
    ec0  <- mean(predict(b1, newdata = transform(db, treatment = 0), type = "response") *
                 predict(b2, newdata = transform(db, treatment = 0), type = "response"))
    ec1  <- mean(predict(b1, newdata = transform(db, treatment = 1), type = "response") *
                 predict(b2, newdata = transform(db, treatment = 1), type = "response"))
    boot[b] <- ec0 - ec1
  }, error = function(e) { boot[b] <<- NA_real_ })
}
boot <- boot[!is.na(boot)]
ci   <- quantile(boot, c(0.025, 0.975))
cat(sprintf("Bootstrap 95%% CI: [$%.0f, $%.0f]\n", ci[1], ci[2]))

# ── Tweedie GLM as one-model sensitivity check (requires tweedie package) ─────────────────
# library(tweedie)
# p_est <- tweedie.profile(cost ~ treatment + age, data = df, p.vec = seq(1.1, 1.9, 0.1))$p.max
# m_tw  <- glm(cost ~ treatment + age, data = df, family = tweedie(var.power = p_est, link.power = 0))
# cat(sprintf("Tweedie combined treatment coef: %.3f (exp = %.3f)\n",
#             coef(m_tw)["treatment"], exp(coef(m_tw)["treatment"])))