Quantile Regression

In plain language

Quantile regression is a statistical method that estimates how patient characteristics or treatment affect any chosen percentile of an outcome — such as the median cost or the 90th- percentile length of stay — rather than just the average. This matters enormously for healthcare costs, because a drug might look affordable for most patients (unchanged median) while quietly exploding the bills of the sickest 5% (a much higher 90th percentile). Unlike standard regression, quantile regression makes no assumption about the shape of the cost distribution and handles extreme high-cost patients without deleting or capping them — though it cannot directly replace a mean-cost estimate when budget impact is what the payer needs.

What quantile regression is and why the mean is sometimes the wrong target

Ordinary least squares (OLS) regression estimates the conditional mean — E[Y | X] — the average outcome for patients with a given covariate profile. For outcomes like total healthcare costs or hospital length of stay (LOS), the mean is dominated by a small fraction of patients with extreme values, and a single estimator on the mean may hide clinically and economically important heterogeneity across the distribution. Quantile regression, introduced by Koenker and Bassett (1978), solves this by estimating the conditional τ-th quantile — Q_τ(Y | X) — for any chosen probability level τ ∈ (0, 1). Setting τ = 0.5 yields median regression (also called LAD, least absolute deviations); τ = 0.90 yields the model for the 90th-percentile patient. Fitting quantile regressions at a grid of τ values traces out the entire shape of the covariate effect across the distribution.

The practical motivation in HEOR is stark. Suppose a specialty drug has no effect on the median annual cost — the typical patient's spending is unchanged — but dramatically increases the 95th- percentile cost because a subset of patients experience serious adverse events requiring intensive hospitalisation. An OLS regression on the mean picks up this tail shift and returns a significant positive mean cost difference, but gives no information about where in the distribution the effect is concentrated. A payer managing population risk cares whether a new therapy shifts 5% of patients into catastrophically expensive territory; a single mean estimate cannot answer that. Quantile regression at τ = 0.50 and τ = 0.90 (or 0.95) can.

The check-function intuition

Quantile regression minimises a weighted sum of absolute residuals — specifically, a loss function called the check function (or pinball loss) ρ_τ(u) that penalises positive residuals (predictions too low) by a factor τ and negative residuals (predictions too high) by a factor (1 − τ). At τ = 0.5 positive and negative errors are weighted equally, and minimising the total gives the conditional median — the value that splits predicted errors evenly above and below the fitted line. At τ = 0.90, an under-prediction is penalised nine times more heavily than an over-prediction, so the fitted line is pulled up to capture the 90th percentile. No Gaussian, gamma, or any other distributional assumption is imposed. The estimator is semi-parametric: it assumes a linear conditional quantile function but leaves the full conditional distribution unspecified. This makes it valid for the heavy-tailed, non-normal cost and LOS distributions typical of claims and EHR data.

Robustness without deleting data — contrast with trimming and winsorisation

Quantile regression at the median is inherently robust to outliers in the outcome. Because the check function weights residuals by their sign rather than squaring them, a single catastrophically expensive patient contributes one large positive residual with bounded weight τ, rather than a squared residual that can dominate the OLS objective. Crucially, this robustness is achieved without removing or winsorising the outlier. The patient remains in the dataset at their true cost; the estimator simply does not let that cost dominate the fit.

This is a meaningful methodological distinction from trimming and winsorisation. Trimming an extreme cost observation changes the estimand: the trimmed mean estimates average cost in a hypothetical population where tail events are excluded — not average cost in the actual population. Winsorisation caps the value, introducing deliberate downward bias in the estimated mean. Quantile regression at the median, by contrast, estimates the conditional median in the full population — a well-defined quantity that is not distorted by the magnitude of extreme values, because the median depends only on rank order, not on how far above it the most extreme observations lie. For upper quantiles (τ = 0.75, 0.90), the estimand is similarly the conditional quantile at that probability level in the full population — not in a tail-trimmed sub-population.

Inference: rank-based and bootstrap standard errors

Two types of standard errors are available for quantile regression estimates in practice:

