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.
What a limited dependent variable is. Ordinary regression assumes the thing being explained can, in principle, take any value on a continuum. A great many outcomes cannot. Some are a yes or no — did this woman take paid work, did this child wheeze. Some are a rank with no arithmetic behind it — a film rated three stars is not the average of two and four. Some are a choice from a menu, where picking one brand means not picking nine others. And some are numbers that hit a wall: hours worked cannot fall below zero, so the 43% of women in the Mroz data who did not work all record exactly 0, and the number tells you they did not work rather than by how much they did not. In each case the outcome is limited — restricted to a set of values that ordinary regression is entitled to run straight past. Fit a straight line to a 0/1 outcome and it will happily predict a probability of 1.4; fit one to hours worked and it will predict −300 hours. The models here exist to respect the restriction rather than ignore it.
Probit or logit — and why both. Two models recur throughout, and a reader is entitled to ask what separates them. Both say the same thing: behind the yes/no you observe sits an unobserved continuous quantity — call it the strength of the inclination — and the answer is yes when it clears a threshold. They differ only in the distribution assumed for the random part of that quantity. Probit assumes it is Normal; logit assumes the logistic, which is very slightly heavier in the tails. The two curves are almost indistinguishable once put on a common scale (a logit coefficient runs about 1.6× the corresponding probit one), and on nearly any dataset they fit equally well and give the same conclusions. What separates them in practice is not fit but convenience. The logit gives interpretable odds ratios and needs no integral to evaluate; the probit inherits the Normal's multivariate machinery, which is what makes correlated equations tractable at all — there is no usable multivariate logistic. That is why the correlated-choice models below (multivariate and multinomial probit) are probits, while the models about a single decision lean on the logit. The collection implements both so that the difference is a demonstrated fact rather than a claim.
Why these models always need a normalisation. A phrase that recurs on every page is that something is fixed, imposed, or identified only up to scale — the error variance set to 1, or the first diagonal of a covariance matrix pinned at 1. This is not a technicality but the central limitation of the data. If the latent inclination were doubled and the threshold doubled with it, every observed yes and no would be exactly the same, so no quantity of data can tell the two apart. The outcome reveals only whether the latent quantity cleared its threshold, never by how much, so its units are unknowable and only ratios of coefficients are estimable. Fixing the scale by hand is what makes the remaining parameters mean something. It is also why coefficients cannot be compared across a probit and a logit without rescaling, why the robust-probit page has to keep three versions of the same coefficient straight, and why willingness-to-pay — a ratio, and so scale-free — is the quantity that can be carried between models.
Odds, since the logit reports them. The odds of an event are its probability divided by the probability of it not happening: a 20% chance is odds of 0.25, or 1 to 4 against. A logit coefficient is a change in log-odds, so exponentiating it gives a multiplier on the odds — an odds ratio. When the binary-logit page reports that each child under six multiplies the odds of labour-force participation by 0.23, it means the odds fall to under a quarter of what they were, all else equal. Odds ratios are convenient because they do not depend on the baseline rate, and misleading for the same reason: a large odds ratio on a rare outcome can still be a very small change in probability.
The one idea underneath all sixteen
Every model in this section rests on the same move. Behind a yes or a no sits an unobserved continuous quantity — an inclination, a utility, a latent score — which is Normal around x′β, and the answer is yes whenever it clears a threshold. That is the left panel, and it explains where a probability comes from: the chance of a yes is just the shaded area, here 0.79.
The green strip is the part that makes these models computable. Reading the picture forwards gives a likelihood that is awkward to sample from directly. Reading it backwards — having observed a yes, the latent value must have been positive, so redraw it from just the shaded part — turns the problem into an ordinary regression on quantities you have just invented. That single step is Albert–Chib data augmentation, and every sampler below is a variation on it: the ordered probit adds more thresholds, the multivariate probit draws a whole vector of correlated latent values at once, and the Tobit augments the censored value itself rather than a utility.
The right panel settles the question of why two models exist for one job. Probit assumes the random part of that latent quantity is Normal, logit assumes the logistic — and once the logistic is rescaled by the factor of about 1.6 that separates their spreads, the two curves never differ by more than 0.017 in probability anywhere. They are, for practical purposes, the same model. What separates them is convenience: the logit reports odds ratios and needs no integral, while the probit inherits the multivariate Normal, which is the only reason correlated-choice models are tractable at all.
How the sixteen examples relate
Four kinds of restricted outcome, each with one trunk model that establishes the sampler and a set of extensions that relax one assumption apiece. Read a trunk first; the branches below it assume it. Three entries are cross-references — their own homes are the count and regression sections — but they belong to the same idea and are listed where they fit.
Yes / no
Did she take paid work? Did this child wheeze?
Binary ProbitAlbert–Chib augmentation — the sampler everything else here is built on Robust probit (robit)a t-link, so mislabelled cases stop dictating the fit Mixture-of-normals probitstops assuming the link and estimates its shape instead Binary LogitPólya–Gamma augmentation — the logit’s answer to Albert–Chib Hierarchical binary logiteach state gets its own slope — the rich-state/poor-state paradox Random-effects panel logitrepeated visits per patient, so the responses are correlated Hierarchical Poissonthe same random-intercept idea, counted outcomes instead of binaryRanked
Normal, mild or severe? Three stars or four?
Sequential Probitone latent score and a set of cut-points — and a test of whether a single slope can govern every step Ordered probit, heteroscedasticadds a per-group spread, separating a divisive film from a mediocre oneChoice from a menu
Which of ten margarines? Which whiskies did this household buy?
Multivariate Probitseveral yes/no answers at once, with the correlations between them estimated Multinomial Probitexactly one option chosen; correlated errors relax IIA within a choice Hierarchical MNLrelaxes IIA the other way — across people, by letting tastes vary Mixture heterogeneityasks whether that population of tastes is one bell curve or severalHitting a wall
Hours worked that pile up at zero; a day’s fishing that caught nothing.
Tobitcensoring — the value exists but was not recorded past a limit Zero-inflated countstwo kinds of zero: never at risk, or at risk and unlucky Mixture-of-normals regressionthe continuous sibling of the mixture-link probit aboveRead across the row rather than down it and a second pattern appears: the same three moves — make it robust, make it hierarchical, make the distribution flexible — are applied to each trunk in turn. Almost nothing here is a new idea; it is one augmentation trick and three ways of loosening it.
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.