Fitting Ordinal Mixed Models with ordinal

Why Model Pain Ratings as Ordered Categories?

A pain rating such as none, mild, moderate, or severe contains more information than a binary outcome, because the categories have a natural order. But the gaps between adjacent categories are not necessarily equal: moving from “none” to “mild” need not represent the same change in pain as moving from “moderate” to “severe.” A cumulative-link mixed model uses the ordering without treating the category labels as equally spaced numbers, and a random effect can account for repeated ratings from the same patient.

This tutorial fits that model with ordinal::clmm. The examples use hypothetical pain data. The model’s interpretation depends on how the response is ordered and on the software’s cumulative-logit sign convention, so those are stated explicitly before interpreting output.

Definition: Let $Y_{ij}$ be the ordered response for observation $i$ in subject $j$, taking one of $K$ ordered categories. A cumulative-link mixed model with a subject random intercept can be written $$ \operatorname{logit}\{P(Y_{ij}\leq k\mid u_j)\} =\theta_k-\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}-u_j,\qquad k=1,\ldots,K-1. $$ Here $\theta_1<\cdots<\theta_{K-1}$ are thresholds, $\mathbf{x}_{ij}$ is the predictor vector, $\boldsymbol{\beta}$ is the vector of fixed-effect coefficients, and $u_j\sim N(0,\sigma_u^2)$ is a subject-specific random intercept. Conditional on the random effects and predictors, observations are modeled as independent. The model also assumes that each predictor has the same coefficient across all cumulative splits; this is the proportional-odds assumption.

With the logit link, the model expresses cumulative probabilities $P(Y\leq k)$, such as the chance of rating “moderate” or lower. Because the formula uses $\theta_k-\mathbf{x}^{\mathsf T}\boldsymbol{\beta}-u$, a positive coefficient shifts probability toward higher response categories. In particular, for a one-unit increase in a predictor, the conditional odds of being above any threshold rather than at or below it are multiplied by $\exp(\beta)$, holding other predictors and the random effect fixed.

Set Up the Response and Fit the Model

Suppose a hypothetical study records pain at baseline, week 2, and week 4 for 60 patients assigned to control or treatment. Each record contains a patient identifier, a treatment group, a visit, and a pain rating. The four response levels are none, mild, moderate, and severe. The response must be an ordered factor, with levels in the intended order; an ordinary factor or a numeric code may not communicate the intended ordinal structure to the model.

An illustrative data preparation and model call are:

R: dat$pain <- ordered(dat$pain, levels = c("none", "mild", "moderate", "severe"))
R: dat$treatment <- relevel(factor(dat$treatment), ref = "control")
R: fit <- ordinal::clmm(pain ~ treatment + visit + (1 | patient), data = dat, link = "logit", Hess = TRUE)

The random intercept allows patients to differ in their underlying tendency to report higher pain. The fixed effects describe the treatment and visit associations conditional on that patient effect. This specification assumes the same treatment and visit slopes across patients; it does not include patient-specific slopes. Check that each row is one rating, that patient identifiers match the repeated-measures structure, and that no response level has been accidentally reordered or dropped.

Worked Example: Reading Thresholds and a Treatment Coefficient

Suppose an illustrative clmm summary reports thresholds $\hat{\theta}_1=-0.50$, $\hat{\theta}_2=1.00$, and $\hat{\theta}_3=2.30$, and a treatment coefficient $\hat{\beta}_{\text{treatment}}=-0.70$ with standard error $0.24$. Treatment is coded relative to control. These numbers are hypothetical software output, not results from a real study.

The thresholds correspond to the boundaries none|mild, mild|moderate, and moderate|severe. They are cut points on the model’s latent cumulative-logit scale. They are not treatment effects and are not the log odds of each individual category.

The treatment coefficient is negative. For treatment versus control, the conditional cumulative odds ratio for being above any pain threshold is $\exp(-0.70)=0.4966$. Thus, for patients with the same visit and random intercept, the fitted odds of being above a given threshold are about $0.497$ times as large in the treatment group as in control. This is an association in the model, not by itself evidence of a causal treatment effect.

