Bayesian Sequential Probit

Models

Two competing models for ordered categorical outcomes, estimated via data-augmentation Gibbs sampling and compared against MLE and bambi/PyMC (HMC). Applied to the Ashford (1959) pneumoconiosis data: 371 British coal miners classified into three ordered severity categories (normal, mild, severe) across 8 exposure-time groups.

Cumulative Ordered Probit — Albert & Chib (1993)

The cumulative model (Albert & Chib 1993) introduces a single latent utility $z_i^* \sim N(\beta_1 \tilde{x}_i, 1)$ and assigns $Y_i = j$ whenever $z_i^*$ falls between cut-points $\gamma_{j-1}$ and $\gamma_j$. A single slope $\beta_1$ governs all transitions simultaneously — the proportional-odds assumption. The posterior precision $\bar{B}^{-1} = X_c'X_c + B_0^{-1}$ is constant across iterations; cut-points are drawn via a uniform slice sampler that automatically preserves their ordering.

Sequential Ordinal Probit — Albert & Chib (2001)

The sequential model (Albert & Chib 2001) decomposes the ordinal outcome into $J-1$ successive binary decisions. At stage $k$, among subjects who reached that stage, each individual exits ($Y_i = k-1$) or continues ($Y_i \geq k$) via an independent binary probit with its own coefficient vector $\alpha_k$. Crucially, stages are conditionally independent given data — the Gibbs sampler reduces to $J-1$ separate Albert & Chib (1993) binary probit runs, each on a different at-risk subset.

Key Finding

The Stage-1 exposure slope $\alpha_{11} = 0.804$ is 2.6× larger than the Stage-2 slope $\alpha_{21} = 0.313$ (posterior probability $P(\alpha_{11} > \alpha_{21}) = 0.984$). The cumulative model's single $\beta_1 = 0.796$ is a weighted average of the two stages — it correctly estimates disease onset risk but overstates the effect on progression from mild to severe. The Stage-2 credible interval [−0.08, 0.72] includes zero: with only 82 at-risk miners, the data cannot conclusively establish a progression effect. All three inference methods (Gibbs, MLE, bambi HMC) agree on parameter estimates to within 0.016.

LPML (log pseudo-marginal likelihood) and LOO elpd (leave-one-out expected log predictive density) are both leave-one-out predictive scores on the same scale — higher (less negative) is better — the first from the Gibbs sampler, the second from the bambi/HMC fit via PSIS-LOO; the two agree closely (−210.24 vs −210.11). They are reported as elpd, not as the LOO-IC = −2·elpd information-criterion form, so they line up directly with the Gibbs LPML. AIC = −2·log-lik + 2·(parameters) is the in-sample counterpart. Differences of a nat or two are negligible.

Notebook

Downloads

Gibbs Sampler — Source Code

"""
seq_ord_probit_gibbs.py

Two Bayesian ordinal regression models via data-augmentation Gibbs sampling:

  1. Cumulative ordered probit  (Albert & Chib 1993)
       P(Y <= j | x) = Phi(gamma_j - x'beta)
       Single latent variable z* ~ N(x'beta, 1).
       Cut-points gamma_1 < ... < gamma_{J-1} drawn via uniform slice.

  2. Sequential ordinal probit  (Albert & Chib 2001)
       At stage k (k = 1,...,J-1), conditional on reaching stage k,
       a subject exits (Y = k-1) vs continues (Y >= k) via a binary probit:
           z_ik ~ N(x_i'alpha_k, 1),   s_ik = 1(z_ik > 0)
       J-1 independent binary probits — different alpha at every stage.

Joint probabilities (sequential, J = 3):
       P(Y=0|x) = Phi(-x'alpha_1)
       P(Y=1|x) = Phi(x'alpha_1) * Phi(-x'alpha_2)
       P(Y=2|x) = Phi(x'alpha_1) * Phi(x'alpha_2)

Design matrices
---------------
  Cumulative : X_c — NO intercept column (cut-points play the role of intercepts)
  Sequential : X_s — typically [1, covariates] (each stage has its own intercept)

References
----------
  Albert, J. H., & Chib, S. (1993).  JASA 88(422), 669-679.
  Albert, J. H., & Chib, S. (2001).  JASA 96(454), 415-425.
"""

import time
import numpy as np
from scipy.special import ndtr, ndtri


# ── truncated-Normal helper ────────────────────────────────────────────────────

