Ordinal Probit Models — Cumulative vs Sequential

Pneumoconiosis data (Ashford 1959) — 371 coal miners, 8 exposure groups, J=3 severity levels

Data: pneumo.csv — exposure time (years) × {normal, mild, severe} counts. Source: Ashford (1959); classic example dataset in VGAM.

Section	Model	Inference
Section 1	Cumulative ordered probit (Albert & Chib 1993)	Gibbs sampler + statsmodels OrderedModel MLE
Section 2	Sequential ordinal probit (Albert & Chib 2001)	Gibbs sampler + sm.Probit × 2 MLE
Section 3	bambi / PyMC (HMC)	cumulative, sratio parallel (non-parallel cs() removed in bambi 0.18)
Section 4	Full comparison	LPML vs log-lik vs LOO-IC, slopes, key findings

Key question: Does exposure time affect getting sick (normal→disease) differently from progression (mild→severe)? The cumulative model imposes one effect; the sequential model estimates both freely.

---

How these models differ from ord_probit_gibbs.py

	ord_probit_gibbs	Cumulative (Section 1)	Sequential (Section 2)
Paper	Liddell & Kruschke (2019)	Albert & Chib (1993)	Albert & Chib (2001)
Purpose	Item quality/dispersion from ratings	Ordered regression on covariates	Stage-specific regression
Covariates	None — groups are the variable	$x_i'\beta$ (single effect)	$x_i'\alpha_k$ per stage
Variance $\sigma^2$	Estimated per group	Fixed = 1 (probit ID)	Fixed = 1 (probit ID)
Parameters	$\mu_g, \sigma_g, \tau_k$	$\beta, \gamma_1, \gamma_2$	$\alpha_1, \alpha_2$ (2D each)
Identification	Outer $\tau$ pinned at 1.5 / K−0.5	No intercept in $X_c$	Intercept absorbed in $\alpha_{k0}$
MLE equivalent	—	statsmodels.OrderedModel	Two sm.Probit fits
Bayesian package	—	bambi cumulative / brms	bambi sratio (parallel only in 0.18+) / brms
R equivalent	—	VGAM::cumulative / ordinal::clm	VGAM::sratio / brms sratio

References:

• Albert & Chib (1993), JASA 88, 669-679 — data-augmentation Gibbs, cumulative probit
• Albert & Chib (2001), JASA 96, 415-425 — sequential ordinal probit

import sys, warnings
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from scipy.special import ndtr
import statsmodels.api as sm
from statsmodels.miscmodels.ordinal_model import OrderedModel as OrdModel

warnings.filterwarnings('ignore')
sys.path.insert(0, r'C:\Users\user\project')

from seq_ord_probit_gibbs import (
    cumulative_probit_gibbs, sequential_probit_gibbs,
    cum_pred_probs, seq_pred_probs,
    compute_lpml, posterior_summary,
)

try:
    import bambi as bmb
    import arviz as az
    BAMBI_OK = True
    print(f'bambi {bmb.__version__} loaded')
except ImportError:
    BAMBI_OK = False
    print('bambi not installed — §3 cells will be skipped')
    print('Install with: conda install -c conda-forge bambi')

print('statsmodels', sm.__version__)
print('All imports OK')

bambi 0.18.0 loaded
statsmodels 0.14.6
All imports OK

DATA_PATH = r'C:\Users\user\project\pneumo.csv'
df = pd.read_csv(DATA_PATH)

# ── expand grouped counts to individual observations ─────────────────────────
rows = []
for _, row in df.iterrows():
    xt = row['exposure.time']
    for j, col in enumerate(['normal', 'mild', 'severe']):
        rows.extend([(xt, j)] * int(row[col]))

obs = pd.DataFrame(rows, columns=['exposure', 'y'])
n   = len(obs)
y   = obs['y'].values
x   = obs['exposure'].values

# standardise exposure for better Gibbs mixing
x_mean, x_sd = x.mean(), x.std()
x_std = (x - x_mean) / x_sd

# design matrices
X_c = x_std.reshape(-1, 1)                          # cumulative: no intercept
X_s = np.column_stack([np.ones(n), x_std])          # sequential: intercept + x

print(f'n={n}  J=3  x_mean={x_mean:.1f} yr  x_sd={x_sd:.1f} yr')
print()
print('Observed counts by group:')
header = f'{"exposure":>10}  {"normal":>8}  {"mild":>8}  {"severe":>8}  {"total":>8}  {"p(mild+sev)":>12}'
print(header)
for _, row in df.iterrows():
    tot = row['normal'] + row['mild'] + row['severe']
    print(f"{row['exposure.time']:10.1f}  {row['normal']:8.0f}  {row['mild']:8.0f}  {row['severe']:8.0f}  {tot:8.0f}  {(row['mild']+row['severe'])/tot:12.3f}")

n=371  J=3  x_mean=23.4 yr  x_sd=14.0 yr

Observed counts by group:
  exposure    normal      mild    severe     total   p(mild+sev)
       5.8        98         0         0        98         0.000
      15.0        51         2         1        54         0.056
      21.5        34         6         3        43         0.209
      27.5        35         5         8        48         0.271
      33.5        32        10         9        51         0.373
      39.5        23         7         8        38         0.395
      46.0        12         6        10        28         0.571
      51.5         4         2         5        11         0.636

---

Data: Ashford (1959) — Pneumoconiosis in British Coal Miners

Study design: Cross-sectional survey of 371 coal miners in South Wales collieries. Each miner's cumulative exposure time (years worked in coal mining) was recorded, and chest radiographs were graded by a radiologist into three ordered severity categories.

Disease categories (ordered):

• Normal (Y=0): No radiological signs of pneumoconiosis
• Mild (Y=1): Early radiological changes (category 1 or 2 in ILO classification)
• Severe (Y=2): Advanced disease (category 3+)

Key features of the data:

• Exposure ranges from 5.8 to 51.5 years — a nearly 10-fold range
• The disease gradient is stark: 0% diseased at ≤5.8 yr, 64% diseased at 51.5 yr
• Groups are unequal in size: n=98 at 5.8 yr (low exposure, mostly healthy) vs n=11 at 51.5 yr (small high-exposure cohort)
• Only 82 miners (22%) have any disease (mild or severe) — Stage-2 estimates are based entirely on this subset
• The mild category is rare throughout: at most 21% of any group, making it the hardest to fit

