Gamma Distribution for Cost and Skewed Outcomes

In plain language

The gamma distribution is a statistical model for positive, right-skewed numbers like healthcare costs, where the spread of costs grows proportionally as the average cost grows. Fitting a gamma generalized linear model (GLM) with a log link lets an analyst estimate how patient characteristics or treatment multiply the average expected cost, producing results directly on the dollar scale — no mathematical back-transformation required. It is the modern standard for healthcare cost regression in real-world evidence studies, replacing the older practice of log-transforming costs and running ordinary regression, which recovers the geometric mean rather than the arithmetic mean that budget models and payers actually need. One important limitation: the gamma model does not accommodate patients with zero cost, so datasets with many zero-cost patients need a two-part model instead.

What the gamma distribution is and why cost data fits it

The gamma distribution is a two-parameter continuous probability distribution defined on strictly positive real values (0, ∞). Parameterised by a shape parameter α (alpha) and a scale parameter θ (theta), its mean is μ = αθ and its variance is Var(Y) = αθ² = μ²/α. This last expression is the key: the variance grows as the square of the mean. Equivalently, the coefficient of variation CV = SD/μ = 1/√α is a constant that does not depend on the mean level. As average costs rise — because sicker patients with more comorbidities drive both higher spending and proportionally higher spread — the variance rises with them. This is precisely the structural pattern that empirical healthcare cost data shows in claims, EHR, and registry datasets. Ordinary least squares (OLS) imposes none of this; the gamma distribution captures it by design.

The gamma family includes several familiar special cases. When α = 1 the gamma reduces to the exponential distribution. The chi-squared distribution with k degrees of freedom is a gamma with α = k/2 and θ = 2. The inverse-Gaussian (Wald) distribution is a related model with variance growing as μ³ — appropriate for very heavy-tailed costs when the Modified Park test (see below) returns λ̂ near 3 rather than 2.

Why OLS fails for cost data

Ordinary least squares makes two structural assumptions that healthcare cost distributions routinely violate. First, OLS assumes homoscedastic residuals — that the residual variance is constant across the range of fitted values. In cost data the opposite is true: patients with high predicted costs show far greater spread than low-cost patients, a textbook case of heteroscedasticity. Second, OLS imposes no non-negativity constraint; for patients with favourable covariates, a fitted OLS model can predict negative costs, which is not a meaningful quantity. Third, the heavy right tail of cost distributions (a small fraction of catastrophically expensive patients) inflates the OLS estimate of the mean and renders confidence intervals unreliable at any sample size.

The most common fix — logging the outcome and running OLS on log(cost) — removes heteroscedasticity and skewness from the model but creates a new problem: log-OLS estimates E[log(Y)|X], not E[Y|X]. Back-transforming by exponentiating the fitted values recovers the geometric mean, not the arithmetic mean. Because arithmetic means aggregate to population totals — which is exactly what a payer budget-impact model or HTA submission requires — the geometric mean is the wrong estimand. Duan's (1983) smearing estimator applies a nonparametric correction factor to recover the arithmetic mean, but this works correctly only when the log-scale residuals are themselves homoscedastic, which can also fail. The gamma GLM eliminates both the retransformation step and the residual-heteroscedasticity risk.

The gamma GLM with log link as the modern default

A GLM with a gamma variance function and a log link models E[Y|X] directly as:

E[Y|X] = μ = exp(β₀ + β₁X₁ + … + βₚXₚ)

with variance function Var(Y|X) = φ · μ² where φ = 1/α is the dispersion parameter. The log link ensures that all fitted values are strictly positive regardless of the covariate values, and the mean is estimated on the original dollar scale without any retransformation step.

The gamma distributional assumption is used for efficiency (optimal weighting of observations) rather than as a hard requirement for consistency. Manning and Mullahy (2001) show that the gamma GLM estimator of E[Y|X] is consistent even when the variance function is misspecified, provided the mean function (log link) is correctly specified. Sandwich (robust) standard errors protect inference against variance misspecification. This robustness property makes the gamma GLM with log link a reliable workhorse rather than a fragile distributional bet.

Interpreting the output

Consider a gamma GLM with log link comparing an index treatment to a comparator on total annual healthcare cost, adjusting for age, sex, Charlson Comorbidity Index, and baseline cost. The model returns a treatment coefficient of 0.405 with 95% CI [0.18, 0.63].

(1) Formal statistical interpretation. exp(0.405) = 1.50. This is the model-estimated ratio of the mean total cost for treated patients to the mean total cost for comparator patients, holding the four adjustment covariates fixed at their observed values. The confidence interval [exp(0.18), exp(0.63)] = [1.20, 1.88] means that values of the true mean cost ratio between 1.20 and 1.88 are compatible with the observed data at the conventional significance threshold — it does NOT mean there is a 95% probability that the true ratio lies in this interval. A frequentist CI is not a credible interval; it is a property of the procedure across hypothetical repeated samples.