A Wald calculation gives $z=-0.70/0.24=-2.917$, with a two-sided normal-approximation $p$-value of about $0.0035$. A 95% Wald interval for the coefficient is $-0.70\pm1.96(0.24)=(-1.1704,-0.2296)$. Exponentiating the endpoints gives an approximate interval $(0.310,0.795)$ for the conditional cumulative odds ratio. Under the model and its assumptions, this interval is compatible with lower cumulative odds of higher pain in the treatment group; it does not assign a 95% probability to the fixed parameter.

From Thresholds to Category Probabilities

Thresholds and fixed effects can be combined to obtain fitted cumulative probabilities for a specified predictor pattern and random effect. Category probabilities are differences between adjacent cumulative probabilities. For the first category, the probability is the first cumulative probability; for the last category, it is one minus the final cumulative probability.

Key result: Define $\eta_{ij}=\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+u_j$ and $F(a)=1/(1+\exp(-a))$. Then $P(Y_{ij}\leq k\mid u_j)=F(\theta_k-\eta_{ij})$. The category probabilities are $P(Y=1)=F(\theta_1-\eta)$, $P(Y=r)=F(\theta_r-\eta)-F(\theta_{r-1}-\eta)$ for $2\leq r\leq K-1$, and $P(Y=K)=1-F(\theta_{K-1}-\eta)$.

This follows directly from the cumulative probabilities. The event $Y=r$, for an interior category, is the event $Y\leq r$ minus the event $Y\leq r-1$. Therefore its probability is $P(Y\leq r)-P(Y\leq r-1)$. The first category has no preceding category, and the final category is the complement of $Y\leq K-1$. Since the thresholds increase and $F$ is increasing, the resulting category probabilities are nonnegative and sum to one.

Worked Example: Calculate Pain-Category Probabilities

Use the illustrative thresholds $-0.50, 1.00, 2.30$ from above. For a control patient at the reference visit with random intercept $u=0$, take $\eta=0$. The cumulative probabilities are $F(-0.50)=0.3775$, $F(1.00)=0.7311$, and $F(2.30)=0.9089$.

The four category probabilities, calculated from the unrounded cumulative probabilities, are:

None: $F(-0.50)=0.3775$.
Mild: $F(1.00)-F(-0.50)=0.3535$.
Moderate: $F(2.30)-F(1.00)=0.1778$.
Severe: $1-F(2.30)=0.0911$.

The rounded values sum to $0.9999$, which is 1 up to rounding. For a treatment patient with the same visit and $u=0$, the treatment coefficient gives $\eta=-0.70$. The cumulative probabilities are $F(0.20)=0.5498$, $F(1.70)=0.8455$, and $F(3.00)=0.9526$. Thus the category probabilities are approximately $0.5498$, $0.2957$, $0.1070$, and $0.0474$. The fitted probability of severe pain is lower for this pattern, consistent with the negative treatment coefficient. These are conditional predictions for $u=0$, not probabilities averaged over the distribution of patient effects.

What the Random Intercept Changes

The patient random intercept shifts all cumulative logits for that patient together. It does not create a separate treatment coefficient for every patient. A positive $u_j$ raises $\eta_{ij}$, reducing each $P(Y\leq k)$ and shifting probability toward higher categories; a negative value shifts probability toward lower categories.

Worked Example: A Patient With a Positive Random Intercept

Use the same thresholds for a control patient at the reference visit, but suppose the fitted patient-specific intercept is $u=0.60$. Then $\eta=0.60$. The cumulative probabilities are $F(-0.50-0.60)=F(-1.10)=0.2497$, $F(1.00-0.60)=F(0.40)=0.5987$, and $F(2.30-0.60)=F(1.70)=0.8455$.

The category probabilities are therefore none $=0.2497$, mild $=0.5987-0.2497=0.3489$, moderate $=0.8455-0.5987=0.2468$, and severe $=1-0.8455=0.1545$, using unrounded cumulative values for the differences. Compared with $u=0$, this patient has less fitted probability in the lowest category and more in the higher categories. The calculation illustrates a conditional, patient-specific prediction; it does not mean the model estimates a distinct fixed effect for this patient.

Assumptions, Diagnostics, and Failure Modes

