Bayesian quantile regression, and when it lies

#The setting

Birth weight is one of the most studied numbers in medicine, and almost nobody cares about its average. A baby of 3.3 kg and one of 3.5 kg are both fine. What matters are the tails: babies under 2,500 g ("low birth weight", about 7% of US births) carry most of the risk of death and disability, and babies over 4,000 g ("macrosomia") bring their own problems at delivery. A clinician's questions are about quantiles: does smoking shift the whole distribution down, or mainly the lower tail? Does diabetes before pregnancy make babies bigger, or just more variable?

Ordinary regression answers a different question - how the mean moves - and assumes every quantile moves by the same amount. Quantile regression (Koenker & Bassett 1978) fits each quantile directly. Koenker & Hallock (2001, Journal of Economic Perspectives) made US birth weights the textbook example; we redo it on 2023 births.

Doing it the Bayesian way looks easy: there is a distribution, the asymmetric Laplace, whose likelihood is maximised exactly at the quantile-regression estimate, and PyMC has it. It is used in many papers and packages (bayesQR, brms's asym_laplace). But nobody believes birth weights are asymmetric-Laplace distributed - it is a working likelihood - and a posterior built on a wrong likelihood can put its centre in the right place and its width in the wrong one. This notebook measures how wrong, fixes it two ways, and then asks whether we should model the whole distribution instead.

part	question	tool
A	What does a quantile minimise?	the check loss, linear programming, a Frisch-Newton solver
B	Do the effects differ across the distribution?	classical quantile regression on 2023 births, bootstrap, the quantile process
C	Can a likelihood give us quantile regression?	the asymmetric Laplace in PyMC; mode vs spread; the units trap
D	How wrong are its credible intervals?	a simulation with known truth, coverage, theory
E	How do we fix them?	sandwich (Yang-Wang-He) adjustment; a Gibbs posterior with a bootstrap-calibrated learning rate
F	What if the fitted quantiles cross?	birth-weight-for-gestational-age centiles; rearrangement; monotone increments
G	Or model the whole distribution?	sinh-arcsinh regression; derived quantile effects; held-out pinball loss, interval score, PIT

E41 built growth charts from a distributional model (LMS / Box-Cox t); here we start from the quantiles and see what that approach can and cannot give.

import logging
import time
import warnings

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import pytensor.tensor as pt
import scipy.sparse as sparse
from scipy import integrate, special, stats
from scipy.interpolate import BSpline
from scipy.optimize import linprog

from pymc_challenges import data

RANDOM_SEED = 83
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-variat")
logging.getLogger("pymc").setLevel(logging.WARNING)
warnings.filterwarnings("ignore", category=RuntimeWarning, module="arviz")
warnings.filterwarnings("ignore", message=".*sample_stats.*")
pd.set_option("display.width", 170)
pd.set_option("display.max_columns", 20)
T_START = time.time()

BLUE, ORANGE, AQUA, PURPLE = "#2a78d6", "#eb6834", "#1baf7a", "#8e5bd0"
INK, MUTED, LIGHT = "#0b0b0b", "#8a8984", "#d9d8d3"


def check_loss(r, tau):
    """The check (pinball) loss rho_tau(r) = r (tau - 1{r < 0}), elementwise."""
    return r * (tau - (r < 0))


def _step_bound(x, dx):
    out = np.full_like(x, 1e20)
    neg = dx < 0
    out[neg] = -x[neg] / dx[neg]
    return out


def rq_fn(X, y, tau, tol=1e-10, max_it=60):
    """Quantile regression by the Frisch-Newton interior-point method (Portnoy & Koenker 1997).

    A NumPy port of Morillo & Koenker's `lp_fnm`: it solves the DUAL linear programme
    max y'a  s.t.  X'a = (1 - tau) X'1,  0 <= a <= 1   and returns the primal coefficients."""
    n, p = X.shape
    A, c, u = X.T, -y, np.ones(n)
    x = (1 - tau) * u
    b = A @ x
    s = u - x
    yy = np.linalg.lstsq(A.T, c, rcond=None)[0]
    r = c - A.T @ yy
    r = r + 0.001 * (r == 0)
    z = np.where(r > 0, r, 0.0)
    w = z - r
    gap, scale = c @ x - yy @ b + w @ u, np.abs(c).sum() + 1.0
    for _ in range(max_it):
        if gap <= tol * scale:
            break
        q = 1 / (z / x + w / s)
        r = z - w
        AQ = A * np.sqrt(q)
        M = AQ @ AQ.T                                   # p x p normal equations
        rhs = np.sqrt(q) * r
        dy = np.linalg.solve(M, AQ @ rhs)
        dx = q * (A.T @ dy - r)
        ds, dz, dw = -dx, -z * (dx / x + 1), -w * (-dx / s + 1)
        fp = min(0.9995 * min(_step_bound(x, dx).min(), _step_bound(s, ds).min()), 1)
        fd = min(0.9995 * min(_step_bound(w, dw).min(), _step_bound(z, dz).min()), 1)
        if min(fp, fd) < 1:                             # Mehrotra-type corrector step
            mu = z @ x + w @ s
            g = (z + fd * dz) @ (x + fp * dx) + (w + fd * dw) @ (s + fp * ds)
            mu = mu * (g / mu) ** 3 / (2 * n)
            dxdz, dsdw = dx * dz, ds * dw
            xi = mu * (1 / x - 1 / s)
            rhs = rhs + np.sqrt(q) * (dxdz - dsdw - xi)
            dy = np.linalg.solve(M, AQ @ rhs)
            dx = q * (A.T @ dy + xi - r - dxdz + dsdw)
            ds = -dx
            dz = mu / x - z - z * dx / x - dxdz
            dw = mu / s - w - w * ds / s - dsdw
            fp = min(0.9995 * min(_step_bound(x, dx).min(), _step_bound(s, ds).min()), 1)
            fd = min(0.9995 * min(_step_bound(w, dw).min(), _step_bound(z, dz).min()), 1)
        x, s = x + fp * dx, s + fp * ds
        yy, w, z = yy + fd * dy, w + fd * dw, z + fd * dz
        gap = c @ x - yy @ b + w @ u
    return -yy


def rq_linprog(X, y, tau):
    """The same problem as a primal LP for scipy's HiGHS: min tau 1'u + (1 - tau) 1'v, X b + u - v = y."""
    n, p = X.shape
    A = sparse.hstack([sparse.csr_matrix(X), sparse.eye(n), -sparse.eye(n)]).tocsc()
    cost = np.r_[np.zeros(p), np.full(n, tau), np.full(n, 1 - tau)]
    res = linprog(cost, A_eq=A, b_eq=y, bounds=[(None, None)] * p + [(0, None)] * (2 * n), method="highs")
    return res.x[:p]


def pinball(y, q, tau):
    return np.mean(check_loss(y - q, tau))

#A. What does a quantile minimise?

The mean is the number $m$ that minimises the average squared distance $\frac1n\sum_i (y_i - m)^2$; the median minimises the average absolute distance. The $\tau$-quantile minimises an asymmetric absolute distance, the check or pinball loss

$$\rho_\tau(r) = r\,\bigl(\tau - \mathbb 1\{r < 0\}\bigr) = \begin{cases} \tau\, r & r \ge 0 \\ (\tau - 1)\, r & r < 0.\end{cases}$$

For $\tau = 0.1$ an observation below the candidate costs 0.9 per gram and one above costs 0.1, so the best candidate sits where 10% of the data are below it: moving it up by a small step gains $0.9 \times$ (share below) and loses $0.1 \times$ (share above), which balance exactly at share below $= 0.1$. Quantile regression replaces the constant by a linear predictor:

$$\hat\beta(\tau) = \arg\min_\beta \sum_i \rho_\tau(y_i - x_i^\top\beta).$$

The objective is piecewise linear, so this is a linear programme (write each residual as $u_i - v_i$ with $u_i, v_i \ge 0$). The solution passes exactly through $p$ of the data points, and no derivative of the loss exists there - a first hint that likelihood-style curvature arguments will need care.

The data first. load only parses the raw codes; unknown values are coded 9, 99, 99.9 or 9999.

data.describe("natality_births")
raw = data.load("natality_births")
print(f"{len(raw):,} singleton births; birth weight unknown for {np.sum(raw.dbwt == 9999)}")
early = raw.oegest_comb < 28
print(f"births before 28 weeks: {early.sum()}; of these, 3rd-trimester cigarettes 'unknown' (99): "
      f"{np.mean(raw.cig_3[early] == 99):.0%}; 1st-trimester cigarettes unknown: {np.mean(raw.cig_1[early] == 99):.0%}")

natality_births: A simple random sample of 20,000 of the 3,491,735 SINGLETON births of 2023 (raw codes; 9 / 99 / 99.9 / 999 / 9999 = unknown, empty = not reported): dob_mm, mager (mother's age), mrace6 (1 White, 2 Black, 3 AIAN, 4 Asian, 5 NHOPI, 6 more than one race), dmar (1 married, 2 unmarried; not reported for California), meduc (1-8 education, 9 unknown), tbo_rec (total birth order), previs (prenatal visits), cig_0..cig_3 (cigarettes per day before pregnancy and in trimesters 1-3), m_ht_in (inches), bmi (pre-pregnancy), pwgt_r (pre-pregnancy weight, lb), wtgain (lb), rf_pdiab / rf_gdiab / rf_ghype (pre-pregnancy diabetes, gestational diabetes, gestational hypertension: Y/N/U), sex (M/F), combgest and oegest_comb (gestation in weeks, combined and obstetric estimate), dbwt (birth weight, grams).
  source: National Center for Health Statistics (NCHS), National Vital Statistics System, Natality public use file 2023 (all 3,605,081 US birth certificates of 2023; codes as in the NCHS 'User Guide to the 2023 Natality Public Use File'), read from NBER's CSV conversion (the URL; 188 MB). US government work, public domain; NCHS asks users to cite the source and not to attempt to identify anyone. BUILT by tools/build_e83_natality.py (a seeded random subsample) - keep the cached file.
  url:    https://data.nber.org/nvss/natality/csv/2023/natality2023us.zip
20,000 singleton births; birth weight unknown for 17
births before 28 weeks: 107; of these, 3rd-trimester cigarettes 'unknown' (99): 76%; 1st-trimester cigarettes unknown: 1%

A trap worth a sentence: "cigarettes in the 3rd trimester" is unknown for three quarters of the babies born before 28 weeks - there was no third trimester. Dropping incomplete records on that variable would silently delete the smallest babies, i.e. the lower tail we care about. We define smoker as any smoking in any trimester (first or second trimester must be known), and drop the few records with other unknown items.

ok = ((raw.dbwt < 9999) & (raw.cig_1 < 99) & (raw.cig_2 < 99) & (raw.bmi < 99) & (raw.wtgain < 99)
      & (raw.rf_pdiab != "U") & (raw.meduc < 9) & (raw.tbo_rec < 9) & (raw.oegest_comb < 99))
d = raw[ok].reset_index(drop=True).copy()
d["bw"] = d.dbwt / 1000                                     # kilograms
d["boy"] = (d.sex == "M").astype(float)
d["black"] = (d.mrace6 == 2).astype(float)
d["age"] = (d.mager - 30) / 10                              # per decade, 0 = 30 years
d["first"] = (d.tbo_rec == 1).astype(float)
d["college"] = (d.meduc >= 6).astype(float)                 # bachelor's degree or more
d["smoker"] = ((d.cig_1 > 0) | (d.cig_2 > 0) | ((d.cig_3 > 0) & (d.cig_3 < 99))).astype(float)
d["bmi5"] = (d.bmi - 25) / 5                                # per 5 BMI units, 0 = 25
d["gain10"] = (d.wtgain - 30) / 10                          # per 10 lb gained, 0 = 30 lb
d["pre_diab"] = (d.rf_pdiab == "Y").astype(float)
d["gest_hyp"] = (d.rf_ghype == "Y").astype(float)
COVS = ["boy", "black", "age", "first", "college", "smoker", "bmi5", "gain10", "pre_diab", "gest_hyp"]
COEF = ["intercept"] + COVS
N_TRAIN = 6000
perm = rng.permutation(len(d))
tr, te = d.iloc[perm[:N_TRAIN]].reset_index(drop=True), d.iloc[perm[N_TRAIN:]].reset_index(drop=True)
X, y = np.column_stack([np.ones(N_TRAIN)] + [tr[c].to_numpy() for c in COVS]), tr.bw.to_numpy()
Xte, yte = np.column_stack([np.ones(len(te))] + [te[c].to_numpy() for c in COVS]), te.bw.to_numpy()
print(f"{len(d):,} complete records ({ok.mean():.1%}); training {len(tr):,}, held out {len(te):,}")
print(f"low birth weight (< 2.5 kg): {np.mean(d.bw < 2.5):.1%}; macrosomia (> 4 kg): {np.mean(d.bw > 4):.1%}")
print(pd.DataFrame({"share / mean": tr[COVS].mean(), "count (binary)": tr[COVS].sum()}).round(3).T.to_string())

18,936 complete records (94.7%); training 6,000, held out 12,936
low birth weight (< 2.5 kg): 6.9%; macrosomia (> 4 kg): 7.1%
                     boy    black      age     first   college   smoker      bmi5   gain10  pre_diab  gest_hyp
share / mean       0.511    0.153   -0.051     0.328     0.364    0.031     0.542   -0.077     0.013     0.096
count (binary)  3068.000  918.000 -308.200  1971.000  2182.000  186.000  3254.040 -461.400    76.000   576.000

Before any regression: the loss curves for three quantile levels on the 6,000 training weights, and a check that the minimiser is the sample quantile.

cand = np.linspace(0.5, 5.0, 901)
TAU3 = (0.1, 0.5, 0.9)
fig, axes = plt.subplots(1, 3, figsize=(14, 4))
ax = axes[0]
r_ = np.linspace(-2, 2, 201)
for tau, c in zip(TAU3, (ORANGE, INK, BLUE)):
    ax.plot(r_, check_loss(r_, tau), color=c, lw=2, label=f"tau = {tau}")
ax.set(xlabel="residual r = y - q", ylabel="loss", title="the check (pinball) loss")
ax.legend()
ax = axes[1]
for tau, c in zip(TAU3, (ORANGE, INK, BLUE)):
    loss = np.array([pinball(y, q, tau) for q in cand])
    ax.plot(cand, loss, color=c, lw=2)
    ax.axvline(cand[np.argmin(loss)], color=c, ls=":", lw=1)
    print(f"tau = {tau}: argmin of the average check loss {cand[np.argmin(loss)]:.3f} kg, "
          f"sample quantile {np.quantile(y, tau):.3f} kg")
ax.set(xlabel="candidate value q (kg)", ylabel="average check loss", xlim=(1.5, 5),
       title="average loss over 6,000 birth weights")
ax = axes[2]
ax.hist(y, bins=np.arange(0.2, 6.0, 0.1), density=True, color=LIGHT, edgecolor="white")
xx = np.linspace(0.2, 6, 400)
ax.plot(xx, stats.norm.pdf(xx, y.mean(), y.std()), color=INK, lw=1.5, label="normal, same mean and sd")
for tau, c in zip(TAU3, (ORANGE, INK, BLUE)):
    ax.axvline(np.quantile(y, tau), color=c, ls=":", lw=1)
ax.axvline(2.5, color=ORANGE, lw=1.2, label="2.5 kg (low birth weight)")
ax.set(xlabel="birth weight (kg)", ylabel="density", title="training birth weights")
ax.legend(fontsize=8)
ax_in = ax.inset_axes([0.62, 0.55, 0.35, 0.3])
ax_in.hist(y, bins=np.arange(0.2, 6.0, 0.1), density=True, color=LIGHT, edgecolor="white")
ax_in.plot(xx, stats.norm.pdf(xx, y.mean(), y.std()), color=INK, lw=1)
ax_in.set(xlim=(0.3, 2.3), ylim=(0, 0.12), yticks=[])
ax_in.tick_params(labelsize=7)
fig.suptitle("A quantile is the minimiser of the check loss (dotted: 10%, 50%, 90% quantiles)", fontsize=11);

tau = 0.1: argmin of the average check loss 2.640 kg, sample quantile 2.637 kg
tau = 0.5: argmin of the average check loss 3.300 kg, sample quantile 3.300 kg
tau = 0.9: argmin of the average check loss 3.915 kg, sample quantile 3.912 kg

The minimisers agree with the sample quantiles to the grid spacing (5 g). The histogram shows why tails need their own model: birth weight has a long left tail (the inset: babies born early), far heavier than a normal distribution with the same mean and sd, and a shorter right tail.

Now the regression. We solve the linear programme two ways: with scipy's general LP solver (HiGHS) and with a short Frisch-Newton interior-point method written above (Portnoy & Koenker 1997; the "fn" method of R's quantreg, recommended there for large problems). Interior-point methods walk through the inside of the feasible set; each step is one weighted least-squares solve with a $p \times p$ system, and a few dozen steps suffice whatever the sample size. We will need thousands of fits for bootstraps, so speed matters.

for tau in (0.1, 0.9):
    t0 = time.time()
    b_lp = rq_linprog(X, y, tau)
    t_lp = time.time() - t0
    t0 = time.time()
    b_fn = rq_fn(X, y, tau)
    t_fn = time.time() - t0
    print(f"tau = {tau}: HiGHS {t_lp * 1000:.0f} ms, Frisch-Newton {t_fn * 1000:.0f} ms; largest coefficient "
          f"difference {np.abs(b_lp - b_fn).max() * 1000:.4f} g; check loss difference "
          f"{(check_loss(y - X @ b_fn, tau).sum() - check_loss(y - X @ b_lp, tau).sum()) * 1000:.2g} g")

tau = 0.1: HiGHS 678 ms, Frisch-Newton 13 ms; largest coefficient difference 0.0000 g; check loss difference 2.3e-07 g

tau = 0.9: HiGHS 673 ms, Frisch-Newton 13 ms; largest coefficient difference 0.0000 g; check loss difference 1.3e-07 g

Same answer (coefficients agree to better than 0.0001 g), roughly 50 times faster.

#B. Effects that differ across the distribution

We fit seven quantile levels, $\tau \in \{0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95\}$, with the ten covariates above, plus ordinary least squares (OLS) for the mean. For uncertainty we use the classical tool, the pairs bootstrap: resample births with replacement, refit, repeat 200 times. Plotting each coefficient against $\tau$ gives the quantile process: a flat line means the covariate shifts the whole distribution; a slope means it stretches or skews it.

TAUS = np.array([0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95])
K = len(TAUS)
beta_hat = np.array([rq_fn(X, y, t) for t in TAUS])               # K x p
ols = np.linalg.lstsq(X, y, rcond=None)[0]
ols_se = np.sqrt(np.diag(np.linalg.inv(X.T @ X)) * np.sum((y - X @ ols) ** 2) / (N_TRAIN - X.shape[1]))
N_BOOT = 200
rng_b = np.random.default_rng(RANDOM_SEED + 1)
t0 = time.time()
boot = np.empty((N_BOOT, K, X.shape[1]))
for bb in range(N_BOOT):
    i = rng_b.integers(0, N_TRAIN, N_TRAIN)
    boot[bb] = [rq_fn(X[i], y[i], t) for t in TAUS]
print(f"{N_BOOT * K:,} bootstrap fits in {time.time() - t0:.0f} s")
tab = pd.DataFrame(beta_hat.T * 1000, index=COEF, columns=[f"q{t:g}" for t in TAUS]).round(0)
tab["OLS"] = (ols * 1000).round(0)
print("coefficients in grams (intercept: a 30-year-old, BMI 25, 30 lb gain, girl, all other 0)")
print(tab.to_string())
print("bootstrap standard errors (g)")
print(pd.DataFrame(boot.std(0).T * 1000, index=COEF, columns=[f"q{t:g}" for t in TAUS]).round(0).to_string())

1,400 bootstrap fits in 19 s
coefficients in grams (intercept: a 30-year-old, BMI 25, 30 lb gain, girl, all other 0)
            q0.05    q0.1   q0.25    q0.5   q0.75    q0.9   q0.95     OLS
intercept  2390.0  2670.0  2994.0  3297.0  3583.0  3837.0  4010.0  3267.0
boy         101.0    73.0    68.0    94.0   114.0   127.0   119.0    93.0
black      -267.0  -236.0  -175.0  -170.0  -170.0  -147.0  -155.0  -188.0
age        -143.0  -102.0   -53.0   -13.0    29.0    34.0    28.0   -27.0
first      -120.0  -147.0  -111.0  -108.0   -90.0   -70.0   -93.0  -110.0
college     157.0   164.0   108.0    79.0    59.0    28.0    35.0    84.0
smoker     -110.0  -200.0  -157.0  -232.0  -298.0  -237.0  -212.0  -214.0
bmi5         17.0    25.0    43.0    55.0    73.0    82.0   106.0    54.0
gain10      107.0    95.0    78.0    81.0    84.0    81.0    96.0    89.0
pre_diab     15.0  -155.0    -6.0   -31.0    29.0   114.0    75.0    10.0
gest_hyp   -420.0  -291.0  -235.0  -224.0  -230.0  -203.0  -212.0  -235.0
bootstrap standard errors (g)
           q0.05   q0.1  q0.25  q0.5  q0.75  q0.9  q0.95
intercept   36.0   33.0   19.0  13.0   15.0  23.0   24.0
boy         42.0   33.0   18.0  16.0   17.0  21.0   26.0
black       74.0   49.0   27.0  18.0   23.0  29.0   34.0
age         33.0   31.0   19.0  17.0   17.0  20.0   29.0
first       44.0   34.0   18.0  17.0   17.0  23.0   27.0
college     52.0   32.0   22.0  17.0   17.0  23.0   29.0
smoker      83.0   88.0   55.0  36.0   47.0  60.0  108.0
bmi5        18.0   11.0    8.0   5.0    6.0   8.0   12.0
gain10      15.0   10.0    6.0   5.0    6.0   8.0   12.0
pre_diab   171.0  139.0   88.0  65.0  108.0  82.0  112.0
gest_hyp    83.0   52.0   37.0  25.0   32.0  30.0   51.0

SHOW = ["age", "pre_diab", "bmi5", "smoker", "gest_hyp", "college"]
LABEL = {"age": "mother 10 years older", "pre_diab": "diabetes before pregnancy", "bmi5": "BMI + 5",
         "smoker": "smoked in pregnancy", "gest_hyp": "gestational hypertension",
         "college": "bachelor's degree or more", "boy": "boy", "black": "Black mother",
         "first": "first birth", "gain10": "10 lb more weight gain", "intercept": "intercept"}
fig, axes = plt.subplots(2, 3, figsize=(14, 7.5), sharex=True)
for ax, cv in zip(axes.ravel(), SHOW):
    j = COEF.index(cv)
    lo, hi = np.quantile(boot[:, :, j], [0.05, 0.95], axis=0) * 1000
    ax.fill_between(TAUS, lo, hi, color=BLUE, alpha=0.25, lw=0, label="bootstrap 90%")
    ax.plot(TAUS, beta_hat[:, j] * 1000, "o-", color=BLUE, lw=2, label="quantile regression")
    ax.axhspan((ols[j] - 1.645 * ols_se[j]) * 1000, (ols[j] + 1.645 * ols_se[j]) * 1000, color=MUTED, alpha=0.25, lw=0)
    ax.axhline(ols[j] * 1000, color=MUTED, lw=1.5, ls="--", label="OLS (mean) with 90%")
    ax.axhline(0, color=INK, lw=0.6)
    ax.set_title(LABEL[cv], fontsize=10)
axes[0, 0].legend(fontsize=8)
for ax in axes[1]:
    ax.set_xlabel("quantile level tau")
for ax in axes[:, 0]:
    ax.set_ylabel("effect on the tau-quantile (g)")
fig.suptitle("The quantile process: classical quantile regression on 6,000 US births of 2023", fontsize=11);

Reading the panels (blue bands: bootstrap 90% intervals; grey: the OLS estimate with its 90% interval):

• Mother's age is the textbook case of an effect that a mean regression hides. At the median ten more years change almost nothing (-13 g, interval includes zero). In the lower tail they cost 100-140 g (5% and 10% quantiles); in the upper tail they add about 30 g. OLS reports a single -27 g, an average of opposite movements: babies of older mothers are more spread out - more small ones, a few more large ones - not uniformly lighter.
• BMI stretches the distribution upwards: +17 g per 5 BMI units at the 5% quantile, +55 g at the median and about +105 g at the 95% quantile. This is the "shifts the upper tail more than the median" pattern; the mean (+54 g) misses it completely.
• A degree matters most at the bottom (+160 g at the 5-10% quantiles, about +30 g at 90-95%), and gestational hypertension is a nearly constant -200 to -235 g except in the lower tail (-290 g at 10%, -420 g at 5%).
• Smoking lowers every quantile, by roughly 100-300 g; with 186 smokers in the training data the bands are too wide to say whether the effect varies with $\tau$.
• Diabetes before pregnancy is where one might expect a widening (larger babies, but also more complications): the estimates swing from -155 g at the 10% quantile to +114 g at the 90% - but with 76 diabetic mothers the bands are 200-500 g wide, and only the one at 90% excludes zero (just: it starts at about +10 g). Rare binary covariates in extreme quantiles are data-hungry.

#C. The asymmetric Laplace as a working likelihood

To make quantile regression Bayesian we need a likelihood. The asymmetric Laplace distribution (ALD) with quantile level $\tau$, location $\mu$ and scale $\sigma$ has density

$$p(y \mid \mu, \sigma, \tau) = \frac{\tau(1-\tau)}{\sigma} \exp\Bigl\{-\frac{\rho_\tau(y - \mu)}{\sigma}\Bigr\},$$

so with $\mu_i = x_i^\top\beta$ the log-likelihood is $-\frac1\sigma\sum_i\rho_\tau(y_i - x_i^\top\beta)$ plus terms without $\beta$: maximising it over $\beta$ is quantile regression, and $\mu$ is the $\tau$-quantile of the ALD (Yu & Moyeed 2001). PyMC has it as pm.AsymmetricLaplace(mu, b, q=tau); careful, its b is a rate in a symmetric parameterisation: the check-loss form above corresponds to $b = \sqrt{\tau(1-\tau)}/\sigma$. We verify that, and that $P(y < \mu) = \tau$, before using it.

q_, sig_, mu_ = 0.9, 0.3, 0.1
yv = np.array([-1.0, 0.0, 0.5, 2.0])
lp_pymc = pm.logp(pm.AsymmetricLaplace.dist(mu=mu_, b=np.sqrt(q_ * (1 - q_)) / sig_, q=q_), yv).eval()
lp_hand = np.log(q_ * (1 - q_) / sig_) - check_loss(yv - mu_, q_) / sig_
draws_ = pm.draw(pm.AsymmetricLaplace.dist(mu=mu_, b=np.sqrt(q_ * (1 - q_)) / sig_, q=q_), 200_000,
                 random_seed=RANDOM_SEED)
print("PyMC logp:", lp_pymc.round(6), "\ncheck loss:", lp_hand.round(6), f"\nP(y < mu) by simulation: {np.mean(draws_ < mu_):.4f}")

PyMC logp: [-1.570639 -1.237306 -2.403973 -6.903973] 
check loss: [-1.570639 -1.237306 -2.403973 -6.903973] 
P(y < mu) by simulation: 0.9000

One model holds all seven quantile levels side by side (they share nothing, so this is seven separate fits run together). Priors: $\beta \sim N(0, 5)$ in kilograms - effectively flat - and a half-normal on each $\sigma_\tau$.

Y_rep = np.repeat(y[:, None], K, axis=1)
t0 = time.time()
with pm.Model(coords={"coef": COEF, "tau": TAUS}) as m_ald:
    beta = pm.Normal("beta", 0.0, 5.0, dims=("coef", "tau"))
    sigma = pm.HalfNormal("sigma", 1.0, dims="tau")
    pm.AsymmetricLaplace("y", mu=pt.as_tensor(X) @ beta, b=np.sqrt(TAUS * (1 - TAUS)) / sigma, q=TAUS,
                         observed=Y_rep)
    idata_ald = pm.sample(random_seed=RANDOM_SEED, progressbar=False)
t_ald = time.time() - t0
s_ald = az.summary(idata_ald, var_names=["beta"], round_to=4)
print(f"ALD, 7 quantile levels: {t_ald:.0f} s; divergences {int(idata_ald.sample_stats['diverging'].sum())}; "
      f"max r_hat {s_ald['r_hat'].max():.3f}; min bulk ESS {s_ald['ess_bulk'].min():.0f}; "
      f"tuning steps {idata_ald.posterior.attrs['tuning_steps']}")
post_ald = idata_ald.posterior["beta"].to_numpy().reshape(-1, len(COEF), K)        # S x p x K
sig_ald = idata_ald.posterior["sigma"].to_numpy().reshape(-1, K)

ALD, 7 quantile levels: 23 s; divergences 0; max r_hat 1.005; min bulk ESS 2205; tuning steps 400

The posterior mode is the classical estimate. With a flat prior and $\sigma$ profiled out, the log posterior of $\beta$ is $-n\log\sum_i\rho_\tau(y_i - x_i^\top\beta)$, a decreasing function of the check loss, so the linear-programming solution must beat every posterior draw:

rows = []
for k, tau in enumerate(TAUS):
    L_draws = check_loss(y[:, None] - X @ post_ald[:, :, k].T, tau).sum(0)
    L_hat = check_loss(y - X @ beta_hat[k], tau).sum()
    sd = post_ald[:, :, k].std(0)
    rows.append({"tau": tau, "check loss at LP (kg)": L_hat, "smallest over 4000 draws": L_draws.min(),
                 "draws beating LP": int(np.sum(L_draws < L_hat - 1e-9)),
                 "max |post. mean - LP| / post. sd": np.max(np.abs(post_ald[:, :, k].mean(0) - beta_hat[k]) / sd),
                 "coef.": COEF[np.argmax(np.abs(post_ald[:, :, k].mean(0) - beta_hat[k]) / sd)]})
print(pd.DataFrame(rows).round(3).to_string(index=False))

 tau  check loss at LP (kg)  smallest over 4000 draws  draws beating LP  max |post. mean - LP| / post. sd    coef.
0.05                400.772                   400.841                 0                             0.458    black
0.10                622.485                   622.563                 0                             0.699    first
0.25               1000.613                  1000.766                 0                             0.296   smoker
0.50               1174.427                  1174.518                 0                             0.206  college
0.75                912.136                   912.244                 0                             0.287    first
0.90                500.562                   500.641                 0                             0.516 gest_hyp
0.95                295.952                   295.987                 0                             0.335      age

At every level the linear-programming solution has a smaller check loss than all 4,000 draws: it is the posterior mode (the $N(0, 5)$ prior is far too weak to move it). The posterior mean is a different point estimate, up to about 0.7 posterior sd away (the "first birth" coefficient at $\tau = 0.1$), because the posterior of a quantile-regression coefficient is not symmetric when a handful of births near the fitted line decide it. So far, so good: the working likelihood puts its centre in the right place.

But the spread is another matter. Compare the posterior sd of each coefficient with its bootstrap standard error, which is an honest (if approximate) frequentist measure of how much the estimate would move in a new sample:

sd_ald = post_ald.std(0).T                                   # K x p
sd_boot = boot.std(0)
ratio_ald = pd.DataFrame(sd_ald / sd_boot, index=[f"tau {t:g}" for t in TAUS], columns=COEF)
print("posterior sd / bootstrap sd, ALD with sigma learned")
print(ratio_ald.round(2).to_string())

posterior sd / bootstrap sd, ALD with sigma learned
          intercept   boy  black   age  first  college  smoker  bmi5  gain10  pre_diab  gest_hyp
tau 0.05       0.46  0.43   0.35  0.45   0.45     0.45    0.42  0.40    0.40      0.37      0.36
tau 0.1        0.52  0.51   0.50  0.52   0.56     0.60    0.56  0.54    0.52      0.49      0.53
tau 0.25       0.66  0.75   0.69  0.67   0.80     0.69    0.71  0.67    0.72      0.70      0.66
tau 0.5        0.87  0.76   0.90  0.77   0.84     0.81    0.92  0.84    0.76      0.79      0.77
tau 0.75       0.75  0.71   0.76  0.72   0.77     0.78    0.73  0.76    0.77      0.63      0.69
tau 0.9        0.54  0.55   0.58  0.59   0.59     0.59    0.59  0.59    0.57      0.54      0.58
tau 0.95       0.50  0.48   0.49  0.44   0.51     0.47    0.35  0.48    0.41      0.50      0.44

Every posterior sd is too small: 0.35-0.46 of the bootstrap sd at $\tau = 0.05$, 0.76-0.92 at the median, 0.35-0.51 at $\tau = 0.95$. A nominal 90% credible interval in the tails is less than half as wide as it should be. Within one quantile level the ratio is about the same for every coefficient - a hint that it is a property of the residual distribution at that quantile rather than of the covariate, which Part D makes precise.

#The units trap

Many implementations fix $\sigma = 1$ ("it cancels in the argmax"). It does cancel in the argmax - and nowhere else. With $\sigma$ fixed the posterior is $\propto \exp\{-\sum_i\rho_\tau(r_i)\}$, whose curvature is proportional to the units of $y$. Fit the median regression twice, once in kilograms and once in grams, both with $\sigma = 1$:

unit_rows = {}
for unit, fac in (("kg", 1.0), ("g", 1000.0)):
    with pm.Model(coords={"coef": COEF}):
        b_ = pm.Normal("beta", 0.0, 5.0 * fac, dims="coef")
        pm.AsymmetricLaplace("y", mu=pt.as_tensor(X) @ b_, b=np.sqrt(0.25) / 1.0, q=0.5, observed=y * fac)
        idu = pm.sample(random_seed=RANDOM_SEED, progressbar=False)
    unit_rows[f"sigma = 1 {unit}"] = idu.posterior["beta"].to_numpy().reshape(-1, len(COEF)).std(0) * 1000 / fac
k50 = list(TAUS).index(0.5)
unit_rows["sigma learned"] = sd_ald[k50] * 1000
unit_rows["bootstrap"] = sd_boot[k50] * 1000
print("posterior sd of median-regression coefficients, in grams")
print(pd.DataFrame(unit_rows, index=COEF).loc[["intercept", "age", "smoker", "bmi5", "pre_diab"]].round(2).to_string())

posterior sd of median-regression coefficients, in grams
           sigma = 1 kg  sigma = 1 g  sigma learned  bootstrap
intercept         25.30         0.76          11.34      13.11
age               27.58         1.19          12.91      16.83
smoker            73.03         3.25          33.43      36.39
bmi5              10.59         0.34           4.48       5.32
pre_diab         125.62         4.11          51.21      64.75

Same data, same point estimate, posterior sds that differ by a factor of roughly 20-35 between the two units (the square-root rule, $\sqrt{1000} \approx 32$, assumes a smooth, near-Gaussian posterior; one about a gram wide need not be). With $\sigma$ = 1 kg the sds are 1.6-2 times the bootstrap's; with $\sigma$ = 1 g they are 11-17 times too small. Learning $\sigma$ removes the arbitrariness, but, as the previous table showed, not the error.

#D. How wrong? A simulation with a known truth

On real data we can only compare with the bootstrap, itself an approximation. To measure coverage we need a truth. We simulate 200 datasets of $n = 200$ from a heteroscedastic, skewed model,

$$y = 1 + 2x + (0.5 + 1.5x)\,\varepsilon,\qquad x \sim U(0, 1),\qquad \varepsilon \sim \frac{\text{Gamma}(3) - 3}{\sqrt3},$$

whose $\tau$-quantile is exactly linear: $Q_\tau(y \mid x) = (1 + 0.5\,q_\tau) + (2 + 1.5\,q_\tau)\,x$ with $q_\tau$ the $\tau$-quantile of $\varepsilon$. For each dataset and $\tau \in \{0.5, 0.9\}$ we compute four intervals for the intercept and slope:

1. ALD posterior with $\sigma$ learned (flat prior on $\beta$, Jeffreys prior $1/\sigma$; $\sigma$ then integrates out analytically and the posterior of $\beta$ is $\propto L(\beta)^{-n}$, $L$ the total check loss);
2. the same draws adjusted by the sandwich formula of Yang, Wang & He (2016, International Statistical Review) - below;
3. a Gibbs posterior $\propto \exp\{-w\,L(\beta)\}$ with learning rate $w$ calibrated by the bootstrap - below;
4. the bootstrap percentile interval (200 resamples) as a frequentist reference.

Posteriors here are sampled with a vectorised random-walk Metropolis sampler (4 chains per dataset, all 200 datasets at once, with a pilot run to set the proposal) because 1,600 NUTS runs would take too long; the model has two parameters, so this is safe, and we check it against NUTS on one dataset.

Why the ALD width is wrong - the sandwich. For large $n$ the classical estimator is normal with covariance $\frac{\tau(1-\tau)}{n} D_1^{-1} D_0 D_1^{-1}$, where $D_0 = E[xx^\top]$ and $D_1 = E[f(Q_\tau \mid x)\,xx^\top]$ involves the conditional density at the quantile. The ALD posterior instead has covariance $\frac{\sigma}{n} D_1^{-1}$: it gets the "bread" $D_1$ right and the "meat" wrong. The two agree only if $\sigma f(Q_\tau\mid x) = \tau(1-\tau)$ everywhere - i.e. if the data really are ALD. Yang, Wang & He's fix follows by solving for the true covariance in terms of the posterior one, $\Sigma_{\rm adj} = \frac{n\,\tau(1-\tau)}{\sigma^2}\,\Sigma_{\rm post}\,\hat D_0\,\Sigma_{\rm post}$ - it needs nothing but the draws, $\hat\sigma$ and $X^\top X / n$. For our simulation model we can evaluate both formulas exactly, which predicts what the simulation should find.

def eps_q(t):
    return (stats.gamma(3).ppf(t) - 3) / np.sqrt(3)


def eps_pdf(e):
    return stats.gamma(3).pdf(np.sqrt(3) * e + 3) * np.sqrt(3)


SIM_TAUS = (0.5, 0.9)
theory = {}
for tau in SIM_TAUS:
    qe = eps_q(tau)
    rho_mean = integrate.quad(lambda e: check_loss(e - qe, tau) * eps_pdf(e), -np.sqrt(3), 40)[0]
    sig_star = rho_mean * integrate.quad(lambda x: 0.5 + 1.5 * x, 0, 1)[0]         # pseudo-true ALD scale
    D0 = np.array([[1, 0.5], [0.5, 1 / 3]])
    D1 = eps_pdf(qe) * np.array([[integrate.quad(lambda x: x ** (a + b) / (0.5 + 1.5 * x), 0, 1)[0]
                                  for b in range(2)] for a in range(2)])
    D1i = np.linalg.inv(D1)
    true_cov, ald_cov = tau * (1 - tau) * D1i @ D0 @ D1i, sig_star * D1i
    theory[tau] = np.sqrt(np.diag(ald_cov) / np.diag(true_cov))
    print(f"tau = {tau}: predicted (ALD posterior sd) / (true sd): intercept {theory[tau][0]:.2f}, slope {theory[tau][1]:.2f}")

tau = 0.5: predicted (ALD posterior sd) / (true sd): intercept 1.00, slope 0.80
tau = 0.9: predicted (ALD posterior sd) / (true sd): intercept 0.66, slope 0.52

def rwm(Xs, Ys, tau, w=None, init=None, prop_cov=None, rng=None, chains=4, n_steps=2500, burn=500):
    """Random-walk Metropolis for all datasets at once. w=None: ALD posterior with sigma integrated out
    (log p = -n log L); otherwise the Gibbs posterior log p = -w L. Returns draws (R, chains*kept, p)."""
    R_, n_, p_ = Xs.shape

    def logp(b):
        L = check_loss(Ys[:, None, :] - np.einsum("rnp,rcp->rcn", Xs, b), tau).sum(-1)
        return -n_ * np.log(L) if w is None else -w[:, None] * L

    def run(b, chol, steps, keep):
        lp = logp(b)
        out, acc = [], 0.0
        for t in range(steps):
            prop = b + np.einsum("rij,rcj->rci", chol, rng.normal(size=b.shape))
            lpp = logp(prop)
            a = np.log(rng.random(lp.shape)) < lpp - lp
            b, lp = np.where(a[..., None], prop, b), np.where(a, lpp, lp)
            if t >= steps - keep:
                out.append(b)
                acc += a.mean()
        return b, np.stack(out, 2), acc / keep

    chol = np.linalg.cholesky(prop_cov * 2.38 ** 2 / p_ * 0.25)
    b0 = init[:, None, :] + np.einsum("rij,rcj->rci", chol, rng.normal(size=(R_, chains, p_)))
    b0, pilot, _ = run(b0, chol, 600, 400)                                  # pilot: learn the scale
    pf = pilot.reshape(R_, -1, p_)
    pc = np.einsum("rsi,rsj->rij", pf - pf.mean(1, keepdims=True), pf - pf.mean(1, keepdims=True)) / pf.shape[1]
    chol = np.linalg.cholesky(pc * 2.38 ** 2 / p_ + 1e-12 * np.eye(p_))
    _, draws, acc = run(b0, chol, n_steps, n_steps - burn)
    return draws.reshape(R_, -1, p_), acc


R_SIM, N_SIM, B_SIM = 200, 200, 100
rng_s = np.random.default_rng(RANDOM_SEED + 2)
xs = rng_s.uniform(0, 1, size=(R_SIM, N_SIM))
Xs = np.stack([np.ones_like(xs), xs], -1)
Ys = 1 + 2 * xs + (0.5 + 1.5 * xs) * (rng_s.gamma(3, size=(R_SIM, N_SIM)) - 3) / np.sqrt(3)
LEVELS = np.array([0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99])
cover, widths, sim_draws = {}, {}, {}
t0 = time.time()
for tau in SIM_TAUS:
    truth = np.array([1 + 0.5 * eps_q(tau), 2 + 1.5 * eps_q(tau)])
    bhat = np.array([rq_fn(Xs[r], Ys[r], tau) for r in range(R_SIM)])
    bs = np.empty((R_SIM, B_SIM, 2))
    for r in range(R_SIM):
        for bb in range(B_SIM):
            i = rng_s.integers(0, N_SIM, N_SIM)
            bs[r, bb] = rq_fn(Xs[r, i], Ys[r, i], tau)
    bcov = np.einsum("rbi,rbj->rij", bs - bs.mean(1, keepdims=True), bs - bs.mean(1, keepdims=True)) / (B_SIM - 1)
    dr_ald, acc_ald = rwm(Xs, Ys, tau, init=bhat, prop_cov=bcov, rng=rng_s)
    Sig = np.einsum("rsi,rsj->rij", dr_ald - dr_ald.mean(1, keepdims=True), dr_ald - dr_ald.mean(1, keepdims=True)) / dr_ald.shape[1]
    sig_hat = check_loss(Ys - np.einsum("rnp,rp->rn", Xs, bhat), tau).sum(-1) / N_SIM
    D0h = np.einsum("rni,rnj->rij", Xs, Xs) / N_SIM
    Sadj = N_SIM * tau * (1 - tau) / sig_hat[:, None, None] ** 2 * Sig @ D0h @ Sig
    # learning rate: w = (1 / sigma) x geometric mean over coefficients of (posterior var / bootstrap var)
    vr = np.diagonal(Sig, axis1=1, axis2=2) / np.diagonal(bcov, axis1=1, axis2=2)
    w = np.exp(np.log(vr).mean(1)) / sig_hat
    dr_gibbs, acc_g = rwm(Xs, Ys, tau, w=w, init=bhat, prop_cov=bcov, rng=rng_s)
    sim_draws[tau] = (bhat, dr_ald, Sadj, dr_gibbs, bs)
    for name in ("ALD posterior", "ALD + sandwich (YWH)", "Gibbs posterior, bootstrap w", "bootstrap percentile"):
        cv, wd = [], []
        for lev in LEVELS:
            a = (1 - lev) / 2
            if name == "ALD posterior":
                lo, hi = np.quantile(dr_ald, [a, 1 - a], axis=1)
            elif name == "ALD + sandwich (YWH)":
                half = stats.norm.ppf(1 - a) * np.sqrt(np.diagonal(Sadj, axis1=1, axis2=2))
                lo, hi = dr_ald.mean(1) - half, dr_ald.mean(1) + half
            elif name == "Gibbs posterior, bootstrap w":
                lo, hi = np.quantile(dr_gibbs, [a, 1 - a], axis=1)
            else:
                lo, hi = np.quantile(bs, [a, 1 - a], axis=1)
            cv.append(((lo <= truth) & (truth <= hi)).mean(0))
            wd.append(np.median(hi - lo, 0))
        cover[(tau, name)], widths[(tau, name)] = np.array(cv), np.array(wd)
    emp_sd = bhat.std(0)
    print(f"tau = {tau}: acceptance ALD {acc_ald:.2f}, Gibbs {acc_g:.2f}; empirical sd of the estimates "
          f"(intercept, slope) {emp_sd.round(3)}; median ALD posterior sd {np.median(np.sqrt(np.diagonal(Sig, axis1=1, axis2=2)), 0).round(3)} "
          f"-> ratio {(np.median(np.sqrt(np.diagonal(Sig, axis1=1, axis2=2)), 0) / emp_sd).round(2)} (theory {theory[tau].round(2)})")
print(f"simulation: {time.time() - t0:.0f} s")
rows = []
for (tau, name), cv in cover.items():
    i90 = list(LEVELS).index(0.9)
    rows.append({"tau": tau, "method": name, "90% coverage: intercept": cv[i90, 0], "slope": cv[i90, 1],
                 "median width: intercept": widths[(tau, name)][i90, 0], "slope ": widths[(tau, name)][i90, 1]})
print(pd.DataFrame(rows).round(3).to_string(index=False))

tau = 0.5: acceptance ALD 0.35, Gibbs 0.35; empirical sd of the estimates (intercept, slope) [0.14  0.364]; median ALD posterior sd [0.133 0.281] -> ratio [0.95 0.77] (theory [1.  0.8])

tau = 0.9: acceptance ALD 0.35, Gibbs 0.35; empirical sd of the estimates (intercept, slope) [0.286 0.752]; median ALD posterior sd [0.186 0.383] -> ratio [0.65 0.51] (theory [0.66 0.52])
simulation: 42 s
 tau                       method  90% coverage: intercept  slope  median width: intercept  slope 
 0.5                ALD posterior                    0.905  0.795                    0.440   0.920
 0.5         ALD + sandwich (YWH)                    0.900  0.875                    0.443   1.180
 0.5 Gibbs posterior, bootstrap w                    0.920  0.865                    0.506   1.074
 0.5         bootstrap percentile                    0.885  0.875                    0.447   1.187
 0.9                ALD posterior                    0.775  0.645                    0.606   1.247
 0.9         ALD + sandwich (YWH)                    0.895  0.890                    0.923   2.407
 0.9 Gibbs posterior, bootstrap w                    0.945  0.875                    1.132   2.324
 0.9         bootstrap percentile                    0.905  0.900                    0.903   2.519

A check of the Metropolis sampler against NUTS on the first simulated dataset at $\tau = 0.9$ (PyMC, pm.AsymmetricLaplace with a flat-ish prior on $\beta$ and a Jeffreys-like prior on $\log\sigma$):

with pm.Model():
    b_ = pm.Normal("beta", 0.0, 100.0, shape=2)
    log_s = pm.Flat("log_sigma")
    pm.AsymmetricLaplace("y", mu=pt.as_tensor(Xs[0]) @ b_, b=np.sqrt(0.09) / pt.exp(log_s), q=0.9, observed=Ys[0])
    id_chk = pm.sample(random_seed=RANDOM_SEED, progressbar=False)
nuts_d = id_chk.posterior["beta"].to_numpy().reshape(-1, 2)
_, dr0, _, _, _ = sim_draws[0.9]
print("NUTS  mean", nuts_d.mean(0).round(3), "sd", nuts_d.std(0).round(3))
print("RWM   mean", dr0[0].mean(0).round(3), "sd", dr0[0].std(0).round(3))

NUTS  mean [1.69 3.77] sd [0.2   0.333]
RWM   mean [1.688 3.767] sd [0.198 0.327]

fig, axes = plt.subplots(2, 2, figsize=(12, 8.5), sharex=True, sharey=True)
styles = {"ALD posterior": (ORANGE, "o"), "ALD + sandwich (YWH)": (BLUE, "s"),
          "Gibbs posterior, bootstrap w": (AQUA, "D"), "bootstrap percentile": (MUTED, "^")}
band = 1.96 * np.sqrt(LEVELS * (1 - LEVELS) / R_SIM)
for row, tau in enumerate(SIM_TAUS):
    for col, par in enumerate(("intercept", "slope")):
        ax = axes[row, col]
        ax.fill_between(LEVELS, LEVELS - band, LEVELS + band, color=LIGHT, lw=0, label="Monte Carlo 95% range")
        ax.plot([0.45, 1], [0.45, 1], color=INK, lw=0.8)
        for name, (c, mk) in styles.items():
            ax.plot(LEVELS, cover[(tau, name)][:, col], marker=mk, color=c, lw=1.5, ms=5, label=name)
        ax.set_title(f"tau = {tau}: {par}", fontsize=10)
        ax.set(xlim=(0.47, 1.0), ylim=(0.2, 1.01))
fig.supxlabel("nominal level of the interval", fontsize=10)
fig.supylabel("share of the 200 datasets whose interval covers the truth", fontsize=10)
axes[0, 0].legend(fontsize=8, loc="upper left")
fig.suptitle("Coverage of credible intervals for quantile-regression coefficients (n = 200, known truth)", fontsize=11);

What the simulation shows:

• The ALD posterior (orange) is too narrow, and the sandwich theory predicted by how much: its sd was 0.95 and 0.77 of the true sampling sd (intercept, slope) at $\tau = 0.5$ and 0.65 and 0.51 at $\tau = 0.9$, against 1.00, 0.80, 0.66 and 0.52 from the formulas. Nominal 90% intervals covered 90% and 80% at the median and only 78% and 65% at $\tau = 0.9$. The damage grows towards the tails, and is larger for the slope, whose information comes from the high-$x$ end where the noise is largest. The intercept at the median is right only by coincidence - there $\sigma f(Q_\tau\mid x)$ happens to average out to about $\tau(1-\tau)$.
• Sample size does not help: the ratio in the formulas does not depend on $n$. More data make the ALD intervals narrower and keep them equally overconfident.
• The sandwich (blue) and the bootstrap (grey) lie on the diagonal within Monte Carlo error at every level (89-90% at nominal 90%).
• The Gibbs posterior (green) with one bootstrap-calibrated learning rate per dataset fixes the slope (87.5% at 90%) but at $\tau = 0.9$ over-covers the intercept (94.5% at 90%, 71% at 50%): a single number rescales the whole posterior, and here the intercept and slope needed different factors (0.66 against 0.52). A learning rate is a scalar fix for a matrix-shaped problem; calibrate it on the quantity you will report, or use the sandwich.
• NUTS and the random-walk sampler agree on the first dataset (means within 0.003, sds within 0.006).

#E. Two fixes on the real data

Fix 1, the sandwich (Yang, Wang & He 2016). Keep the ALD draws, replace their covariance by $\Sigma_{\rm adj}$. Cheap and automatic; the adjusted interval is a normal approximation around the posterior mean.

Fix 2, a generalised (Gibbs) posterior. Drop the pretence of a likelihood. Bissiri, Holmes & Walker (2016, JRSS B) showed that updating beliefs with a loss is coherent: the posterior is $\pi(\beta \mid y) \propto \pi(\beta)\exp\{-w\,\sum_i\rho_\tau(y_i - x_i^\top\beta)\}$. The ALD with $\sigma$ fixed is the special case $w = 1/\sigma$; the learning rate $w$ is not given by the model and must be chosen. Syring & Martin (2019, Biometrika) choose it so that credible sets have their nominal frequentist coverage, estimated by the bootstrap. We use a simpler moment-matching version: posterior variance scales like $1/w$, so we set $w_\tau = \frac1{\hat\sigma_\tau}\times$ (geometric mean over coefficients of posterior variance / bootstrap variance) and refit - it is exactly what we did in the simulation. One number per $\tau$ can rescale the posterior but not reshape it, which is why the simulation's intercept and slope did not come out equally well.

sig_hat_k = sig_ald.mean(0)
D0_hat = X.T @ X / N_TRAIN
sd_ywh = np.empty_like(sd_ald)
for k, tau in enumerate(TAUS):
    S = np.cov(post_ald[:, :, k].T)
    sd_ywh[k] = np.sqrt(np.diag(N_TRAIN * tau * (1 - tau) / sig_hat_k[k] ** 2 * S @ D0_hat @ S))
w_tau = np.exp(np.log(sd_ald ** 2 / sd_boot ** 2).mean(1)) / sig_hat_k
print(pd.DataFrame({"learning rate w": w_tau, "ALD's 1 / sigma": 1 / sig_hat_k}, index=[f"tau {t:g}" for t in TAUS]).round(2).T.to_string())
t0 = time.time()
with pm.Model(coords={"coef": COEF, "tau": TAUS}) as m_gibbs:
    beta = pm.Normal("beta", 0.0, 5.0, dims=("coef", "tau"))
    resid = y[:, None] - pt.as_tensor(X) @ beta
    pm.Potential("loss", -(w_tau * (resid * (TAUS - (resid < 0)))).sum())
    idata_gibbs = pm.sample(random_seed=RANDOM_SEED, progressbar=False)
s_g = az.summary(idata_gibbs, var_names=["beta"], round_to=4)
print(f"Gibbs posterior: {time.time() - t0:.0f} s; divergences {int(idata_gibbs.sample_stats['diverging'].sum())}; "
      f"max r_hat {s_g['r_hat'].max():.3f}; min bulk ESS {s_g['ess_bulk'].min():.0f}")
post_gibbs = idata_gibbs.posterior["beta"].to_numpy().reshape(-1, len(COEF), K)
sd_gibbs = post_gibbs.std(0).T
summ = pd.DataFrame({"ALD": np.median(sd_ald / sd_boot, 1), "ALD + sandwich": np.median(sd_ywh / sd_boot, 1),
                     "Gibbs, calibrated w": np.median(sd_gibbs / sd_boot, 1)}, index=[f"tau {t:g}" for t in TAUS])
print("median over the 11 coefficients of (posterior sd / bootstrap sd)")
print(summ.round(2).to_string())
for name, sdm in (("sandwich", sd_ywh), ("Gibbs", sd_gibbs)):
    rr = sdm / sd_boot
    kmin, jmin = np.unravel_index(np.argmin(rr), rr.shape)
    kmax, jmax = np.unravel_index(np.argmax(rr), rr.shape)
    print(f"{name}: lowest ratio {rr.min():.2f} ({COEF[jmin]}, tau {TAUS[kmin]:g}), highest {rr.max():.2f} "
          f"({COEF[jmax]}, tau {TAUS[kmax]:g}); share of the 77 ratios within 0.8-1.25: {np.mean((rr > 0.8) & (rr < 1.25)):.0%}")

                 tau 0.05  tau 0.1  tau 0.25  tau 0.5  tau 0.75  tau 0.9  tau 0.95
learning rate w      2.52     2.71      2.93     3.42      3.53     3.95      4.26
ALD's 1 / sigma     14.95     9.63      5.99     5.10      6.57    11.97     20.25

Gibbs posterior: 4 s; divergences 0; max r_hat 1.004; min bulk ESS 1448
median over the 11 coefficients of (posterior sd / bootstrap sd)
           ALD  ALD + sandwich  Gibbs, calibrated w
tau 0.05  0.42            0.91                 1.06
tau 0.1   0.52            0.99                 0.97
tau 0.25  0.69            0.94                 1.01
tau 0.5   0.81            0.98                 0.99
tau 0.75  0.75            1.00                 1.02
tau 0.9   0.58            0.99                 1.02
tau 0.95  0.48            1.07                 1.05
sandwich: lowest ratio 0.70 (smoker, tau 0.05), highest 1.23 (college, tau 0.05); share of the 77 ratios within 0.8-1.25: 96%
Gibbs: lowest ratio 0.83 (smoker, tau 0.95), highest 1.37 (smoker, tau 0.05); share of the 77 ratios within 0.8-1.25: 99%

SHOW_E = ["age", "pre_diab", "smoker", "bmi5"]
fig, axes = plt.subplots(1, 4, figsize=(15, 4.6))
z90 = stats.norm.ppf(0.95)
xpos = np.arange(K)
for ax, cv in zip(axes, SHOW_E):
    j = COEF.index(cv)
    specs = [("ALD posterior", ORANGE, -0.24, np.quantile(post_ald[:, j, :], [0.05, 0.95], axis=0)),
             ("ALD + sandwich", BLUE, -0.08, np.array([post_ald[:, j, :].mean(0) - z90 * sd_ywh[:, j],
                                                        post_ald[:, j, :].mean(0) + z90 * sd_ywh[:, j]])),
             ("Gibbs, calibrated w", AQUA, 0.08, np.quantile(post_gibbs[:, j, :], [0.05, 0.95], axis=0)),
             ("bootstrap", MUTED, 0.24, np.quantile(boot[:, :, j], [0.05, 0.95], axis=0))]
    for name, c, off, (lo, hi) in specs:
        ax.vlines(xpos + off, lo * 1000, hi * 1000, color=c, lw=3.5, label=name)
    ax.plot(xpos, beta_hat[:, j] * 1000, "o", color=INK, ms=3, label="LP estimate")
    ax.axhline(0, color=INK, lw=0.5)
    ax.set(title=LABEL[cv], xlabel="quantile level tau", xticks=xpos, xticklabels=[f"{t:g}" for t in TAUS])
axes[0].set_ylabel("effect on the quantile (g)")
axes[0].legend(fontsize=7, loc="lower right")
fig.suptitle("90% intervals for quantile effects: raw ALD posterior vs two fixes vs the bootstrap", fontsize=11);

rows = []
for cv, tau in (("smoker", 0.05), ("pre_diab", 0.9), ("age", 0.95), ("bmi5", 0.05), ("age", 0.5)):
    j, k = COEF.index(cv), list(TAUS).index(tau)
    m_ = post_ald[:, j, k].mean()
    ivs = {"ALD": np.quantile(post_ald[:, j, k], [0.05, 0.95]), "sandwich": m_ + np.array([-z90, z90]) * sd_ywh[k, j],
           "Gibbs": np.quantile(post_gibbs[:, j, k], [0.05, 0.95]), "bootstrap": np.quantile(boot[:, k, j], [0.05, 0.95])}
    rows.append({"effect": f"{cv} at tau {tau:g}", "LP (g)": beta_hat[k, j] * 1000,
                 **{nm: f"{lo * 1000:5.0f} to {hi * 1000:4.0f}" for nm, (lo, hi) in ivs.items()}})
print("90% intervals (g)")
print(pd.DataFrame(rows).round(0).to_string(index=False))

90% intervals (g)
             effect  LP (g)           ALD      sandwich         Gibbs     bootstrap
 smoker at tau 0.05  -110.0  -160 to  -44  -197 to   -6  -309 to   66  -210 to   33
pre_diab at tau 0.9   114.0    45 to  192    11 to  217   -18 to  272    11 to  241
    age at tau 0.95    28.0    12 to   53   -15 to   79   -12 to   79   -14 to   78
   bmi5 at tau 0.05    17.0     7 to   30   -11 to   48   -12 to   46   -10 to   46
     age at tau 0.5   -13.0   -31 to   11   -38 to   18   -36 to   16   -38 to   16

Both fixes do their job on the real data. The median ratio of posterior sd to bootstrap sd is 0.91-1.07 for the sandwich and 0.97-1.06 for the calibrated Gibbs posterior, against 0.42-0.81 for the raw ALD. The learning rates (2.5-4.3) are far below the ALD's implicit $1/\sigma_\tau$ (5-20): the ALD was learning from the data 1.5-6 times too fast, most of all in the tails. Individual coefficients scatter more (sandwich 0.70-1.23, Gibbs 0.83-1.37). The extremes are mostly the 186 smokers at the 5% and 95% levels (and the degree coefficient at 5% for the sandwich), where the bootstrap itself is noisy - neither side of the ratio is exact there.

Does it change conclusions? The table picks five effects. The raw ALD says that at the 5% quantile smoking clearly lowers birth weight (-160 to -44 g), that older mothers have heavier babies at the 95% quantile (+12 to +53 g), that BMI raises the 5% quantile (+7 to +30 g) and that pre-pregnancy diabetes raises the 90% quantile (+45 to +192 g). With honest intervals the age and BMI effects include zero under all three methods; the smoking interval includes zero under the bootstrap and the Gibbs posterior and only just excludes it under the sandwich (-6 g); the diabetes interval starts at +11 g (sandwich, bootstrap) or -18 g (Gibbs). At the median, where the ALD is least wrong, all four agree about age. The tails - where quantile regression is supposed to earn its keep - are exactly where the raw posterior overstates what the data say.

#F. Quantile crossing: birth weight for gestational age

Doctors judge a newborn's size for its gestational age: "small for gestational age" usually means below the 10th centile of birth weight among babies born at the same week. Those centile curves are a quantile regression of birth weight on gestational age. We fit them with a cubic B-spline in the obstetric estimate of gestation (24-42 weeks; knots at 30, 34, 37 and 39 weeks) on the training births, first one quantile level at a time, at 19 levels from 5% to 95%.

Nothing forces separately fitted quantile curves to stay in order. Where data are plentiful they do; where data are thin they can cross - a "10th centile" above the "15th" - which is logically impossible for a distribution.

G_LO, G_HI, KNOTS = 24.0, 42.0, [30, 34, 37, 39]


def bbasis(x, k=3):
    t = np.r_[[G_LO] * (k + 1), KNOTS, [G_HI] * (k + 1)]
    nb = len(t) - k - 1
    x = np.clip(np.asarray(x, float), G_LO, G_HI - 1e-9)
    return np.column_stack([BSpline(t, np.eye(nb)[j], k)(x) for j in range(nb)])


keep_g = (tr.oegest_comb >= G_LO) & (tr.oegest_comb <= G_HI)
g_tr, yg = tr.oegest_comb[keep_g].to_numpy(float), tr.bw[keep_g].to_numpy()
keep_gt = (te.oegest_comb >= G_LO) & (te.oegest_comb <= G_HI)
g_te, yg_te = te.oegest_comb[keep_gt].to_numpy(float), te.bw[keep_gt].to_numpy()
Bm, Bm_te = bbasis(g_tr), bbasis(g_te)
ggrid = np.linspace(G_LO, G_HI, 181)
Bg = bbasis(ggrid)
print(f"training births at 24-42 weeks: {len(g_tr):,} (dropped {int((~keep_g).sum())}); "
      f"below 28 weeks: {int(np.sum(g_tr < 28))}, below 32 weeks: {int(np.sum(g_tr < 32))}")
TG19 = np.round(np.arange(0.05, 0.951, 0.05), 2)
th19 = np.array([rq_fn(Bm, yg, t) for t in TG19])
Q19 = Bg @ th19.T                                                     # grid x 19
cross_grid = (np.diff(Q19, axis=1) < 0).any(1)
Q19_tr = Bm @ th19.T
print(f"grid points (0.1 week apart) where some pair of the 19 curves is out of order: {cross_grid.sum()} "
      f"(weeks {ggrid[cross_grid].min():.1f}-{ggrid[cross_grid].max():.1f})")
print(f"training births whose 19 fitted quantiles are not in order: {(np.diff(Q19_tr, axis=1) < 0).any(1).sum()}")
edges = np.flatnonzero(np.diff(np.r_[0, cross_grid.astype(int), 0]))
print("stretches with crossing (weeks):", [f"{ggrid[a]:.1f}-{ggrid[b - 1]:.1f}" for a, b in zip(edges[::2], edges[1::2])])
odd = (g_tr < 32) & (yg > 2.5)
print(f"training births before 32 weeks: {int(np.sum(g_tr < 32))}; at 42 weeks: {int(np.sum(g_tr == 42))}. "
      "Those before 32 weeks weighing over 2.5 kg:")
print(tr[keep_g][odd][["oegest_comb", "combgest", "dbwt"]].rename(
    columns={"oegest_comb": "obstetric estimate (wk)", "combgest": "combined estimate (wk)", "dbwt": "grams"}).to_string(index=False))
Q19_sorted = np.sort(Q19, axis=1)                                     # rearrangement