(2) Practical interpretation for a decision-maker. "After adjusting for patient age, sex, comorbidity burden, and prior-year cost, treated patients' expected annual costs are approximately 50% higher than comparator patients' expected costs, with data consistent with 20% to 88% higher costs. This is an adjusted association from observational data; whether it reflects a causal cost increment attributable to the treatment depends on the extent to which the measured covariates capture the full confounding structure — in particular, residual confounding by disease severity or unmeasured indication factors may still be present."

Note what the coefficient does NOT tell you: it says nothing about median costs, it says nothing about absolute dollar differences (to get those, multiply the reference group's mean cost by 0.50), and it does not establish causality in an unadjusted observational design.

The Modified Park test for choosing the variance family

Manning and Mullahy (2001) proposed the Modified Park test to guide selection among variance families. The test fits a preliminary log-OLS model, computes the squared residuals on the log scale, regresses log(squared residuals) on log(fitted values from the OLS), and examines the slope coefficient λ̂:

• λ̂ ≈ 0: constant variance; OLS on levels is appropriate.
• λ̂ ≈ 1: variance proportional to μ; Poisson GLM or quasi-Poisson.
• λ̂ ≈ 2: variance proportional to μ²; gamma GLM is the natural choice.
• λ̂ ≈ 3: variance proportional to μ³; inverse-Gaussian GLM.

In practice, healthcare cost data in pharmacy and medical claims almost always returns λ̂ in the range 1.5–2.5, landing firmly in the gamma zone. The Poisson GLM with log link is a useful alternative when λ̂ is near 1 or when zeros are present (the Poisson distribution is defined at zero, making it more flexible than the gamma for zero-inclusive distributions when used as a quasi-likelihood estimator). The Modified Park test has limited power in small samples and should be combined with residual plots (plot raw residuals against fitted values) for a fuller distributional diagnosis.

The zero problem

The gamma distribution is defined on strictly positive values — it assigns zero probability to Y = 0. In real claims datasets, patients with zero total cost during a measurement window are common: newly enrolled members who used no services, patients with only well-visit care outside the study window, or patients whose costs fall below the deductible and are not reflected in claims. Applying a gamma GLM directly to a zero-inclusive distribution yields biased estimates because the zeros are structurally incompatible with the model's support.

The standard solution for zero-heavy cost outcomes is a two-part model (hurdle model): Part 1 fits a logistic regression for P(Y > 0 | X), and Part 2 fits a gamma GLM for E[Y | Y > 0, X]. The overall mean is E[Y|X] = P(Y > 0|X) × E[Y | Y > 0, X]. When only a small fraction of patients have zero costs (< 5–10%), the bias from ignoring zeros is modest and a gamma GLM with a small constant added (e.g., Y + 0.01) or a Poisson quasi-likelihood approach may be acceptable as a sensitivity analysis. For serious analysis of zero-heavy cost distributions, the two-part model is the principled choice and is discussed further in the healthcare-costs-pppm-pppy-pmpm concept.

Pros, cons, and trade-offs

Pros: - Directly estimates E[Y|X] on the natural dollar scale; exponentiated coefficients are mean cost ratios with no retransformation step, retransformation bias, or smearing correction. - The variance function V(μ) = φμ² matches the empirical pattern of claims-based cost data (rising spread with rising mean; approximately constant coefficient of variation). - Point estimates are consistent under mild variance misspecification when the log link is correctly specified (Manning & Mullahy 2001); sandwich SEs protect inference. - Accommodates covariate adjustment, inverse-probability-of-treatment weighting (IPTW), stratification, and marginal standardisation in a single GLM framework. - Implemented natively in statsmodels (Python), base R, and SAS PROC GENMOD — no specialty packages required. - Marginal standardisation (averaging model-predicted costs under counterfactual treatment assignments) is straightforward and produces absolute risk/cost differences as well as ratios.

Cons: - Does not handle zeros — requires a two-part model for cost distributions with a non-trivial zero mass. - Exponentiated coefficients are mean ratios, not mean differences in dollars; to report absolute cost increments, multiply by the reference group's marginal mean. - Standard asymptotic standard errors assume correct variance specification; use sandwich SEs as the default in observational cost analyses. - Gamma GLM can overfit in very small samples (< 30 per group); bootstrapped confidence intervals are advisable when n is small. - Modified Park test has low power at small n; lambda_hat estimates are noisy and should be combined with residual plots, not used as the sole criterion.

When to use

Use the gamma GLM with log link when: - The outcome is strictly positive continuous — per-patient total cost, medical charges, length of stay expressed in days when always > 0, per-member-per-month cost. - The distribution is right-skewed with variance that increases with the mean — the standard pattern for pharmacy and medical claims costs. - The target estimand is the mean cost (not the median), because means aggregate to population totals needed for budget-impact modelling, value-framework submissions, and payer negotiations. - Covariate adjustment for measured confounders, IPTW weighting from a propensity model, or standardisation across a reference population is required alongside the cost regression. - The Modified Park test returns λ̂ near 2, or a histogram/QQ plot of log-residuals shows approximately homoscedastic spread consistent with a gamma variance function. - The dataset has few or no zero-cost patients (< 5–10%); if zeros are common, a two-part model wrapper around the gamma GLM is the preferred architecture.

When NOT to use

• Zeros in the outcome: if 10% or more of patients have zero total cost in the study window,
• Bounded outcomes: proportions (0–1), adherence rates, and other outcomes bounded above
• Count outcomes: hospital admissions, emergency department visits, prescription fills, and
• When the estimand is median cost: if the policy question concerns the patient in the middle
• When causal language is required but confounding is uncontrolled: a gamma GLM on raw
• Small n with extreme outliers: in cohorts of fewer than 20 patients per group, a single
• Log-normal comparison: when a log-normal model fits better (the log-residuals are normally

Worked example

Scenario

A health outcomes analyst is studying total medical costs in a 12-month follow-up window for eight patients: four who received an index drug (treated) and four who received a comparator (untreated). The analyst wants to compute the mean cost for each group, confirm the mean ratio, and understand how a gamma GLM with log link represents this relationship so they can interpret the model coefficient correctly.

Dataset

Total 12-month medical cost (USD) per patient. The treated patients' values are exactly twice the untreated patients' values, making the mean ratio 2.0 and the coefficient of variation identical across groups — a clean illustration of the gamma CV property.

Steps

• Sum untreated costs: 1000+1500+2000+3500 = 8000. Mean untreated cost = 8000/4 = 2000.

• Sum treated costs: 2000+3000+4000+7000 = 16000. Mean treated cost = 16000/4 = 4000.

• Mean cost ratio (treated to untreated) = 4000/2000 = 2.0. Treated patients cost twice as much on average.

• In a gamma GLM with log link, the treatment coefficient represents the log of this mean ratio. The natural log of 2.0 is approximately 0.693, so the model returns a treatment coefficient near 0.693. Exponentiating confirms: exp(0.693) is approximately 2.0.

• Coefficient of variation check: for the untreated group, the deviations from the mean of 2000 are -1000, -500, 0, and +1500. The sample variance is ((-1000)^2+(-500)^2+0^2+1500^2)/(4-1) = (1000000+250000+0+2250000)/3 = 3500000/3, giving SD approximately 1080. CV = 1080/2000 approximately 0.54. For the treated group, all values are exactly double, so SD doubles to approximately 2160, mean doubles to 4000, and CV = 2160/4000 approximately 0.54. The CV is the same in both groups, exactly as the gamma variance function predicts.

• If the true coefficient were exactly ln(2) = 0.693, and if we reported it rounded to 0.405 in a real study (because of covariate adjustment shifting the estimate), then exp(0.405) = 1.50 — the model estimates treated patients cost about 50% more on average after covariate adjustment. This illustrates that the exp() of any log-link coefficient is always a mean cost ratio, never a mean cost difference in dollars.

Result

Mean untreated cost = 8000/4 = 2000. Mean treated cost = 16000/4 = 4000. Mean ratio = 4000/2000 = 2.0. In a gamma GLM with log link the treatment coefficient is approximately ln(2) near 0.693, and exp(0.693) is approximately 2.0, meaning treated patients' expected costs are twice as high. The coefficient of variation is approximately 0.54 in both groups, consistent with the gamma variance assumption that CV is constant across mean levels.

Runnable example

import numpy as np
import statsmodels.api as sm

# ── Simulate cost data: binary treatment (0=untreated, 1=treated) ──────────────────────────
rng = np.random.default_rng(seed=0)
n = 400
treatment = np.repeat([0, 1], n // 2)                       # 200 untreated, 200 treated
age       = rng.normal(60, 10, n)
# True model: log(E[cost]) = 7.0 + 0.405 * treatment + 0.02 * age
mu_true = np.exp(7.0 + 0.405 * treatment + 0.02 * age)
phi     = 0.5                                                # dispersion; shape = 1/phi = 2
cost    = rng.gamma(shape=1.0 / phi, scale=mu_true * phi)   # E[cost] = mu_true

# ── Design matrix ─────────────────────────────────────────────────────────────────────────
X = sm.add_constant(np.column_stack([treatment, age]))       # [intercept, treatment, age]

# ── Fit gamma GLM with log link ───────────────────────────────────────────────────────────
fit = sm.GLM(
    cost, X,
    family=sm.families.Gamma(link=sm.families.links.Log())
).fit(cov_type="HC3")                                        # sandwich SEs by default

print(fit.summary())

# ── Mean ratios: exp(coef) and 95% CI on ratio scale ─────────────────────────────────────
names = fit.model.exog_names
coefs = fit.params
cis   = fit.conf_int()
print("\nMean ratios (exp of coefficients):")
for nm, b, lo, hi in zip(names, coefs, cis.iloc[:, 0], cis.iloc[:, 1]):
    print(f"  {nm:12s}  ratio={np.exp(b):.3f}  95% CI [{np.exp(lo):.3f}, {np.exp(hi):.3f}]")

# ── Marginal standardisation: mean predicted cost under each treatment assignment ─────────
X_trt = sm.add_constant(np.column_stack([np.ones(n),  age]))
X_unt = sm.add_constant(np.column_stack([np.zeros(n), age]))
mean_trt = fit.predict(X_trt).mean()
mean_unt = fit.predict(X_unt).mean()
print(f"\nMarginalised mean cost (treated):   ${mean_trt:,.0f}")
print(f"Marginalised mean cost (untreated): ${mean_unt:,.0f}")
print(f"Marginalised mean ratio: {mean_trt / mean_unt:.3f}")
print(f"Marginalised mean difference: ${mean_trt - mean_unt:,.0f}")

# ── Modified Park test: guide variance family selection ──────────────────────────────────
# Step 1: preliminary log-OLS to get fitted values on the log scale
log_cost = np.log(cost)
b_ols, _, _, _ = np.linalg.lstsq(X, log_cost, rcond=None)
log_fitted = X @ b_ols
log_sq_resid = np.log((log_cost - log_fitted) ** 2)
valid = np.isfinite(log_sq_resid)
# Step 2: regress log(squared residuals) on log(fitted) -> slope = lambda_hat
park_X      = sm.add_constant(log_fitted[valid])
park_result = sm.OLS(log_sq_resid[valid], park_X).fit()
lambda_hat  = park_result.params[1]
print(f"\nModified Park test: lambda_hat = {lambda_hat:.2f}")
print("  ~0 -> OLS on levels  |  ~1 -> Poisson  |  ~2 -> Gamma  |  ~3 -> Inverse-Gaussian")

set.seed(0)
n         <- 400
treatment <- rep(0:1, each = n / 2)
age       <- rnorm(n, mean = 60, sd = 10)

# True model: log(E[cost]) = 7.0 + 0.405 * treatment + 0.02 * age
mu_true <- exp(7.0 + 0.405 * treatment + 0.02 * age)
phi     <- 0.5           # dispersion; shape = 1/phi = 2
cost    <- rgamma(n, shape = 1 / phi, scale = mu_true * phi)  # E[cost] = mu_true

df <- data.frame(cost = cost, treatment = treatment, age = age)

# ── Fit gamma GLM with log link ───────────────────────────────────────────────────────────
fit <- glm(cost ~ treatment + age, data = df, family = Gamma(link = "log"))
summary(fit)

# ── Mean ratios: exp(coef) and 95% CI on ratio scale ─────────────────────────────────────
cat("\nMean ratios (exp of coefficients):\n")
round(exp(coef(fit)), 3)
round(exp(confint(fit)), 3)   # profile-likelihood CIs; use sandwich SEs for robustness

# ── Marginal standardisation ──────────────────────────────────────────────────────────────
df_trt <- df; df_trt$treatment <- 1
df_unt <- df; df_unt$treatment <- 0
mean_trt <- mean(predict(fit, newdata = df_trt, type = "response"))
mean_unt <- mean(predict(fit, newdata = df_unt, type = "response"))
cat(sprintf("\nMarginalised mean cost (treated):   $%.0f\n", mean_trt))
cat(sprintf("Marginalised mean cost (untreated): $%.0f\n", mean_unt))
cat(sprintf("Marginalised mean ratio:      %.3f\n", mean_trt / mean_unt))
cat(sprintf("Marginalised mean difference: $%.0f\n", mean_trt - mean_unt))

# ── Modified Park test ────────────────────────────────────────────────────────────────────
log_cost   <- log(cost)
b_ols      <- coef(lm(log_cost ~ treatment + age, data = df))
log_fitted <- cbind(1, treatment, age) %*% b_ols
log_sq_res <- log((log_cost - as.numeric(log_fitted))^2)
valid      <- is.finite(log_sq_res)
park_fit   <- lm(log_sq_res[valid] ~ log_fitted[valid])
lambda_hat <- coef(park_fit)[2]
cat(sprintf("\nModified Park test: lambda_hat = %.2f\n", lambda_hat))
cat("  ~0 OLS  |  ~1 Poisson  |  ~2 Gamma  |  ~3 Inverse-Gaussian\n")