Interpretation requires the response categories to be genuinely ordered and the ordering to be correctly encoded. The cumulative-logit model assumes a common predictor effect across the $K-1$ cumulative splits. For example, the treatment coefficient is constrained to be the same for “above none” versus “none or lower” and for “severe” versus “moderate or lower.” This is the proportional-odds assumption. The normal random-intercept distribution and the specified dependence structure are also modeling assumptions, not facts guaranteed by using mixed-model software.

Review the model fit and data before relying on estimates. Check convergence messages and whether the Hessian is usable for standard errors. Look for very sparse or empty response categories, especially within treatment or visit groups; sparse data can produce unstable thresholds or large standard errors. Plot observed ratings by visit and treatment, and compare them with fitted probabilities. Also inspect patient-level patterns: a random-intercept model may be inadequate if patients show substantially different changes over time. Such patterns motivate reconsidering the random-effects structure rather than interpreting a single average slope as universal.

A useful diagnostic question is whether the treatment or visit association plausibly remains constant across cumulative thresholds. Large discrepancies in observed-versus-fitted cumulative patterns can indicate that the proportional-odds constraint is too restrictive. This check is substantive as well as graphical: a treatment could plausibly affect the chance of any pain differently from its effect on the chance of severe pain. The standard errors and p-values from the fitted model do not test every assumption automatically.

Verify the response and rows.
Confirm the category order, one record per rating, and correct subject identifiers. Check missingness and category counts by predictor groups.

Fit the stated conditional model.
Use an ordered response, fixed effects that match the study question, and random effects that match the repeated-measures design.

Review estimates and fit.
Check convergence, threshold ordering, uncertainty, observed category patterns, and whether common cumulative slopes are plausible.

Report effects in context.
State the reference group, conditional cumulative odds ratio, uncertainty, and whether any probabilities are conditional on a random effect or averaged over it.

Common Mistakes and Alternatives

Reading thresholds as predictor effects. Thresholds locate the category boundaries on the latent scale; the predictor coefficients describe shifts relative to those boundaries.
Reversing the sign interpretation. Under the convention $\operatorname{logit}\{P(Y\leq k)\}=\theta_k-\eta$, a positive coefficient shifts toward higher categories. Always check the response order and model convention.
Using integer scores as if the gaps were equal. Regressing scores 0, 1, 2, and 3 treats their spacing as meaningful. A cumulative-link model uses order without imposing that equal-spacing assumption.
Calling a conditional prediction a population-average probability. A prediction at $u=0$ describes a subject with random intercept zero. Population-average probabilities require accounting for variation in $u$.

A cumulative-link mixed model is appropriate when the response is ordered and repeated or clustered observations require a random effect. A binary mixed model discards distinctions among multiple ordered categories. A multinomial mixed model does not use the ordering in the same way and estimates a more flexible set of category comparisons, often with more parameters. Treating category scores as continuous is simpler but adds an equal-spacing interpretation. None of these alternatives is universally preferable: the response scale and scientific question determine which assumptions are defensible.

The model also has limits. The proportional-odds restriction may be unsuitable, estimates can be unstable with sparse categories or few clusters, and random-effects predictions are uncertain, especially for subjects with little data. A fitted model does not resolve design limitations or establish causation. Report the response ordering, link, fixed and random effects, reference levels, and the conditional nature of the reported odds ratios or predictions.

Check Your Understanding

Answer each question, then compare with the brief solutions.

For the stated cumulative-logit convention, what does a positive treatment coefficient imply about higher pain categories?
If the coefficient is $-0.40$, what is the conditional cumulative odds ratio for being above a pain threshold?
Why is a threshold not the log odds of a particular response category?
What does a prediction at random intercept $u=0$ describe?
Name one diagnostic concern if a response category is empty in a treatment-by-visit group.

Solutions. (1) It shifts fitted probability toward higher categories, increasing the odds of being above each threshold. (2) $\exp(-0.40)=0.6703$, so the conditional odds are multiplied by about $0.670$. (3) A threshold defines a boundary for a cumulative probability $P(Y\leq k)$; an individual-category probability is obtained from differences between cumulative probabilities. (4) It describes the conditional prediction for a subject whose random intercept is zero, not an average over all subjects. (5) Sparse or empty cells can make thresholds and effects unstable and can undermine reliable model fitting or inference.