• Rank-based standard errors (the Koenker-Bassett sandwich, computed by default in R's
• Bootstrap standard errors (nonparametric xy-pair bootstrap, which re-samples patient rows

Neither approach requires distributional assumptions beyond standard regularity conditions (continuity of the conditional distribution at the quantile, bounded moments). This is a stronger guarantee than OLS standard errors, which require homoscedasticity or the additional heteroscedasticity-robust sandwich correction.

When the mean is still the right target

Quantile regression is powerful, but it does not replace mean-based analysis when the decision problem requires means. Budget-impact models, cost-effectiveness analyses (ICER numerators), and payer negotiations about total drug budget all depend on the arithmetic mean cost per patient because means aggregate to population totals. If median drug cost is $2,000 but mean drug cost is $15,000, the payer's actual budget exposure is proportional to the mean. Median regression cannot be directly used to project budget impact — it answers a different question. For budget-impact and value-framework submissions, a gamma GLM with log link (see the gamma-distribution entry) remains the primary analytic tool. Quantile regression at the median is not interchangeable with a mean estimate; they are separate estimands.

Quantile regression is best positioned as: - Alongside a mean model, to characterise where effects concentrate and whether the mean summarises the distribution adequately. - As the primary analysis when the policy question directly concerns a specific quantile: "Does this intervention reduce catastrophic cost (i.e., the 90th-percentile cost)?" or "What is the median LOS for high-risk patients on this care pathway?" - As a robustness check when outlier sensitivity is a concern; if the median regression and the mean regression give qualitatively different answers, the difference is informative about the role of the tail.

Censored quantile regression

When outcomes are right-censored (costs or LOS censored at the end of the observation window, or time-to-event outcomes), censored quantile regression — implemented via the Powell (1986) estimator and available in quantreg::crq() in R — extends the framework to survival-type data, analogously to how the Cox model extends OLS regression to censored time-to-event outcomes.

Pros, cons, and trade-offs

Pros: - Estimates any conditional quantile of the outcome distribution; reveals heterogeneous covariate effects that a single mean estimate cannot capture — for example, a drug that reduces median cost while inflating 90th-percentile cost. - No distributional assumption on the outcome; valid for the heavy-tailed, zero-inclusive, non-normal distributions typical of healthcare cost and LOS data. - Robust to outliers in the outcome without deleting or winsorising them; the estimand remains the conditional quantile in the full population. - Fitting at multiple τ values (e.g., 0.25, 0.50, 0.75, 0.90, 0.95) provides a complete picture of distributional shift — richer than any single mean estimate. - Implemented natively in Python (statsmodels QuantReg), R (quantreg::rq), and SAS (PROC QUANTREG); no specialty packages required.

Cons: - Estimates conditional quantiles, not means; cannot substitute for a mean-based estimator when budget impact or ICER calculations are the goal. - Less efficient than OLS when the outcome is genuinely normally distributed; OLS achieves lower variance for the mean under normality, while quantile regression sacrifices some efficiency for distributional flexibility. - Covariate adjustment and propensity score weighting (weighted quantile regression) are feasible but require careful specification, and bootstrap standard errors are computationally intensive in large datasets. - Quantile crossing — fitted quantile curves for different τ values that cross each other at some covariate values, violating the mathematical monotonicity requirement — can occur in finite samples and requires joint quantile estimation or rearrangement to correct. - Interpretation of upper-quantile coefficients is less familiar to clinical and payer audiences than a mean cost ratio and requires explicit communication in the analysis report.

When to use

Use quantile regression when: - The research question concerns a specific quantile of the cost or LOS distribution — the median, the 90th percentile, or the probability of catastrophic cost above a threshold. - The treatment effect may be heterogeneous across the distribution: a drug may reduce median cost while concentrating high costs in a small subgroup, a pattern that OLS on the mean cannot identify. - The outcome distribution is sufficiently non-normal and heavy-tailed that distributional assumptions in GLMs are unattractive and winsorisation would materially change the estimand. - A distributional-robustness check is desired alongside a primary mean analysis; divergence between median and mean regression signals the presence of an influential tail. - The decision question is specifically about payer risk management at the catastrophic-cost tier — "What fraction of the benefit-risk trade-off occurs in the top 10% of spenders?"

When NOT to use

Do not use quantile regression as the primary analysis when: - The decision-relevant estimand is the arithmetic mean cost (budget impact, ICER numerator, PMPM premium calculation); the mean and the median are different quantities, and substituting one for the other is a methodological error for budget-relevant questions. Route instead to a gamma GLM with log link, a two-part model, or a bootstrap mean difference. - The sample is very small (n < 50 per group); quantile regression standard errors are unreliable at extreme quantiles (τ near 0 or 1) with limited data, and bootstrap confidence intervals can be wide and unstable. - Outcomes are binary or count-valued and an approximating continuous-quantile model is inappropriate; use logistic regression for binary outcomes and Poisson or negative-binomial models for counts. - The analysis requires causal language but confounding is uncontrolled; an unadjusted quantile regression on observational data estimates a biased conditional quantile. Pair with propensity score weighting (weighted quantile regression) or pre-specify as a descriptive sensitivity.

Interpreting the output

Consider a study comparing 12-month total healthcare costs for patients initiating Drug A versus Drug B in a commercial claims database, covariate-adjusted for age, sex, Charlson Comorbidity Index, and prior-year cost. Quantile regression is fit at τ = 0.50 and τ = 0.90. Results: the τ = 0.50 coefficient on Drug A vs Drug B is −500 (95% CI −900 to −100); the τ = 0.90 coefficient is +8,200 (95% CI +1,500 to +14,900).

(1) Formal interpretation. The τ = 0.50 coefficient of −500 estimates that, conditional on the four adjustment covariates held fixed at their observed values, the 50th percentile of the 12-month cost distribution is $500 lower for Drug A patients than for Drug B patients. This is a statement about the conditional median, not the mean cost difference. A 95% CI of [−900, −100] means that values of the true conditional median difference between −$900 and −$100 are compatible with the observed data under the repeated-sampling interpretation of a frequentist confidence interval — it does NOT mean there is 95% probability the true median difference falls in that range. The τ = 0.90 coefficient of +8,200 estimates that the 90th-percentile cost is $8,200 higher for Drug A patients, with a 95% CI of [+1,500, +14,900] excluding zero. These two estimates come from separate quantile-regression fits and can differ in sign — that divergence is the signal quantile regression is designed to detect.

(2) Practical interpretation. The typical Drug A patient (at the median) costs $500 less than the typical Drug B patient, but the highest-cost 10% of Drug A patients cost $8,200 more than their counterparts on Drug B. Drug A appears to compress costs for most patients while concentrating expense in a high-cost subgroup — possibly patients experiencing serious adverse events or disease progression requiring rescue therapy. A mean regression might show a modest positive or null mean cost difference that conceals this distributional story entirely. A payer managing population risk needs both numbers: the median tells them what the majority of members will experience; the 90th-percentile coefficient quantifies the tail-budget exposure. However, for a formal budget-impact model, the payer must use the arithmetic mean cost estimate — obtainable from a gamma GLM — not the median, because population budget = mean cost × member count.

Worked example

Scenario

A health economist is comparing 12-month costs for five patients in a claims database: four patients had routine care and one had a catastrophic complication requiring a transplant. The economist wants to understand the difference between the mean and median cost summary, and what happens to each when the catastrophic patient's cost doubles — mimicking how a more expensive rescue therapy would affect the statistics reported to a payer.

Dataset

Twelve-month total cost (thousands of USD) for five patients. Patient P5 is a catastrophic outlier; the other four represent typical utilization. All values are used in both the mean and median calculation.

Steps

• Sort the five costs in ascending order: 1, 2, 3, 4, 90.

• The median is the middle (3rd) value of the sorted list: median = 3 (thousand USD).

• Compute the mean: sum = 1+2+3+4+90 = 100. Mean = 100/5 = 20 (thousand USD).

• Now suppose a more expensive rescue therapy doubles P5's cost from 90 to 180 (thousand USD). The sorted list becomes 1, 2, 3, 4, 180.

• New mean: sum = 1+2+3+4+180 = 190. New mean = 190/5 = 38 (thousand USD). The mean nearly doubled, from 20 to 38.

• New median: the middle (3rd) value of [1, 2, 3, 4, 180] is still 3. The median is unchanged = 3 (thousand USD).

Result

Baseline: median = 3 (thousand USD), mean = 100/5 = 20. After doubling P5's cost: median = 3 (unchanged), mean = 190/5 = 38. Doubling the extreme value nearly doubled the mean but left the median completely unchanged. A mean regression detects the catastrophic-cost shift; a median regression does not. A drug that causes 1 in 5 patients to have twice the catastrophic cost appears to nearly double average spending — yet the typical (median) patient's cost is unaffected. Both facts are true and relevant: the payer's budget depends on the mean; the clinical experience of most patients is captured by the median.

Runnable example

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

# ── Simulate cost data: binary treatment, asymmetric distributional effect ──────────────────
rng = np.random.default_rng(seed=42)
n = 600
treatment = np.repeat([0, 1], n // 2)
age       = rng.normal(60, 10, n)
# True model: median cost lower for treated; upper tail higher for treated (heterogeneous effect)
noise = rng.exponential(scale=8_000, size=n)
cost  = (5_000
         + 100 * age
         - 600 * treatment          # reduces median
         + 12_000 * treatment * (noise > np.quantile(noise, 0.85))   # inflates 90th pct
         + noise)
df = pd.DataFrame({"cost": cost, "treatment": treatment, "age": age})

# ── 1. Median regression (tau = 0.50, rank-based SEs) ────────────────────────────────────
mod = smf.quantreg("cost ~ treatment + age", df)
res_50 = mod.fit(q=0.5)                       # default: rank-based standard errors
print("=== Median regression (tau=0.50) ===")
print(res_50.summary())
ci_50 = res_50.conf_int()
print(f"Treatment: coef={res_50.params['treatment']:,.0f}  "
      f"95% CI [{ci_50.loc['treatment', 0]:,.0f}, {ci_50.loc['treatment', 1]:,.0f}]")

# ── 2. 90th-percentile regression (tau = 0.90, bootstrap SEs) ────────────────────────────
res_90 = mod.fit(q=0.9, vcov="iid")           # use "robust" or bootstrap for extreme quantiles
print("\n=== 90th-percentile regression (tau=0.90) ===")
ci_90 = res_90.conf_int()
print(f"Treatment: coef={res_90.params['treatment']:,.0f}  "
      f"95% CI [{ci_90.loc['treatment', 0]:,.0f}, {ci_90.loc['treatment', 1]:,.0f}]")
print("\nNote: For tau>=0.90, prefer bootstrap SEs. statsmodels quantreg supports")
print("      vcov='robust' (sandwich) or bootstrap via res.conf_int(p=0.05, q=0.90).")

# ── 3. Quantile grid: trace treatment coefficient across tau ──────────────────────────────
print("\n=== Treatment coefficient across quantile grid ===")
print(f"{'tau':>6}  {'coef':>10}  {'CI_low':>10}  {'CI_high':>10}")
for tau in [0.10, 0.25, 0.50, 0.75, 0.90]:
    r   = mod.fit(q=tau, vcov="iid")
    ci  = r.conf_int()
    b   = r.params["treatment"]
    lo  = ci.loc["treatment", 0]
    hi  = ci.loc["treatment", 1]
    print(f"{tau:>6.2f}  {b:>10,.0f}  {lo:>10,.0f}  {hi:>10,.0f}")
print("\nA rising coefficient from tau=0.50 to tau=0.90 indicates the treatment")
print("concentrates its cost effect in the upper tail (heterogeneous distributional effect).")

library(quantreg)

# ── Simulate cost data (same structure as Python example) ────────────────────────────────
set.seed(42)
n         <- 600
treatment <- rep(0:1, each = n / 2)
age       <- rnorm(n, mean = 60, sd = 10)
noise     <- rexp(n, rate = 1 / 8000)
cost      <- (5000
              + 100 * age
              - 600 * treatment
              + 12000 * treatment * (noise > quantile(noise, 0.85))
              + noise)
df <- data.frame(cost = cost, treatment = treatment, age = age)

# ── 1. Median regression: tau = 0.50, rank-based standard errors ─────────────────────────
fit_50 <- rq(cost ~ treatment + age, data = df, tau = 0.50)
cat("=== Median regression (tau = 0.50) ===\n")
print(summary(fit_50, se = "rank"))          # rank-based SEs: default, reliable near median
# Coefficient interpretation: the treatment coefficient estimates the shift in the
# conditional MEDIAN cost (not the mean) for treated vs untreated patients with same age.

# ── 2. 90th-percentile regression: tau = 0.90, bootstrap SEs ────────────────────────────
fit_90 <- rq(cost ~ treatment + age, data = df, tau = 0.90)
cat("\n=== 90th-percentile regression (tau = 0.90) ===\n")
print(summary(fit_90, se = "boot", R = 1000, bsmethod = "xy"))
# se = "boot", bsmethod = "xy": nonparametric xy-pair bootstrap; preferred for tau >= 0.90

# ── 3. Quantile grid: fit all tau values at once and plot ─────────────────────────────────
taus    <- c(0.10, 0.25, 0.50, 0.75, 0.90)
fit_all <- rq(cost ~ treatment + age, data = df, tau = taus)
cat("\n=== Treatment coefficient across quantile grid ===\n")
grid_summary <- summary(fit_all, se = "rank")
for (i in seq_along(taus)) {
  tbl <- grid_summary[[i]]$coefficients
  b   <- tbl["treatment", "Value"]
  lo  <- tbl["treatment", "Lower bd"]
  hi  <- tbl["treatment", "Upper bd"]
  cat(sprintf("tau=%.2f  coef=%8.0f  95%%CI [%8.0f, %8.0f]\n", taus[i], b, lo, hi))
}

# ── 4. Censored quantile regression for right-censored LOS outcomes ───────────────────────
# crq() in quantreg implements the Powell (1986) censored quantile estimator.
# Example structure (replace LOS and censor_flag with real columns):
# fit_crq <- crq(Surv(LOS, 1 - censor_flag, type = "right") ~ treatment + age,
#                data = df, method = "Powell", taus = 0.50)
# summary(fit_crq, taus = 0.50)