Ordinal Regression

Definition

Ordinal regression models the probability of a discrete ordered outcome $Y_i \in \{1,\ldots,J+1\}$ as a function of covariates $\mathbf{x}_i$ using either a single latent continuous variable crossing ordered thresholds (cumulative model) or a sequence of binary decisions at each level (sequential model). The three principal non-nested specifications — cumulative probit, sequential probit, and two-step compound — differ in how they decompose $\Pr(Y_i = j)$ , and cannot be compared by classical likelihood ratio tests; Bayesian model comparison via Bayes factors is required.

Key Ideas

Non-nestedness: the three main ordinal models (cumulative, sequential, two-step) are not special cases of one another; Bayes factors computed from marginal likelihoods provide the only principled comparison tool.
Cut-point ordering constraint: the cumulative model's thresholds $\gamma_1 \leq \cdots \leq \gamma_{J-1}$ are ordered, complicating Markov Chain Monte Carlo (MCMC); the log-spacing reparameterization of Albert-Chib (1997b) lifts them to an unrestricted real vector.
Latent variable structure: all three models admit data augmentation — introducing continuous latent variables $z_i$ that reduce each model to a Normal regression conditional on the augmented data.
Prior specification for Bayes factors: valid Bayes factors require proper priors; Albert-Chib (1997b) develops three elicitation strategies (training sample, imaginary prior sample, Dirichlet-on-multinomial) specifically for this purpose.

How It Works

Three Model Families

Cumulative probit (McCullagh 1980): $\Pr(Y_i \leq j \mid \boldsymbol{\beta}, \boldsymbol{\gamma}) = \Phi(\gamma_j - \mathbf{x}_i'\boldsymbol{\beta}), \quad \gamma_1 \leq \cdots \leq \gamma_{J-1}$ A single latent variable $z_i \sim \mathcal{N}(\mathbf{x}_i'\boldsymbol{\beta}, 1)$ determines category via $Y_i = j \iff \gamma_{j-1} < z_i \leq \gamma_j$ . The ordering constraint on $\boldsymbol{\gamma}$ is the central algorithmic difficulty.

Sequential probit (Tutz 1990, 1991): $\Pr(Y_i = j \mid Y_i \geq j, \boldsymbol{\beta}, \boldsymbol{\gamma}) = \Phi(\gamma_j - \mathbf{x}_i'\boldsymbol{\beta})$ Each observation $Y_i = j$ requires $j_i = \min\{j+1, J\}$ category-specific independent latent variables $z_{ij} \sim \mathcal{N}(\mathbf{x}_i'\boldsymbol{\beta}, 1)$ . The thresholds $\boldsymbol{\gamma}$ are unrestricted. Equivalent to a discrete-time hazard / continuation-ratio model. See Sequential Ordinal Model for a detailed treatment.

Two-step / compound model (Tutz 1989): Partition the $J+1$ categories into $s$ ordered subsets; fit a cumulative probit at the subset level and a separate cumulative probit within each subset. The joint posterior factors into $s+1$ independent components, so inference reduces to $s+1$ parallel applications of the cumulative algorithm.

Cut-Point Reparameterization (Algorithm 1 of Albert-Chib 1997b)

To handle the ordered-threshold constraint in the cumulative and two-step models, transform $\boldsymbol{\gamma}$ to an unrestricted vector $\boldsymbol{\alpha}$ : $\alpha_1 = \log \gamma_1, \qquad \alpha_j = \log(\gamma_j - \gamma_{j-1}), \quad j \geq 2$ $\text{Inverse: } \gamma_j = \sum_{k=1}^{j} \exp(\alpha_k)$ A multivariate- $t$ Metropolis-Hastings proposal $f_T(\boldsymbol{\alpha} \mid \hat{\boldsymbol{\alpha}}, \tau^2 V, \nu)$ is centred at the mode $\hat{\boldsymbol{\alpha}}$ with curvature $V$ from Newton-Raphson on $\ln \pi(\boldsymbol{\alpha} \mid \mathbf{y}, \boldsymbol{\beta})$ . The Jacobian $|\partial\boldsymbol{\gamma}/\partial\boldsymbol{\alpha}| = \prod_j \exp(\alpha_j)$ adjusts the acceptance ratio.

Marginal Likelihood Computation (Chib 1995)