def _rtn(lo, hi, mu, rng):
    """
    Vectorised draw from N(mu_i, 1) truncated to (lo_i, hi_i).
    Uses inverse-CDF; falls back to interval midpoint when mass < 1e-300.
    """
    p_lo = ndtr(lo - mu)
    p_hi = ndtr(hi - mu)
    gap  = p_hi - p_lo
    u    = rng.uniform(size=len(mu))
    p    = p_lo + u * gap
    tiny = gap < 1e-300
    p    = np.where(tiny, 0.5 * (p_lo + p_hi), p)
    return ndtri(np.clip(p, 1e-15, 1.0 - 1e-15)) + mu


# ── cumulative ordered probit ─────────────────────────────────────────────────

def cumulative_probit_gibbs(
    y, X_c,
    b0=None, B0=None,
    R=10000, burn=2000,
    seed=42, nprint=1000,
):
    """
    Albert & Chib (1993) Gibbs for the cumulative ordered probit.

    Model:  z_i ~ N(x_i'beta, 1)
            Y_i = j   iff   gamma[j-1] < z_i <= gamma[j]
            gamma[-1] = -inf,  gamma[J-1] = +inf  (sentinels)
            gamma[0] < ... < gamma[J-2]  — J-1 free cut-points

    X_c must NOT include an intercept column (cut-points absorb the intercept).

    Parameters
    ----------
    y   : (n,) int   labels in {0, ..., J-1}
    X_c : (n, k)     design matrix, no intercept
    b0  : (k,)       prior mean for beta       [default zeros]
    B0  : (k, k)     prior covariance          [default 100*I]
    R, burn : int    total / burn-in iterations
    seed, nprint : int

    Returns
    -------
    dict  'beta'  (R_kept, k) | 'gamma' (R_kept, J-1) | 'z' (R_kept, n)
    """
    rng = np.random.default_rng(seed)
    y   = np.asarray(y, int).ravel()
    X_c = np.asarray(X_c, float)
    n, k = X_c.shape
    J    = int(y.max()) + 1

    if b0 is None: b0 = np.zeros(k)
    if B0 is None: B0 = 100.0 * np.eye(k)
    b0        = np.asarray(b0, float)
    B0_inv    = np.linalg.inv(np.asarray(B0, float))
    B0_inv_b0 = B0_inv @ b0

    XtX  = X_c.T @ X_c
    Bbar = np.linalg.inv(XtX + B0_inv)
    L    = np.linalg.cholesky(Bbar)

    # initial cut-points: equally spaced between ±2
    gamma = np.linspace(-1.5, 1.5, J - 1)
    beta  = np.zeros(k)

    NINF = -1e10
    PINF =  1e10

    R_kept      = R - burn
    beta_draws  = np.empty((R_kept, k))
    gamma_draws = np.empty((R_kept, J - 1))
    z_draws     = np.empty((R_kept, n))
    kept = 0
    t0   = time.time()

    for r in range(1, R + 1):
        gext = np.r_[NINF, gamma, PINF]

        # Block 1: z_i | beta, gamma, y  — truncated Normal
        mu = X_c @ beta
        z  = _rtn(gext[y], gext[y + 1], mu, rng)

        # Block 2: beta | z  — Normal posterior (Bbar is constant)
        bbar = Bbar @ (X_c.T @ z + B0_inv_b0)
        beta = bbar + L @ rng.standard_normal(k)

        # Block 3: gamma_j | z, y  — uniform slice (independent per j)
        # Posterior: Uniform(max(z: y=j), min(z: y=j+1))
        # Ordering is automatically satisfied because z draws respect current gamma.
        for j in range(J - 1):
            lo_j = z[y == j].max()     if (y == j    ).any() else NINF
            hi_j = z[y == j + 1].min() if (y == j + 1).any() else PINF
            if lo_j < hi_j:
                gamma[j] = rng.uniform(lo_j, hi_j)

        if r > burn:
            beta_draws[kept]  = beta
            gamma_draws[kept] = gamma.copy()
            z_draws[kept]     = z
            kept += 1

        if nprint > 0 and r % nprint == 0:
            print(f'  [cumulative]  iter {r:6d}/{R}  ({time.time()-t0:.1f}s)')

    print(f'[cumulative] Done in {time.time()-t0:.1f}s  |  kept: {R_kept}')
    return {'beta': beta_draws, 'gamma': gamma_draws, 'z': z_draws}


# ── sequential ordinal probit ─────────────────────────────────────────────────