Why this is the canonical dataset for ordinal models: Ashford (1959) introduced this data to illustrate that a single continuous exposure variable (time) drives an ordered polytomous response. It has since become the standard benchmark for comparing cumulative vs sequential ordinal regression — see VGAM package documentation, Albert & Chib (2001), and Liddell & Kruschke (2019).

groups = df['exposure.time'].values
counts = df[['normal', 'mild', 'severe']].values.astype(float)
props  = counts / counts.sum(axis=1, keepdims=True)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# raw counts stacked bar
ax = axes[0]
colors = ['#4c9be8', '#f5a623', '#e84c4c']
labels = ['Normal', 'Mild', 'Severe']
bottom = np.zeros(len(groups))
for j, (lab, col) in enumerate(zip(labels, colors)):
    ax.bar(groups, counts[:, j], bottom=bottom, color=col, label=lab, width=2.5, edgecolor='white')
    bottom += counts[:, j]
ax.set_xlabel('Exposure time (years)')
ax.set_ylabel('Count')
ax.set_title('Observed counts')
ax.legend(loc='upper right', fontsize=9)

# proportions
ax = axes[1]
bottom = np.zeros(len(groups))
for j, (lab, col) in enumerate(zip(labels, colors)):
    ax.bar(groups, props[:, j], bottom=bottom, color=col, label=lab, width=2.5, edgecolor='white')
    bottom += props[:, j]
ax.set_xlabel('Exposure time (years)')
ax.set_ylabel('Proportion')
ax.set_title('Observed proportions')
ax.legend(loc='upper right', fontsize=9)

fig.suptitle('Pneumoconiosis data — Ashford (1959)  |  n=371 miners, 8 exposure groups', fontsize=11)
plt.tight_layout()
plt.show()

---

1. Cumulative Ordered Probit (Albert & Chib 1993)

Model: $$z_i^* \sim N(\beta_1 \tilde{x}_i,\, 1), \quad Y_i = j \iff \gamma_{j-1} < z_i^* \leq \gamma_j$$

where $\tilde{x}_i$ is standardised exposure time. Equivalently: $$P(Y_i \leq j \mid x_i) = \Phi(\gamma_j - \beta_1 \tilde{x}_i)$$

• Single slope $\beta_1$ affects ALL transitions equally — the parallel-regressions assumption (the probit analogue of proportional odds, a term that strictly names the cumulative logit model; the shorthand "proportional odds" is used loosely below).
• Design matrix $X_c$: no intercept — cut-points $\gamma_1 < \gamma_2$ absorb the intercept.
• Prior: $\beta_1 \sim N(0, 100)$, cut-points via uniform slice.

3-block Gibbs:

1. $z_i \mid \beta, \gamma, y_i \sim TN(\beta_1 \tilde{x}_i, 1;\, \gamma_{y_i-1}, \gamma_{y_i})$
2. $\beta \mid z \sim N(\bar{b}, \bar{B})$ — standard Normal posterior
3. $\gamma_j \mid z, y \sim \mathrm{Uniform}(\max_{y=j} z_i,\; \min_{y=j+1} z_i)$ — independent per $j$

out_cum = cumulative_probit_gibbs(
    y, X_c,
    b0=np.zeros(1), B0=np.array([[100.0]]),
    R=15000, burn=3000,
    seed=1, nprint=3000,
)

beta_c  = out_cum['beta']    # (12000, 1)
gamma_c = out_cum['gamma']   # (12000, 2)

  [cumulative]  iter   3000/15000  (0.1s)
  [cumulative]  iter   6000/15000  (0.2s)
  [cumulative]  iter   9000/15000  (0.3s)
  [cumulative]  iter  12000/15000  (0.5s)
  [cumulative]  iter  15000/15000  (0.6s)
[cumulative] Done in 0.6s  |  kept: 12000

cum_params = np.column_stack([beta_c[:, 0], gamma_c])
cum_names  = ['beta_1 (exposure, std)', 'gamma_1', 'gamma_2']
print('Cumulative ordered probit — posterior summary')
print(posterior_summary(cum_params, cum_names).to_string())

print()
print('beta_1 on raw scale (per year):',
      f"{beta_c[:, 0].mean() / x_sd:.4f}  (divide std-scale mean by x_sd={x_sd:.2f})")

Cumulative ordered probit — posterior summary
                          mean      sd    q025    q975
beta_1 (exposure, std)  0.7964  0.0960  0.6176  0.9939
gamma_1                 1.0133  0.1038  0.8314  1.2468
gamma_2                 1.5505  0.1198  1.3377  1.7940

beta_1 on raw scale (per year): 0.0567  (divide std-scale mean by x_sd=14.04)

Interpretation — Cumulative model parameters

Exposure slope β₁ = 0.796 (standardised) = 0.057 per year: Each additional year of mining shifts the latent disease propensity $z^*$ by 0.057 standard deviations upward. Over the full range of 45.7 years (5.8 → 51.5 yr), the total shift is $0.057 \times 45.7 \approx 2.6$ standard deviations — a large effect.

Cut-points γ₁ = 1.013, γ₂ = 1.551 (both positive, well above zero):

At the mean exposure ($\tilde{x}=0$, i.e. 23.4 yr): $$P(\text{normal}) = \Phi(1.013) \approx 84\%, \quad P(\text{mild}) = \Phi(1.551) - \Phi(1.013) \approx 9.5\%, \quad P(\text{severe}) \approx 6\%$$

At 5.8 yr ($\tilde{x}=-1.26$): $P(\text{normal}) = \Phi(1.013 + 1.003) = \Phi(2.016) \approx 97.8\%$ — nearly all healthy ✓

At 51.5 yr ($\tilde{x}=+2.01$): $P(\text{normal}) = \Phi(1.013 - 1.600) = \Phi(-0.587) \approx 27.9\%$ — 72% diseased

Proportional-odds assumption: The same β₁ governs both transitions simultaneously. This means, at every exposure level, the ratio of odds for P(Y≤0) to P(Y≤1) is constant. The fitted-probability table below will reveal whether this constraint holds in practice.

