Limited Dependent Variables
Models for outcomes that are discrete or censored rather than continuous, threaded by one idea — a latent variable and data augmentation: introduce the unobserved utility (or censored value) and the awkward probit / logit / Tobit likelihood collapses into an ordinary conjugate regression. Binary choice runs deepest: the probit (Albert–Chib) and logit (Pólya–Gamma) samplers, then their extensions — hierarchical (varying intercepts & slopes), random-effects panel logit, robust (robit / t-link), and a mixture-of-normals probit that estimates the link itself. Ordinal outcomes get the ordered and sequential probit (proportional-odds and its violations). Unordered choice spans the multivariate and multinomial probit and the hierarchical multinomial logit — with a mixture-heterogeneity variant — tracing how IIA is relaxed within a choice versus across people. And the Tobit closes the set with censored Gaussian regression. Each is built from scratch and cross-checked across PyMC and R.
Bayesian Binary Probit
A Python implementation of the Albert & Chib (1993) data-augmentation Gibbs
sampler for the binary probit model. Conditioning on latent utilities
zi reduces the probit likelihood to a standard Normal
regression, making β conjugate and the posterior precision
matrix constant across iterations — computed once via Cholesky, making the
sampler very fast. Implements Albert & Chib (1995) Bayesian residual
diagnostics: CPO (harmonic mean estimator) and latent residuals
ēi — a diagnostic with no classical analogue.
Applied to synthetic data (n=400, 25 label-flipped outliers): Section 1 Gibbs vs MLE,
Section 2 latent-residual and CPO outlier detection, Section 3 a high-leverage case where
CPO detects 4/5 planted outliers while Pearson detects 0/5. Includes an R
comparison via bayesm::rbprobitGibbs.
Bayesian Binary Logit
A Python implementation of the Bayesian binary logit model via two from-scratch
samplers — a tuned random-walk Metropolis (inverse-Fisher proposal,
Roberts–Gelman–Gilks 2.38/√k scaling) and a Pólya–Gamma
data-augmentation Gibbs sampler (Polson, Scott & Windle 2013). Because the
logit posterior has no conjugate form, Pólya–Gamma augmentation introduces a latent
ωi ~ PG(1, xi'β) that makes β conditionally
Gaussian — the exact logit analog of Albert & Chib for probit, and why
bayesm does probit by Gibbs but logit by Metropolis. Both samplers are
cross-checked against bambi, PyMC (NUTS), R's bayesm
(rmnlIndepMetrop), and the MLE — six routes, two languages, in agreement.
Applied to the canonical Mroz (1987) women's labor-force participation data
(n=753): a child under 6 multiplies the odds of participation by
e−1.44 ≈ 0.23, with odds ratios, predicted probabilities,
and average partial effects each carrying a credible interval.
Bayesian Hierarchical Binary Logit
A Python implementation of the multilevel (varying intercept + slope) binary logit
via a from-scratch hierarchical Pólya–Gamma Gibbs sampler — the
multilevel extension of the flat-logit PG sampler. Pólya–Gamma augmentation makes
every level conjugate: state effects (αj, βj)
are Gaussian, the group-level regression Δ matrix-normal, and the
random-effect covariance Vb inverse-Wishart. The model
reproduces the rich-state / poor-state voting paradox (Gelman, Shor,
Bafumi & Park 2005): a positive within-state income slope (richer
individuals vote Republican) coexisting with a negative between-state
contextual effect (richer states vote Democratic) — a resolution of Simpson's paradox
that a pooled logit hides. Validated on synthetic data, then fit to the 2004 National
Election Pool exit poll (n=67,381 across 50 states): Δ01 = +0.21
(all 50 state slopes positive), Δ10 = −0.31. Cross-checked three ways —
from-scratch PG-Gibbs, PyMC (LKJ/NUTS), and R's bayesm
(rhierBinLogit) — agreeing to ~0.002 on Δ.
Bayesian Random-Effects Panel Logit
A longitudinal mixed logistic regression — a binary GLMM with a
subject random intercept shared across each subject's repeated measurements,
modelling the within-subject correlation of panel data. The binary, longitudinal twin of the
Epil Poisson GLMM, fit from scratch by
Pólya–Gamma data-augmentation Gibbs (the same trick as the
binary logit) — essential here because the toenail /
onychomycosis trial (294 patients, up to 7 visits) is famous for a very large
random-effect variance (σ ≈ 4) that defeats naive Metropolis and challenges likelihood
quadrature. Onycholysis declines ~33%/month and terbinafine clears it faster (treatment×time
≈ −0.14); because σ is so large, the subject-specific effects are several times larger than a
marginal/GEE model would report. Cross-checked against PyMC (non-centered NUTS) and R
lme4::glmer (with the Lesaffre–Spiessens quadrature caveat).
Bayesian Robust Binary Probit
A Python implementation of robit regression (Liu 2004) via the Geweke (1993) Normal scale-mixture Gibbs sampler. Replaces the standard probit's normal link with a Student-tν CDF: each observation receives a latent precision weight λi ~ Gamma(ν/2, ν/2) that scales its variance — outliers acquire small λi and are downweighted automatically. The 3-block sampler extends Albert & Chib (1993) with a Gamma block for λi; unlike standard probit, the posterior precision B̅-1 = X'ΛX + B0-1 must be recomputed at every iteration. Applied to: Section 1 synthetic t5 data (standard probit attenuated 23%; robit recovers true β), Section 2 sensitivity to ν (LPML peaks at the true value), Section 3 Finney (1947) vasoconstriction (n=39): λi identifies discordant observations without analyst input, including one missed by Liu (2004).
View example →Bayesian Mixture-of-Normals Probit
A Python implementation of the Geweke & Keane (1997) mixture-of-normals probit, where the binary error distribution is the link function — so a flexible mixture εt ~ Σj pj N(αj, hj−1) estimates the link instead of assuming probit's symmetric Φ. A from-scratch 5-block Gibbs sampler (Albert–Chib latent utilities → component indicators → GLS β → Normal-Inverse-Gamma components → Dirichlet weights) handles both a symmetric scale mixture (generalizing the robit) and a full location-scale mixture for asymmetric links. Applied to the Bliss (1935) beetle dose–mortality data — the textbook case where symmetric probit/logit underfit: the mixture discovers the skew on its own, cutting grouped deviance from 10.1 to 3.40, matching the cloglog link (3.45) statisticians historically chose by hand. Cross-checked in R against estimated flexible links (Gosset t-link, Aranda–Ordaz) — symmetric flexibility isn't enough; the asymmetric family drives λ→0 (cloglog), confirming the same story. The binary sibling of the mixture-of-normals regression.
View example →Bayesian Sequential Probit
A Python implementation of two competing ordinal regression models via Albert & Chib (1993, 2001) data-augmentation Gibbs sampling, compared against MLE and bambi/PyMC (HMC). The cumulative model uses a single latent utility and one slope for all transitions (proportional-odds); the sequential model decomposes the ordinal outcome into J−1 independent binary probit stages, each with its own coefficient vector — reducing to separate Albert & Chib (1993) binary probit Gibbs runs on successive at-risk subsets. Applied to Ashford's (1959) pneumoconiosis data (371 coal miners, 8 exposure groups, J=3 severity levels): the Stage-1 exposure slope (onset: 0.804) is 2.6× larger than Stage-2 (progression: 0.313), with P(α11 > α21) = 0.984 — exposing a proportional-odds violation that the cumulative model cannot detect. All three inference methods agree on estimates to within 0.016.
View example →Bayesian Ordered Probit
A Bayesian heteroscedastic ordered probit model estimated via the Albert & Chib (1993)
4-block Gibbs sampler, with a PyMC / NUTS implementation and an R MASS::polr
comparison. Estimates both latent quality (μ) and polarisation (σ) per movie,
revealing cases where the metric mean is misleading. Applied to Amazon Prime movie ratings
(36 movies, 1–5 stars) from Liddell & Kruschke (2018).
Bayesian Multivariate Probit Models
A Python implementation of the Bayesian multivariate probit Gibbs sampler (Edwards & Allenby 2003; Chib & Greenberg 1998), with worked examples using the Scotch Whisky and Ohio Children's Wheeze Study datasets. Includes an R comparison using the bayesm package. Applied in Meseguer (2024) to study correlation patterns among primary and secondary diagnosis codes in the Social Security Disability programs.
View example →Bayesian Multinomial Probit Models
A Python implementation of the Bayesian multinomial probit Gibbs sampler (McCulloch & Rossi 1994; McCulloch, Polson & Rossi 1999), applied to the IRI margarine scanner-panel dataset (516 households, 10 brands). Three specifications — price-only, brand intercepts + price for the top 4 brands, and the full 10-brand model — are estimated and compared. Includes a PyMC implementation via GHK-augmented NUTS and an R comparison using bayesm.
View example →Bayesian Hierarchical Multinomial Logit
A Python implementation of the random-coefficients multinomial logit via a from-scratch
RW-Metropolis-in-Gibbs sampler — per-respondent part-worths
βi ~ N(Δ'zi, Vβ) updated by random-walk
Metropolis (proposal scaled by each respondent's MNL Hessian) inside a conjugate Gibbs
hierarchy (Δ matrix-normal, Vβ inverse-Wishart) — the structure of
bayesm's rhierMnlRwMixture. The notebook works through IIA:
the hierarchical model relaxes it across people (taste heterogeneity), the
complement to multinomial probit relaxing it within a choice (error covariance Σ).
Validated on synthetic data, then fit to the bayesm margarine scanner panel
(516 households, 4,470 purchases, 10 products): strong price aversion (≈ −8.9), with income
lowering and family size raising price sensitivity. Scale-free willingness-to-pay
matches the multinomial probit fit to the cent
despite a different error structure. Cross-checked against R bayesm
(intercepts correlate 0.91).
Bayesian Hierarchical MNL — Mixture Heterogeneity
An extension of the hierarchical multinomial logit that replaces the single-normal
heterogeneity prior with a mixture of normals
(βi | si=k ~ N(μk, Σk)),
replicating bayesm's rhierMnlRwMixture (Rossi, Allenby & McCulloch
Ch. 5) via a from-scratch RW-Metropolis-in-Gibbs sampler. The point is the mixture as a
flexible taste distribution — letting the population of part-worths be
skewed or fat-tailed — not as a segmentation device. Validated by recovering a known
bimodal price sensitivity on synthetic data, then applied to the bayesm
margarine panel: a 2-component mixture is favored (log marginal density up ≈44), revealing
a non-normal, leptokurtic taste distribution — but with heavily overlapping
components (fat tails, not clean segments). The mean part-worths, willingness-to-pay, and
IIA conclusions of the single-normal hier-MNL are all
unchanged — the mixture is a robustness refinement that matters for tail and individual-level
inference. Cross-checked against R bayesm (ncomp=2 wins; same verdict).
Bayesian Tobit — Censored Gaussian Regression
The Tobit model for outcomes observed only within a known region — a sensor
that saturates, a detection limit, or hours worked piling up at 0 — where ordinary regression
is biased. A latent Gaussian y* ~ N(x'β, σ²) is seen only when uncensored, and
a from-scratch data-augmentation Gibbs sampler (Chib 1992) imputes each censored
value from a truncated normal, turning the awkward density/CDF likelihood into an
ordinary conjugate regression — the same augmentation trick as the
ZIP indicator and the
Albert–Chib probit latent. Handles left, right, and
interval censoring in one framework. Applied to the classic Mroz (1987) labor-supply data
(753 women, hours worked left-censored at 0 for the 43% who didn't work): the Tobit recovers
effects that naive OLS attenuates toward zero (a child under 6 ≈ −900 hours), matching the
AER::tobit benchmark. Cross-checked against PyMC (pm.Censored) and R.