Multilevel Ordinal Regression

Definition

Multilevel ordinal regression extends the proportional odds (or partial proportional odds) model to clustered or longitudinal data by adding subject-specific random effects $\boldsymbol{\theta}_j \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$. For observation $i$ nested in level-2 unit $j$ and cumulative category boundary $c$, the model is $\log[P_{ijc}/(1-P_{ijc})] = \gamma_c - [x_{ij}'\beta + z_{ij}'T\theta_j]$, where $TT'=\Omega$ is the Cholesky-parameterized random-effect covariance. Because there is a single $\beta$ vector (one latent continuum, cut by thresholds), the model inherits the "threshold concept" of the cumulative logit: all category transitions are driven by the same covariate effects, differing only in their cut-points $\gamma_c$. Estimation proceeds by maximum marginal likelihood (MML) via Gauss-Hermite quadrature (GHQ) and Fisher scoring.

Key Ideas

• Proportional odds: One regression coefficient vector $\beta$ governs all $C-1$ cumulative log-odds; the log-odds ratio for any covariate is assumed constant across thresholds (McCullagh 1980). A single intraclass correlation (ICC) $= \hat\sigma^2/(\hat\sigma^2 + \pi^2/3)$ summarises between-subject clustering for the ordinal outcome.
• Partial proportional odds (Hedeker-Mermelstein 1998): A subset $x^*_{ij}$ of $h$ covariates is allowed category-varying coefficients $\beta_c$. Relaxing proportionality for relevant predictors (e.g., treatment whose effect changes across severity levels) improves fit, but regression lines for continuous $x^*$ must eventually cross, potentially producing inadmissible (negative) category probabilities at extreme covariate values.
• Complementary log-log and discrete-time survival: Replacing the logit link with $\log[-\log(1-P_{ijc})]$ maps the model onto the grouped-time proportional hazard (PH) model (Prentice-Gloeckler 1978; Singer-Willett 1993). Coefficients are invariant to interval length under the PH assumption. The ordinal formulation (eq. 6.6 in Hedeker 2007) and the pooled-dichotomous formulation (eq. 6.7) are equivalent when covariates are time-invariant; the pooled form handles time-varying covariates without reformulating the model.
• IRT connection: Under homogeneous item discrimination, the multilevel logistic / ordinal model is the Rasch model (random subjects, fixed items as fixed effects). With item-specific discrimination parameters $a_m$ in the heterogeneous-$T$ formulation ($T=[a_1\cdots a_m]'$), the model is equivalent to Samejima's graded-response item response theory (IRT) model for ordinal items, enabling a unified treatment of multilevel and IRT estimation via quadrature.
• Heterogeneous variance: Including $z_{ij}$ as a $g\times 1$ group-indicator vector yields group-specific random-effect SDs. For binary outcomes, this recovers polychoric and tetrachoric correlations in twin or family designs.
• EAP estimates: After convergence, empirical Bayes (EAP) posterior means $\hat\theta_j = h(Y_j)^{-1}\int\theta_j\ell(\cdot)g(\theta)d\theta$ (eq. 6.19 of Hedeker 2007) provide subject-specific random-effect estimates by a second quadrature pass.

How It Works

1. Conditional likelihood: Given $\theta_j$, the observations within unit $j$ are conditionally independent ordinal draws; the product over time points gives $\ell(\mathbf{y}_j | \theta_j) = \prod_i \prod_c [P_{ijc} - P_{ij,c-1}]^{y_{ijc}}$.
2. Marginal likelihood: $h(\mathbf{y}_j) = \int \ell(\mathbf{y}_j | \theta)g(\theta)d\theta$; approximated by Gauss-Hermite quadrature with $Q^r$ points ($r$ = number of random effects).
3. Fisher scoring: Gradient and observed information are summed over units and quadrature nodes; iterative updates of $(\beta, \boldsymbol{\gamma}, T)$ converge to the MML estimate.
4. Proportionality testing: likelihood-ratio (LR) $\chi^2$ with $h(C-2)$ degrees of freedom tests the null that all $h$ covariates satisfy proportional odds; partial non-proportionality models are nested within the full non-proportional model.
5. Computational alternatives: marginal/penalized quasi-likelihood (MQL/PQL) linearisation can bias variance estimates; Laplace approximation (as in HLM software) is fast and often sufficiently accurate; adaptive quadrature (3 pts/dim) outperforms ordinary GHQ for $r>2$; Markov Chain Monte Carlo (MCMC) is preferred for large $r$.

Why It Matters

• Multilevel proportional odds handles the overwhelming majority of clinical, educational, and social-science ordinal outcomes (Likert scales, symptom severity, disease stage) while correctly accounting for within-cluster correlation — something the marginal generalized estimating equations (GEE) approach cannot do when subject-specific inference is needed.
• The partial proportional odds extension provides a parsimonious compromise between the rigid proportional odds assumption and the fully non-proportional model, maintaining interpretability while accommodating covariate-specific threshold effects.
• The discrete-time survival interpretation unifies event-history and ordinal-data frameworks: the same software (MIXOR) handles both, and the cloglog link makes hazard-ratio interpretation natural.
• The IRT bridge enables multilevel software to estimate psychometric models without dedicated IRT packages, and vice versa: IRT software with random group effects can handle standard multilevel ordinal problems.

Open Questions

• When the number of categories $C$ is large, the partial proportional odds model with $h(C-2)$ extra parameters may overfit; penalised or Bayesian approaches to selecting which covariates are non-proportional are an active area.
• Adaptive quadrature reduces node counts but the $Q^r$ grid still scales exponentially with $r$; for three-level models or models with many random slopes, Laplace or MCMC remains the practical choice.
• Robustness to random-effect (RE) misspecification: normal random effects may be inappropriate when the population distribution of individual trajectories is multimodal or heavy-tailed; nonparametric mixture alternatives exist but are rarely used in practice for ordinal outcomes.

Related

• Ordinal Regression
• Mixed-Effects Multinomial Regression
• Mixed-Effects Regression
• Random Effects Probit
• Generalized Linear Mixed Model
• Sequential Ordinal Model
• Donald Hedeker
• Robert D. Gibbons
• Hedeker (2007) Multilevel Models for Ordinal and Nominal Variables