# predicted probs at the 8 observed exposure times
x_groups     = df['exposure.time'].values
x_groups_std = (x_groups - x_mean) / x_sd
X_c_pred     = x_groups_std.reshape(-1, 1)

p_bar_cum, p_draws_cum = cum_pred_probs(X_c_pred, beta_c, gamma_c)

print('Cumulative model — fitted P(Y=j) at each exposure group:')
print(f'{"exposure":>10}  {"obs_norm":>9}  {"fit_norm":>9}  {"obs_mild":>9}  {"fit_mild":>9}  {"obs_sev":>9}  {"fit_sev":>9}')
for g, row in df.iterrows():
    tot  = row['normal'] + row['mild'] + row['severe']
    on, om, os_ = row['normal']/tot, row['mild']/tot, row['severe']/tot
    fn, fm, fs  = p_bar_cum[g, 0], p_bar_cum[g, 1], p_bar_cum[g, 2]
    print(f"{row['exposure.time']:10.1f}  {on:9.3f}  {fn:9.3f}  {om:9.3f}  {fm:9.3f}  {os_:9.3f}  {fs:9.3f}")

Cumulative model — fitted P(Y=j) at each exposure group:
  exposure   obs_norm   fit_norm   obs_mild   fit_mild    obs_sev    fit_sev
       5.8      1.000      0.976      0.000      0.018      0.000      0.006
      15.0      0.944      0.930      0.037      0.047      0.019      0.023
      21.5      0.791      0.868      0.140      0.083      0.070      0.050
      27.5      0.729      0.782      0.104      0.123      0.167      0.095
      33.5      0.627      0.670      0.196      0.165      0.176      0.165
      39.5      0.605      0.540      0.184      0.197      0.211      0.263
      46.0      0.429      0.395      0.214      0.210      0.357      0.395
      51.5      0.364      0.284      0.182      0.200      0.455      0.517

---

Model-comparison metrics — CPO, LPML, and LOO elpd

Three related predictive scores appear from here on. All are leave-one-out measures — they ask how well the model predicts each observation when that observation is not used to fit it — so they reward fit while penalising overfitting, and all are on the log-predictive-density scale where higher (less negative) is better.

• CPO — Conditional Predictive Ordinate. For observation $i$, the predictive density of $y_i$ given all the other observations, $\text{CPO}_i = p(y_i \mid y_{-i})$. A small CPO flags a point the model predicts poorly. Estimated from the posterior draws by the harmonic mean of the in-sample likelihoods.
• LPML — Log Pseudo-Marginal Likelihood $= \sum_i \log \text{CPO}_i$. The CPOs summed into one leave-one-out score for the whole model; this is the Gibbs-sampler column below.
• LOO elpd — leave-one-out expected log predictive density. The same quantity computed by ArviZ's PSIS-LOO (az.loo) for the bambi/HMC fits. Because it is the log-predictive-density form (not the $\text{LOO-IC} = -2 \times \text{elpd}$ information-criterion form), it is directly comparable to the Gibbs LPML — and indeed the two match closely below (cumulative: LPML $-210.24$ vs LOO elpd $-210.11$).

So LPML and LOO elpd are the same idea from two samplers (Gibbs vs HMC), reported on the same scale; the in-sample MLE log-likelihood and AIC are shown alongside for contrast. A difference of only a nat or two between models is negligible — see the discussion in §4.

# LPML on the full n=371 individual observations
_, p_draws_cum_full = cum_pred_probs(X_c, beta_c, gamma_c)
lpml_cum, cpo_cum   = compute_lpml(y, p_draws_cum_full)
print(f'Cumulative model  LPML = {lpml_cum:.2f}')

Cumulative model  LPML = -210.24

# ── statsmodels OrderedModel — cumulative probit MLE ─────────────────────────
y_cat   = pd.Categorical(y, categories=[0, 1, 2], ordered=True)
res_ord = OrdModel(y_cat, pd.DataFrame({'x_std': x_std}),
                   distr='probit').fit(method='bfgs', disp=False)

beta_mle_cum = res_ord.params['x_std']
thresh_raw   = res_ord.params[1:]                              # log-diff parameterisation
thresh_all   = res_ord.model.transform_threshold_params(thresh_raw)
# transform_threshold_params prepends -inf as sentinel; actual thresholds at [1], [2], ...
gamma_mle_1, gamma_mle_2 = thresh_all[1], thresh_all[2]
loglik_mle_cum = res_ord.llf

print('statsmodels OrderedModel — cumulative probit MLE')
print(res_ord.summary().tables[1])
print()
print(f'Cut-points (actual scale):  gamma_1 = {gamma_mle_1:.4f}   gamma_2 = {gamma_mle_2:.4f}')
print()
print('--- Gibbs vs MLE (cumulative) ---')
print(f'  beta_1   Gibbs = {beta_c[:,0].mean():.4f}   MLE = {beta_mle_cum:.4f}')
print(f'  gamma_1  Gibbs = {gamma_c[:,0].mean():.4f}   MLE = {gamma_mle_1:.4f}')
print(f'  gamma_2  Gibbs = {gamma_c[:,1].mean():.4f}   MLE = {gamma_mle_2:.4f}')
print(f'  log-lik  LPML  = {lpml_cum:.2f}   MLE log-lik = {loglik_mle_cum:.2f}')

statsmodels OrderedModel — cumulative probit MLE
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
x_std          0.7852      0.093      8.481      0.000       0.604       0.967
0/1            1.0055      0.094     10.665      0.000       0.821       1.190
1/2           -0.6493      0.150     -4.328      0.000      -0.943      -0.355
==============================================================================

Cut-points (actual scale):  gamma_1 = 1.0055   gamma_2 = 1.5279

--- Gibbs vs MLE (cumulative) ---
  beta_1   Gibbs = 0.7964   MLE = 0.7852
  gamma_1  Gibbs = 1.0133   MLE = 1.0055
  gamma_2  Gibbs = 1.5505   MLE = 1.5279
  log-lik  LPML  = -210.24   MLE log-lik = -207.17

---

2. Sequential Ordinal Probit (Albert & Chib 2001)

Key idea: Each transition is a separate binary probit with its own parameter vector.