Using the Chib (1995) identity evaluated at a posterior mode $\psi^* = (\boldsymbol{\beta}^*, \boldsymbol{\alpha}^*)$ : $\ln m(\mathbf{y}) = \ln f(\mathbf{y} \mid \psi^*) + \ln \pi(\psi^*) - \ln \pi(\psi^* \mid \mathbf{y})$ The posterior ordinate is decomposed as $\pi(\boldsymbol{\beta}^*, \boldsymbol{\alpha}^* \mid \mathbf{y}) = \pi(\boldsymbol{\beta}^* \mid \mathbf{y}) \cdot \pi(\boldsymbol{\alpha}^* \mid \mathbf{y}, \boldsymbol{\beta}^*)$ . The $\boldsymbol{\beta}^*$ ordinate is estimated by Rao-Blackwellisation over MCMC latent-variable draws; the $\boldsymbol{\alpha}^*$ ordinate by Gaussian kernel smoothing (Scott 1992 bandwidth) from a reduced run fixing $\boldsymbol{\beta} = \boldsymbol{\beta}^*$ .

Prior Elicitation for Bayes Factors

Training sample: split data, fit flat-prior model to the training portion, use posterior moments as hyperparameters for the main analysis.
Imaginary prior sample: construct an artificial dataset implying equal ordinal frequencies and no covariate associations; use its implied posterior as the prior.
Dirichlet-on-multinomial: place Dirichlet priors on response probability vectors $\pi_k$ at chosen covariate scenarios; transform to a prior on $(\boldsymbol{\beta}, \boldsymbol{\gamma})$ via the Jacobian.

Why It Matters

Ordinal outcomes (Likert scales, grades, health ratings, visit counts) are ubiquitous; choosing the wrong structural model (cumulative vs. sequential) leads to systematically biased inference.
The log-spacing reparameterization converts the constrained cut-point sampling problem into an unconstrained one, achieving fast mixing without the specialised algorithms of Cowles (1996) or Chen-Nandram (1996).
The Bayes factor framework allows formal, probabilistic model selection among non-nested alternatives using a single unified MCMC run per model.
Bayesian latent residuals $r_i = z_i - \mathbf{x}_i'\boldsymbol{\beta}$ (prior $r_i \sim \mathcal{N}(0,1)$ ) provide observation-level diagnostics: tail probabilities $\Pr(|r_i| > 2 \mid \mathbf{y}) > 0.05$ flag poorly fitted units.

Failure Modes of Metric Models for Ordinal Data

Liddell-Kruschke (2019) prove that applying metric models (t-test, analysis of variance (ANOVA), regression) to ordinal data can produce all three categories of inferential error. The root cause is geometric: the ordinal mean $\bar{Y}(\mu,\sigma) = \sum_{k} k \bigl[\Phi\bigl((\theta_k - \mu)/\sigma\bigr) - \Phi\bigl((\theta_{k-1}-\mu)/\sigma\bigr)\bigr]$ is a sigmoid function of the latent mean $\mu$ , with shape controlled by the latent standard deviation (SD) $\sigma$ . Two groups with different $\sigma$ trace two distinct sigmoid curves over $\mu$ . Any vertical slice (same $\mu$ , different $\sigma$ ) or horizontal slice (same $\bar{Y}$ , different $\mu$ ) through those curves generates an error.

False alarm (Type I): Groups with equal latent means $\mu_1 = \mu_2$ but unequal latent SDs $\sigma_1 \neq \sigma_2$ produce unequal ordinal means $\bar{Y}_1 \neq \bar{Y}_2$ . A metric model reports a spurious group effect.

Miss (Type II): Groups with unequal latent means $\mu_1 \neq \mu_2$ can have equal ordinal means when the sigmoid curves intersect at those $\mu$ values. A metric model reports no effect.

Inversion: When the two sigmoid curves cross, the ordinal difference $\bar{Y}_1 - \bar{Y}_2$ can have the opposite sign from the latent difference $\mu_1 - \mu_2$ . Extends to factorial interactions (a genuine cross-over becomes non-crossing and vice versa) and trend analyses.

Two things that do not help: (a) averaging multiple correlated Likert items preserves the sigmoidal distortion — all three failure modes persist; (b) equal observed SDs across groups do not imply equal latent SDs.

Empirical prevalence: An audit of the Journal of Personality and Social Psychology (JPSP), Psychological Science (PS), and Journal of Experimental Psychology: General (JEP:G) (2016) found that 100% of 68 articles using Likert data analysed it with metric models.

Open Questions

Sensitivity of Bayes factors to the choice among the three prior elicitation methods; no single method dominates in all settings.
Computational cost of kernel smoothing for the $\boldsymbol{\alpha}^*$ ordinate scales poorly when the number of cut-points is large.
Extension to random-effects or longitudinal ordinal data with subject-specific latent variables — addressed by the maximum marginal likelihood / Gauss-Hermite approach in Multilevel Ordinal Regression (Hedeker 2007), which operates in a frequentist rather than Bayesian framework.
The Liddell-Kruschke (2019) argument conflates two separable issues: (1) using an ordinal model vs. a metric model, and (2) Bayesian vs. frequentist inference. A frequentist proportional-odds logistic regression would solve (1) without requiring MCMC. The relative severity of violations in practice (e.g., small vs. large latent SD heterogeneity) remains an open empirical question.