def sequential_probit_gibbs(
    y, X_s,
    b0=None, B0=None,
    R=10000, burn=2000,
    seed=42, nprint=1000,
):
    """
    Albert & Chib (2001) Gibbs for the sequential ordinal probit.

    Stage k (k = 1,...,J-1):  among subjects with Y_i >= k-1,
        define  s_ik = 1(Y_i >= k)   (binary: continue vs exit)
        z_ik ~ N(x_i'alpha_k, 1),   s_ik = 1(z_ik > 0)

    The J-1 stages share the same prior but have independent parameters.
    Each stage reduces to a binary probit (Albert & Chib 1993).

    Parameters
    ----------
    y   : (n,) int   labels in {0, ..., J-1}
    X_s : (n, k)     design matrix (typically includes intercept column)
    b0  : (k,)       prior mean for each alpha_k  [default zeros]
    B0  : (k, k)     prior covariance             [default 100*I]
    R, burn, seed, nprint : int

    Returns
    -------
    dict
      'alpha'    : list of J-1 arrays (R_kept, k)
      'z'        : list of J-1 arrays (R_kept, n_k)  — latent utilities per stage
      'at_risk'  : list of J-1 bool (n,) masks
    """
    rng = np.random.default_rng(seed)
    y   = np.asarray(y, int).ravel()
    X_s = np.asarray(X_s, float)
    n, k = X_s.shape
    J    = int(y.max()) + 1

    if b0 is None: b0 = np.zeros(k)
    if B0 is None: B0 = 100.0 * np.eye(k)
    b0        = np.asarray(b0, float)
    B0_inv    = np.linalg.inv(np.asarray(B0, float))
    B0_inv_b0 = B0_inv @ b0

    # build at-risk subsets for each stage
    at_risk = []   # bool masks length n
    s_list  = []   # binary outcome at each stage
    X_list  = []   # design submatrix
    for stage in range(1, J):              # stage = 1 .. J-1
        mask = (y >= stage - 1)
        at_risk.append(mask)
        s_list.append((y[mask] >= stage).astype(float))
        X_list.append(X_s[mask])

    # precompute posterior covariance for each stage (constant — sigma^2 = 1 fixed)
    Bbar_list = []
    L_list    = []
    for Xs in X_list:
        Bbar = np.linalg.inv(Xs.T @ Xs + B0_inv)
        Bbar_list.append(Bbar)
        L_list.append(np.linalg.cholesky(Bbar))

    alpha_list = [np.zeros(k) for _ in range(J - 1)]
    z_list     = [np.zeros(mask.sum()) for mask in at_risk]

    R_kept      = R - burn
    alpha_draws = [np.empty((R_kept, k)) for _ in range(J - 1)]
    z_draws_out = [np.empty((R_kept, mask.sum())) for mask in at_risk]
    kept = 0
    t0   = time.time()

    NINF = -1e10
    PINF =  1e10

    for r in range(1, R + 1):

        for q in range(J - 1):
            s    = s_list[q]
            Xs   = X_list[q]
            mu   = Xs @ alpha_list[q]
            lo   = np.where(s == 1, 0.0,  NINF)
            hi   = np.where(s == 1, PINF, 0.0)
            z    = _rtn(lo, hi, mu, rng)
            z_list[q] = z

            bbar = Bbar_list[q] @ (Xs.T @ z + B0_inv_b0)
            alpha_list[q] = bbar + L_list[q] @ rng.standard_normal(k)

        if r > burn:
            for q in range(J - 1):
                alpha_draws[q][kept]    = alpha_list[q]
                z_draws_out[q][kept]    = z_list[q]
            kept += 1

        if nprint > 0 and r % nprint == 0:
            print(f'  [sequential]  iter {r:6d}/{R}  ({time.time()-t0:.1f}s)')

    print(f'[sequential] Done in {time.time()-t0:.1f}s  |  kept: {R_kept}')
    return {'alpha': alpha_draws, 'z': z_draws_out, 'at_risk': at_risk}


# ── predicted probabilities ────────────────────────────────────────────────────

def cum_pred_probs(X_c, beta_draws, gamma_draws):
    """
    Posterior predictive probabilities for the cumulative model.

    Returns
    -------
    p_bar   : (n, J)    posterior mean P(Y=j | x_i)
    p_draws : (R, n, J) full draw array (for CPO / LPML)
    """
    n, R = X_c.shape[0], beta_draws.shape[0]
    J    = gamma_draws.shape[1] + 1
    p_draws = np.empty((R, n, J))

    for r in range(R):
        mu = X_c @ beta_draws[r]        # (n,)
        g  = gamma_draws[r]             # (J-1,)
        # cumulative CDF at each cut-point
        Fj = np.array([ndtr(g[j] - mu) for j in range(J - 1)]).T   # (n, J-1)
        p_draws[r, :, 0]        = Fj[:, 0]
        for j in range(1, J - 1):
            p_draws[r, :, j]    = Fj[:, j] - Fj[:, j - 1]
        p_draws[r, :, J - 1]   = 1.0 - Fj[:, J - 2]

    return p_draws.mean(axis=0), p_draws