Stage 1 (all $n$ subjects — normal vs any disease): $$z_{i1} \sim N(\alpha_{10} + \alpha_{11}\tilde{x}_i,\; 1), \quad s_{i1} = \mathbf{1}(z_{i1} > 0)$$ $s_{i1}=0 \Rightarrow Y_i = 0$ (normal); $s_{i1}=1 \Rightarrow$ enter Stage 2.

Stage 2 (subjects with $Y_i \geq 1$ only — mild vs severe): $$z_{i2} \sim N(\alpha_{20} + \alpha_{21}\tilde{x}_i,\; 1), \quad s_{i2} = \mathbf{1}(z_{i2} > 0)$$ $s_{i2}=0 \Rightarrow Y_i=1$ (mild); $s_{i2}=1 \Rightarrow Y_i=2$ (severe).

Joint probabilities: $$P(Y=0\mid x) = \Phi(-x'\alpha_1), \quad P(Y=1\mid x) = \Phi(x'\alpha_1)\,\Phi(-x'\alpha_2), \quad P(Y=2\mid x) = \Phi(x'\alpha_1)\,\Phi(x'\alpha_2)$$

Stages are conditionally independent given data: two separate Albert & Chib (1993) binary probit Gibbs runs. Prior: $\alpha_k \sim N(0,\, 100\cdot I_2)$ at each stage.

out_seq = sequential_probit_gibbs(
    y, X_s,
    b0=np.zeros(2), B0=100.0 * np.eye(2),
    R=15000, burn=3000,
    seed=2, nprint=3000,
)

alpha1 = out_seq['alpha'][0]   # (12000, 2)  stage 1: normal → disease
alpha2 = out_seq['alpha'][1]   # (12000, 2)  stage 2: mild → severe

at_risk1 = out_seq['at_risk'][0]   # all n
at_risk2 = out_seq['at_risk'][1]   # only y >= 1
print(f'Stage 1 at-risk: {at_risk1.sum()} (all)   Stage 2 at-risk: {at_risk2.sum()} (mild + severe)')

  [sequential]  iter   3000/15000  (0.1s)
  [sequential]  iter   6000/15000  (0.2s)
  [sequential]  iter   9000/15000  (0.4s)
  [sequential]  iter  12000/15000  (0.5s)
  [sequential]  iter  15000/15000  (0.6s)
[sequential] Done in 0.6s  |  kept: 12000
Stage 1 at-risk: 371 (all)   Stage 2 at-risk: 82 (mild + severe)

print('Sequential model — Stage 1  (normal → any disease)')
print(posterior_summary(alpha1, ['alpha_1_0 (intercept)', 'alpha_1_1 (exposure, std)']).to_string())
print()
print('Sequential model — Stage 2  (mild → severe)')
print(posterior_summary(alpha2, ['alpha_2_0 (intercept)', 'alpha_2_1 (exposure, std)']).to_string())

print()
print('--- Exposure slopes on raw scale (per year) ---')
print(f'  Stage 1  alpha_1_1 / x_sd = {alpha1[:, 1].mean() / x_sd:.4f}')
print(f'  Stage 2  alpha_2_1 / x_sd = {alpha2[:, 1].mean() / x_sd:.4f}')

Sequential model — Stage 1  (normal → any disease)
                             mean      sd    q025    q975
alpha_1_0 (intercept)     -1.0165  0.0981 -1.2121 -0.8295
alpha_1_1 (exposure, std)  0.8041  0.1000  0.6118  1.0063

Sequential model — Stage 2  (mild → severe)
                             mean      sd    q025    q975
alpha_2_0 (intercept)     -0.1722  0.2241 -0.6071  0.2658
alpha_2_1 (exposure, std)  0.3131  0.2040 -0.0802  0.7180

--- Exposure slopes on raw scale (per year) ---
  Stage 1  alpha_1_1 / x_sd = 0.0573
  Stage 2  alpha_2_1 / x_sd = 0.0223

Interpretation — Sequential model parameters

Stage 1 — Who gets sick? (n = 371)

$\alpha_{10} = -1.017$, $\alpha_{11} = 0.804$ (std) $= 0.057$ per year.

At mean exposure ($\tilde{x}=0$): $P(\text{any disease}) = \Phi(-1.017) \approx 15\%$. At 51.5 yr: $P(\text{any disease}) = \Phi(-1.017 + 0.804 \times 2.01) = \Phi(0.599) \approx 72\%$.

The Stage-1 slope (0.804) is virtually identical to the cumulative β₁ (0.796) — both models agree completely on the onset transition. The 95% CrI [0.61, 1.01] firmly excludes zero.

Stage 2 — Who progresses to severe? (n = 82 diseased only)

$\alpha_{20} = -0.172$, $\alpha_{21} = 0.313$ (std) $= 0.022$ per year.

At mean exposure: $P(\text{severe} \mid \text{diseased}) = \Phi(-0.172) \approx 43\%$. At 51.5 yr: $P(\text{severe} \mid \text{diseased}) = \Phi(-0.172 + 0.313 \times 2.01) = \Phi(0.457) \approx 68\%$.

Key finding: The Stage-2 slope (0.313) is 2.6× smaller than Stage-1 (0.804), and its 95% CrI [−0.08, 0.72] includes zero. This means:

• Once a miner has any disease, further years of exposure add much less additional risk of progression to severe
• The proportional-odds assumption (common β for both transitions) is violated — the cumulative model is applying a β₁ ≈ 0.796 to Stage 2 where the true effect is closer to 0.31
• $P(\alpha_{11} > \alpha_{21}) = 0.984$: near-certain that Stage-1 effect exceeds Stage-2

Epidemiological interpretation: Total exposure time is primarily a determinant of whether a miner develops pneumoconiosis, not of how quickly mild disease progresses to severe. Individual factors not captured here (dust type, ventilation, individual susceptibility) likely drive progression.

X_s_pred = np.column_stack([np.ones(len(x_groups)), x_groups_std])
p_bar_seq, p_draws_seq = seq_pred_probs(X_s_pred, [alpha1, alpha2])

print('Sequential model — fitted P(Y=j) at each exposure group:')
print(f'{"exposure":>10}  {"obs_norm":>9}  {"fit_norm":>9}  {"obs_mild":>9}  {"fit_mild":>9}  {"obs_sev":>9}  {"fit_sev":>9}')
for g, row in df.iterrows():
    tot  = row['normal'] + row['mild'] + row['severe']
    on, om, os_ = row['normal']/tot, row['mild']/tot, row['severe']/tot
    fn, fm, fs  = p_bar_seq[g, 0], p_bar_seq[g, 1], p_bar_seq[g, 2]
    print(f"{row['exposure.time']:10.1f}  {on:9.3f}  {fn:9.3f}  {om:9.3f}  {fm:9.3f}  {os_:9.3f}  {fs:9.3f}")

Sequential model — fitted P(Y=j) at each exposure group:
  exposure   obs_norm   fit_norm   obs_mild   fit_mild    obs_sev    fit_sev
       5.8      1.000      0.977      0.000      0.016      0.000      0.007
      15.0      0.944      0.931      0.037      0.044      0.019      0.025
      21.5      0.791      0.869      0.140      0.077      0.070      0.055
      27.5      0.729      0.782      0.104      0.116      0.167      0.102
      33.5      0.627      0.669      0.196      0.159      0.176      0.172
      39.5      0.605      0.538      0.184      0.197      0.211      0.265
      46.0      0.429      0.392      0.214      0.227      0.357      0.381
      51.5      0.364      0.280      0.182      0.238      0.455      0.482

_, p_draws_seq_full = seq_pred_probs(X_s, [alpha1, alpha2])
lpml_seq, cpo_seq   = compute_lpml(y, p_draws_seq_full)
print(f'Sequential model  LPML = {lpml_seq:.2f}')

Sequential model  LPML = -211.00

# ── statsmodels sm.Probit × 2 — sequential MLE ───────────────────────────────
# Stage 1: all n subjects, outcome = any disease (y >= 1)
y1_mle = (y >= 1).astype(float)
X1_mle = np.column_stack([np.ones(n), x_std])
res_s1 = sm.Probit(y1_mle, X1_mle).fit(disp=False)

# Stage 2: only at-risk subjects (y >= 1), outcome = severe (y == 2)
mask2  = y >= 1
y2_mle = (y[mask2] == 2).astype(float)
X2_mle = np.column_stack([np.ones(mask2.sum()), x_std[mask2]])
res_s2 = sm.Probit(y2_mle, X2_mle).fit(disp=False)

loglik_mle_seq = res_s1.llf + res_s2.llf

print('Sequential probit MLE — Stage 1 (normal → disease, n=371):')
print(f'  intercept = {res_s1.params[0]:.4f}   slope = {res_s1.params[1]:.4f}   log-lik = {res_s1.llf:.2f}')
print(f'  Gibbs:  alpha_10 = {alpha1[:,0].mean():.4f}   alpha_11 = {alpha1[:,1].mean():.4f}')
print()
print(f'Sequential probit MLE — Stage 2 (mild → severe, n={int(mask2.sum())}):')
print(f'  intercept = {res_s2.params[0]:.4f}   slope = {res_s2.params[1]:.4f}   log-lik = {res_s2.llf:.2f}')
print(f'  Gibbs:  alpha_20 = {alpha2[:,0].mean():.4f}   alpha_21 = {alpha2[:,1].mean():.4f}')
print()
print(f'Sequential total log-lik (stage1 + stage2) = {loglik_mle_seq:.2f}')
print(f'Cumulative MLE log-lik                      = {loglik_mle_cum:.2f}')
print(f'Delta (seq − cum)                           = {loglik_mle_seq - loglik_mle_cum:+.2f}')

Sequential probit MLE — Stage 1 (normal → disease, n=371):
  intercept = -1.0094   slope = 0.7976   log-lik = -151.61
  Gibbs:  alpha_10 = -1.0165   alpha_11 = 0.8041

Sequential probit MLE — Stage 2 (mild → severe, n=82):
  intercept = -0.1771   slope = 0.3115   log-lik = -55.44
  Gibbs:  alpha_20 = -0.1722   alpha_21 = 0.3131

Sequential total log-lik (stage1 + stage2) = -207.05
Cumulative MLE log-lik                      = -207.17
Delta (seq − cum)                           = +0.12

---

3. bambi / PyMC (HMC)

Two bambi models, both using PyMC/NUTS as the backend:

bambi model	Equivalent to	Notes
cumulative("probit")	Section 1 Gibbs cumulative	Same model, HMC vs data-augmentation Gibbs
sratio("probit")	Parallel stopping-ratio — different from Section 2	Single β across all stages (proportional hazards)

Non-parallel sratio (stage-specific slopes, exact equivalent of Section 2 Gibbs) used cs(x_std) in older bambi versions. That syntax was removed in bambi 0.18.0 / formulae 0.6.2. The Section 2 Gibbs sequential sampler is the working non-parallel implementation.

Sign convention: In bambi sratio, positive x_std means higher exposure → more likely to continue (not stop) at each stage → more disease, same direction as our $\alpha_{k1} > 0$.

print('=' * 80)
print(f'{"§4 — Full comparison  (all methods)":^80}')
print('=' * 80)

# ── Gibbs LPML ──
print(f'\nBayesian (Gibbs)   LPML — cumulative  : {lpml_cum:.2f}')
print(f'Bayesian (Gibbs)   LPML — sequential  : {lpml_seq:.2f}  (delta {lpml_seq - lpml_cum:+.2f})')

# ── MLE log-lik ──
print(f'\nFrequentist (MLE)  log-lik — cumulative : {loglik_mle_cum:.2f}')
print(f'Frequentist (MLE)  log-lik — sequential : {loglik_mle_seq:.2f}  (delta {loglik_mle_seq - loglik_mle_cum:+.2f})')

# ── exposure slopes ──
print('\n--- Exposure slope (standardised x) — all methods ---')
print(f'  Cumulative  Gibbs    beta_1    : {beta_c[:,0].mean():.4f}  '
      f'95% CrI [{np.quantile(beta_c[:,0], 0.025):.3f}, {np.quantile(beta_c[:,0], 0.975):.3f}]')
print(f'  Cumulative  MLE      beta_1    : {beta_mle_cum:.4f}  '
      f'95% CI  [{beta_mle_cum - 1.96*res_ord.bse["x_std"]:.3f}, {beta_mle_cum + 1.96*res_ord.bse["x_std"]:.3f}]')
print(f'  Sequential  Gibbs    alpha_11  : {alpha1[:,1].mean():.4f}  '
      f'95% CrI [{np.quantile(alpha1[:,1], 0.025):.3f}, {np.quantile(alpha1[:,1], 0.975):.3f}]')
print(f'  Sequential  Gibbs    alpha_21  : {alpha2[:,1].mean():.4f}  '
      f'95% CrI [{np.quantile(alpha2[:,1], 0.025):.3f}, {np.quantile(alpha2[:,1], 0.975):.3f}]')
print(f'  Sequential  MLE stg1 alpha_11  : {res_s1.params[1]:.4f}  '
      f'95% CI  [{res_s1.params[1] - 1.96*res_s1.bse[1]:.3f}, {res_s1.params[1] + 1.96*res_s1.bse[1]:.3f}]')
print(f'  Sequential  MLE stg2 alpha_21  : {res_s2.params[1]:.4f}  '
      f'95% CI  [{res_s2.params[1] - 1.96*res_s2.bse[1]:.3f}, {res_s2.params[1] + 1.96*res_s2.bse[1]:.3f}]')

diff_draws = alpha1[:, 1] - alpha2[:, 1]
diff_mle   = res_s1.params[1] - res_s2.params[1]
print(f'\nP(alpha_11 > alpha_21)  Gibbs = {(diff_draws > 0).mean():.3f}  |  MLE diff = {diff_mle:+.4f}')

================================================================================
                      §4 — Full comparison  (all methods)                       
================================================================================

Bayesian (Gibbs)   LPML — cumulative  : -210.24
Bayesian (Gibbs)   LPML — sequential  : -211.00  (delta -0.77)

Frequentist (MLE)  log-lik — cumulative : -207.17
Frequentist (MLE)  log-lik — sequential : -207.05  (delta +0.12)

--- Exposure slope (standardised x) — all methods ---
  Cumulative  Gibbs    beta_1    : 0.7964  95% CrI [0.618, 0.994]
  Cumulative  MLE      beta_1    : 0.7852  95% CI  [0.604, 0.967]
  Sequential  Gibbs    alpha_11  : 0.8041  95% CrI [0.612, 1.006]
  Sequential  Gibbs    alpha_21  : 0.3131  95% CrI [-0.080, 0.718]
  Sequential  MLE stg1 alpha_11  : 0.7976  95% CI  [0.608, 0.987]
  Sequential  MLE stg2 alpha_21  : 0.3115  95% CI  [-0.090, 0.713]

P(alpha_11 > alpha_21)  Gibbs = 0.984  |  MLE diff = +0.4862

---

4. Key Findings and Analysis

Finding 1: All inference methods agree on the exposure effect

The cumulative β₁ estimates across three independent inference methods:

Method	Cumulative β₁	Sequential α₁₁	Sequential α₂₁
Gibbs (data augmentation)	0.796	0.804	0.313
MLE (BFGS)	0.785	0.798	0.312
bambi HMC (NUTS)	0.780	—	—

All three methods agree to within 0.016 on the cumulative slope and within 0.006 on both sequential slopes. This confirms the implementations are correct and the results are robust to inference method.

Finding 2: The proportional-odds assumption is violated

$P(\alpha_{11} > \alpha_{21}) = 0.984$ (Gibbs). The exposure slope at Stage 1 (onset, full cohort n=371) is $\alpha_{11} = 0.804$ — nearly 2.6× larger than at Stage 2 (progression, at-risk n=82): $\alpha_{21} = 0.313$. The cumulative model's single $\beta_1 = 0.796$ is a weighted average of these two stage-specific effects: it correctly estimates onset risk but overstates the effect on progression from mild to severe.

The Stage-2 CrI $[-0.08,\, 0.72]$ includes zero: the data do not conclusively establish a progression effect of exposure. This is primarily a sample-size problem — only 82 of 371 miners are at risk at Stage 2.

Finding 3: All model-fit metrics are coherent

Metric	Cumulative	Sequential / sratio	Δ (seq − cum)	Favours
Gibbs LPML	−210.24	−211.00 (sequential)	−0.77	Cumulative
MLE log-lik	−207.17	−207.05 (sequential)	+0.12	Sequential (marginally)
bambi LOO elpd	−210.11	−212.30 (sratio ∥)	−2.19	Cumulative
AIC	420.3	422.1 (sequential)	+1.76	Cumulative

These are not contradictory:

• MLE log-lik rewards more parameters in-sample — the extra sequential parameter always helps slightly. LR test: $\chi^2 = 2 \times 0.12 = 0.24$ on 1 df, $p \approx 0.62$ — no significant improvement.
• LPML and the LOO elpd penalise parameters estimated with high uncertainty. The Stage-2 slope ($\alpha_{21}$, sd=0.20, CrI crosses 0) shifts when any one diseased miner is held out, reducing CPO for Stage-2 observations. The sequential model pays an out-of-sample penalty for the flexibility it cannot yet support with n=82.
• AIC = $-2 \log L + 2k$: $-2 \times 0.12 + 2 \times 1 = +1.76$ — consistent with LPML.

The data are just strong enough to detect the heterogeneity ($P \approx 0.98$) but not strong enough that the sequential model's extra parameter pays off out-of-sample.

Finding 4: Scientific vs predictive value

Both models produce nearly identical fitted probabilities at every exposure group (max difference ≈ 0.04, at the highest-exposure group where the two links diverge most; below 0.02 at every other group). The sequential model's contribution is scientific, not predictive: it shows that exposure time is primarily a determinant of who gets sick (Stage 1, large and well-estimated effect) rather than who progresses (Stage 2, smaller and uncertain effect). The cumulative model structurally cannot make this distinction.

Epidemiological implication: Reducing dust exposure would primarily prevent disease onset. Among miners who already have mild pneumoconiosis, additional exposure years add much less incremental risk of progression to severe — though the evidence is too weak to rule this out definitively. A study focused on the 82 already-diseased miners (or a larger cohort) would be needed to resolve the Stage-2 question.

fig = plt.figure(figsize=(15, 10))
gs  = gridspec.GridSpec(2, 3, figure=fig, hspace=0.4, wspace=0.35)

colors = ['#4c9be8', '#f5a623', '#e84c4c']
labels = ['Normal', 'Mild', 'Severe']

# ── top row: fitted prob curves vs observed ───────────────────────────────────
x_fine     = np.linspace(x_groups.min() - 2, x_groups.max() + 2, 200)
x_fine_std = (x_fine - x_mean) / x_sd

Xc_fine = x_fine_std.reshape(-1, 1)
Xs_fine = np.column_stack([np.ones(len(x_fine)), x_fine_std])

p_fine_cum, _ = cum_pred_probs(Xc_fine, beta_c, gamma_c)
p_fine_seq, _ = seq_pred_probs(Xs_fine, [alpha1, alpha2])

for j, (lab, col) in enumerate(zip(labels, colors)):
    ax = fig.add_subplot(gs[0, j])
    # observed proportions
    tot_g  = df[['normal','mild','severe']].sum(axis=1).values
    obs_j  = df[['normal','mild','severe']].values[:, j] / tot_g
    ax.scatter(x_groups, obs_j, s=50, color=col, zorder=5, label='Observed')
    # fitted curves
    ax.plot(x_fine, p_fine_cum[:, j], 'k-',  lw=2,   label='Cumulative')
    ax.plot(x_fine, p_fine_seq[:, j], 'k--', lw=1.5, label='Sequential')
    ax.set_title(f'P(Y = {j} = {lab})', fontsize=10)
    ax.set_xlabel('Exposure time (years)')
    ax.set_ylabel('Probability')
    ax.set_ylim(-0.05, 1.05)
    ax.legend(fontsize=8)

# ── bottom-left: posterior of exposure slopes ─────────────────────────────────
ax = fig.add_subplot(gs[1, 0])
ax.hist(beta_c[:, 0],  bins=60, density=True, alpha=0.6, color='steelblue', label='Cumulative $\\beta_1$')
ax.hist(alpha1[:, 1],  bins=60, density=True, alpha=0.6, color='darkorange', label='Seq. Stage-1 $\\alpha_{11}$')
ax.hist(alpha2[:, 1],  bins=60, density=True, alpha=0.6, color='firebrick',  label='Seq. Stage-2 $\\alpha_{21}$')
ax.axvline(0, color='black', lw=1, ls=':')
ax.set_xlabel('Exposure slope (standardised)')
ax.set_ylabel('Density')
ax.set_title('Posterior of exposure effect')
ax.legend(fontsize=8)

# ── bottom-middle: posterior of (alpha_11 - alpha_21) ─────────────────────────
ax = fig.add_subplot(gs[1, 1])
ax.hist(diff_draws, bins=60, density=True, color='purple', alpha=0.7)
ax.axvline(0,              color='black', lw=1.5, ls='--', label='0')
ax.axvline(diff_draws.mean(), color='red', lw=1.5, label=f'mean={diff_draws.mean():.3f}')
ax.set_xlabel('$\\alpha_{11} - \\alpha_{21}$  (stage-1 minus stage-2 slope)')
ax.set_ylabel('Density')
ax.set_title('Stage-1 vs Stage-2 exposure effect')
ax.legend(fontsize=9)

# ── bottom-right: CPO scatter ─────────────────────────────────────────────────
ax = fig.add_subplot(gs[1, 2])
scatter_c = ax.scatter(np.log(cpo_cum), np.log(cpo_seq),
                       c=y, cmap='Set1', s=8, alpha=0.5)
mn = min(np.log(cpo_cum).min(), np.log(cpo_seq).min()) - 0.1
mx = max(np.log(cpo_cum).max(), np.log(cpo_seq).max()) + 0.1
ax.plot([mn, mx], [mn, mx], 'k--', lw=1, label='Equal')
ax.set_xlabel('log CPO — Cumulative')
ax.set_ylabel('log CPO — Sequential')
ax.set_title('CPO comparison (above diagonal = seq better)')
cbar = plt.colorbar(scatter_c, ax=ax, ticks=[0, 1, 2])
cbar.set_ticklabels(['Normal', 'Mild', 'Severe'])
ax.legend(fontsize=8)

fig.suptitle('Cumulative vs Sequential Ordered Probit — Pneumoconiosis data', fontsize=12, y=1.01)
plt.savefig(r'C:\Users\user\project\ordinal_comparison.png',
            dpi=120, bbox_inches='tight')
plt.show()

# ── bambi — Bayesian inference via HMC (PyMC/NUTS backend) ───────────────────
if not BAMBI_OK:
    print("bambi not installed. Run:  conda install -c conda-forge bambi")
else:
    import pymc as pm

    def _loo_elpd(r):
        """Return (elpd, se) from ELPDData — handles ArviZ 0.x (elpd_loo) and 1.x (elpd)."""
        elpd = getattr(r, 'elpd_loo', None)
        if elpd is None:
            elpd = getattr(r, 'elpd', float('nan'))
        se = getattr(r, 'se', float('nan'))
        return float(elpd), float(se)

    df_bmb = pd.DataFrame({
        'y'     : pd.Categorical(y, categories=[0, 1, 2], ordered=True),
        'x_std' : x_std,
    })

    # 1. Cumulative probit (parallel) — same model as §1 Gibbs
    print("Fitting bambi cumulative probit ...")
    m_cum_bmb = bmb.Model("y ~ x_std", df_bmb, family="cumulative", link="probit")
    idata_cum_bmb = m_cum_bmb.fit(
        draws=2000, tune=1000, random_seed=20, target_accept=0.9, progressbar=False,
    )

    # 2. sratio parallel — single β across all stages (proportional hazards)
    print("Fitting bambi sratio parallel ...")
    m_sra_bmb = bmb.Model("y ~ x_std", df_bmb, family="sratio", link="probit")
    idata_sra_bmb = m_sra_bmb.fit(
        draws=2000, tune=1000, random_seed=21, target_accept=0.9, progressbar=False,
    )

    # Compute log-likelihood post-hoc (PyMC 6 / bambi 0.18 don't store it by default)
    try:
        with m_cum_bmb.backend.model:
            pm.compute_log_likelihood(idata_cum_bmb)
        with m_sra_bmb.backend.model:
            pm.compute_log_likelihood(idata_sra_bmb)

        loo_cum_bmb = az.loo(idata_cum_bmb)
        loo_sra_bmb = az.loo(idata_sra_bmb)

        elpd_c, se_c = _loo_elpd(loo_cum_bmb)
        elpd_s, se_s = _loo_elpd(loo_sra_bmb)

        print(f"\n--- bambi LOO elpd (higher = better; same scale as the Gibbs LPML) ---")
        print(f"  cumulative (parallel)   elpd = {elpd_c:.2f}   se={se_c:.2f}")
        print(f"  sratio parallel         elpd = {elpd_s:.2f}   se={se_s:.2f}")
        print(f"  (Gibbs LPML: cumulative={lpml_cum:.2f},  sequential={lpml_seq:.2f})")
    except Exception as e:
        print(f"\nLOO not available ({e})")
        print("  LOO comparison already available from Gibbs LPML in §4.")

    # ── Slope comparison: Gibbs vs bambi vs MLE ───────────────────────────────
    b_cum_bmb = idata_cum_bmb.posterior['x_std'].values.ravel()
    b_sra_bmb = idata_sra_bmb.posterior['x_std'].values.ravel()

    print(f"\n--- Cumulative exposure slope (std x) ---")
    print(f"  Gibbs   posterior mean = {beta_c[:,0].mean():.4f}   sd = {beta_c[:,0].std():.4f}")
    print(f"  bambi   posterior mean = {b_cum_bmb.mean():.4f}   sd = {b_cum_bmb.std():.4f}")
    print(f"  MLE                   = {beta_mle_cum:.4f}")

    print(f"\n--- sratio parallel slope (std x) ---")
    print(f"  bambi sratio mean = {b_sra_bmb.mean():.4f}   sd = {b_sra_bmb.std():.4f}")
    print(f"  (Gibbs stage-1 α₁₁ = {alpha1[:,1].mean():.4f},  stage-2 α₂₁ = {alpha2[:,1].mean():.4f})")
    print(f"  sratio parallel β sits between the two Gibbs stage-specific slopes")

Initializing NUTS using jitter+adapt_diag...

Fitting bambi cumulative probit ...

Multiprocess sampling (4 chains in 4 jobs)
NUTS: [threshold, x_std]
Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 5 seconds.
Initializing NUTS using jitter+adapt_diag...

Fitting bambi sratio parallel ...

Multiprocess sampling (4 chains in 4 jobs)
NUTS: [threshold, x_std]
Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 5 seconds.

Output()

Output()

--- bambi LOO elpd (higher = better; same scale as the Gibbs LPML) ---
  cumulative (parallel)   elpd = -210.11   se=14.37
  sratio parallel         elpd = -212.30   se=14.79
  (Gibbs LPML: cumulative=-210.24,  sequential=-211.00)

--- Cumulative exposure slope (std x) ---
  Gibbs   posterior mean = 0.7964   sd = 0.0960
  bambi   posterior mean = 0.7800   sd = 0.0905
  MLE                   = 0.7852

--- sratio parallel slope (std x) ---
  bambi sratio mean = 0.7121   sd = 0.0847
  (Gibbs stage-1 α₁₁ = 0.8041,  stage-2 α₂₁ = 0.3131)
  sratio parallel β sits between the two Gibbs stage-specific slopes

---

Summary — All methods

Parameter estimates (exposure slope, standardised x)

Method	Estimator	Cumulative β₁	Stage-1 α₁₁	Stage-2 α₂₁
Gibbs (AC 1993/2001)	Posterior mean	0.796	0.804	0.313
statsmodels MLE	MLE	0.785	0.798	0.312
bambi cumulative (HMC)	Posterior mean	0.780	—	—
bambi sratio parallel (HMC)	Posterior mean	0.712 (single β, see Section 3)	—	—

All three independent inference methods agree on the cumulative β₁ to within 0.016 and on both sequential slopes to within 0.006.

Model fit

Metric	Cumulative	Sequential / sratio ∥	Δ (seq − cum)	Favours
Gibbs LPML	−210.24	−211.00	−0.77	Cumulative
MLE log-lik	−207.17	−207.05	+0.12	Sequential (marginal)
bambi LOO elpd	−210.11	−212.30 (sratio ∥)	−2.19	Cumulative
AIC	420.3	422.1	+1.76	Cumulative

Gibbs vs MLE vs HMC — what to expect

• Gibbs ≈ MLE at n=371 with a flat prior: posterior mean ≈ MLE, LPML ≈ log-lik
• Gibbs ≈ bambi (same model, different sampler): cumulative Gibbs (0.796) vs bambi cumulative (0.780) are near-identical; small differences are Monte Carlo noise
• sratio parallel ≠ sequential Gibbs: the parallel model forces one β across stages; the Section 2 Gibbs has stage-specific $\alpha_k$. Non-parallel sratio(cs()) was removed in bambi 0.18 — the Section 2 Gibbs sampler is the working non-parallel implementation.

Stage-effect interpretation

• $\alpha_{11} > 0$: exposure drives onset of disease (normal → any)
• $\alpha_{21} > 0$: exposure drives progression (mild → severe)
• If $\alpha_{11} \gg \alpha_{21}$: exposure is primarily an initiating factor
• If $\alpha_{11} \approx \alpha_{21}$: proportional-odds (cumulative model) is adequate
• This dataset: $\alpha_{11} = 0.804 \gg \alpha_{21} = 0.313$, $P(\alpha_{11} > \alpha_{21}) = 0.984$ — exposure initiates disease; progression appears to depend on other factors not measured here

Computational notes

	Gibbs	MLE	bambi/HMC
Convergence	Fast (data augmentation, closed-form blocks)	BFGS, instant	NUTS, slower
Latent z	Returned (enables AC1995 residuals)	Not available	Not available
LPML	Direct from posterior draws	Requires bootstrap	LOO elpd from ArviZ
Stage independence	Explicit (separate Gibbs per stage)	Explicit (separate fits)	Joint (NUTS samples all at once)