def seq_pred_probs(X_s, alpha_draws):
    """
    Posterior predictive probabilities for the sequential model.

    Parameters
    ----------
    X_s         : (n, k)
    alpha_draws : list of J-1 arrays, each (R, k)

    Returns
    -------
    p_bar   : (n, J)
    p_draws : (R, n, J)
    """
    n  = X_s.shape[0]
    S  = len(alpha_draws)      # J - 1
    J  = S + 1
    R  = alpha_draws[0].shape[0]
    p_draws = np.empty((R, n, J))

    for r in range(R):
        # Phi(x'alpha_k) = P(continue at stage k)  ->  shape (n, S)
        phi_cont = np.column_stack(
            [ndtr( X_s @ alpha_draws[q][r]) for q in range(S)]
        )
        phi_exit = 1.0 - phi_cont

        # P(Y=0) = Phi(-x'alpha_1)
        p_draws[r, :, 0] = phi_exit[:, 0]
        # P(Y=j) = prod_{l=1}^{j} Phi(x'alpha_l) * Phi(-x'alpha_{j+1})
        for j in range(1, J - 1):
            p_draws[r, :, j] = phi_cont[:, :j].prod(axis=1) * phi_exit[:, j]
        # P(Y=J-1) = prod_{l=1}^{J-1} Phi(x'alpha_l)
        p_draws[r, :, J - 1] = phi_cont.prod(axis=1)

    return p_draws.mean(axis=0), p_draws


# ── LPML via harmonic-mean CPO ────────────────────────────────────────────────

def compute_lpml(y, p_draws):
    """
    LPML = sum_i log CPO_i,
    CPO_i = 1 / mean_r(1 / p(y_i | theta^r)).

    Parameters
    ----------
    y        : (n,) int
    p_draws  : (R, n, J)

    Returns
    -------
    lpml : float
    cpo  : (n,) array
    """
    y = np.asarray(y, int)
    R, n, J = p_draws.shape
    # lik[i, r] = p(y_i | theta^r)
    lik = p_draws[:, np.arange(n), y].T     # (n, R)
    lik = np.clip(lik, 1e-10, 1.0)
    cpo  = 1.0 / np.mean(1.0 / lik, axis=1)
    lpml = float(np.sum(np.log(np.clip(cpo, 1e-15, 1.0))))
    return lpml, cpo


# ── posterior summary ─────────────────────────────────────────────────────────

def posterior_summary(draws, param_names=None):
    """(R, p) draws -> DataFrame with mean, sd, q025, q975."""
    import pandas as pd
    p  = draws.shape[1]
    if param_names is None:
        param_names = [f'p{j}' for j in range(p)]
    lo = np.quantile(draws, 0.025, axis=0)
    hi = np.quantile(draws, 0.975, axis=0)
    return pd.DataFrame({
        'mean': draws.mean(0).round(4),
        'sd':   draws.std(0).round(4),
        'q025': lo.round(4),
        'q975': hi.round(4),
    }, index=param_names)

References

• Ashford, J. R. (1959). An approach to the analysis of data for semi-quantal responses in biological assay. Biometrics 15(4), 573–581. — the pneumoconiosis coal-miner data used throughout
• Albert, J. H. & Chib, S. (1993). Bayesian analysis of binary and polychotomous response data. Journal of the American Statistical Association 88(422), 669–679. — the data-augmentation Gibbs sampler for the cumulative ordered probit (Section 1)
• Albert, J. H. & Chib, S. (2001). Sequential ordinal modeling with applications to survival data. Biometrics 57(3), 829–836. — the sequential (stage-specific) ordinal probit of Section 2
• McCullagh, P. (1980). Regression models for ordinal data. Journal of the Royal Statistical Society: Series B 42(2), 109–142. — the cumulative-link model and the parallel-regressions (proportional-odds) assumption
• Yee, T. W. (2010). The VGAM package for categorical data analysis. Journal of Statistical Software 32(10), 1–34. — VGAM::cumulative / sratio, the R implementations of both links
• Liddell, T. M. & Kruschke, J. K. (2019). Analyzing ordinal data with metric models: what could possibly go wrong? Journal of Experimental Social Psychology 79, 328–348. — on treating ordinal outcomes properly rather than as metric
• Vehtari, A., Gelman, A. & Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing 27(5), 1413–1432. — the LOO ELPD / PSIS-LOO used for the bambi model comparison