A Longitudinal Ordinal Pain Study

Why Pain Ratings Need an Ordinal Longitudinal Model

A pain score recorded as “none,” “mild,” “moderate,” or “severe” has a natural order, but the gaps between adjacent categories are not known to be equal. If the same patient reports pain at several follow-up visits, those records are also related: a patient’s ratings may tend to be consistently higher or lower than another patient’s ratings. A longitudinal ordinal mixed model addresses both features by modeling ordered categories and accounting for repeated observations within patients.

The case below is hypothetical. Patients are assigned to an active treatment or usual care, and pain is assessed at weeks 2, 4, and 8 after treatment begins. Baseline pain is recorded before treatment. The aim is to estimate the treatment association with pain category over follow-up, adjusting for visit and baseline pain, while allowing each patient to have a different underlying tendency toward higher pain.

Model: Let $Y_{ij}$ be the ordered pain category for patient $i$ at follow-up visit $j$, with categories none, mild, moderate, and severe. A cumulative-logit random-intercept model is $$ \operatorname{logit}\{P(Y_{ij}\leq k\mid b_i,\mathbf{x}_{ij})\} =\theta_k-\eta_{ij}, \qquad k=1,2,3, $$ where $\eta_{ij}=\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+b_i$, and $b_i\sim N(0,\sigma_b^2)$. The $\theta_k$ are ordered thresholds, $\theta_1<\theta_2<\theta_3$. With this sign convention, a positive coefficient shifts probability toward higher (worse) pain categories; a negative coefficient shifts it toward lower categories.

For the proposed analysis, $\mathbf{x}_{ij}$ includes treatment, visit as a categorical predictor, and baseline pain as a categorical predictor. Usual care is the treatment reference category, and week 2 is the visit reference category. The treatment coefficient is assumed common across the three follow-up visits: the model does not include a treatment-by-visit interaction. Thus its interpretation is an adjusted treatment contrast that is constant on the cumulative-logit scale over the modeled visits.

Prepare the Repeated Ratings

Use one row per patient per observed visit. This is long format: patient identifiers repeat across rows, while each row contains one outcome rating and the covariates for that observation. Store the outcome as an ordered factor with its levels explicitly set in the clinical order. Do not alphabetize the labels or treat category numbers as measurements with equal spacing.

The following small excerpt is hypothetical. Each patient contributes a baseline record and three planned follow-up records in the full data set; only baseline and week 2 are shown here to illustrate the layout.

Patient	Treatment	Visit	Baseline pain	Current pain
P01	Usual care	Baseline	Moderate	Moderate
P01	Usual care	Week 2	Moderate	Mild
P02	Active	Baseline	Severe	Severe
P02	Active	Week 2	Severe	Moderate

Worked Example: Checking the Long-Format Records

Suppose the full hypothetical study has 80 patients, each scheduled for one baseline and three follow-up visits. If all visits are observed, there should be $80\times4=320$ rows, with four rows per patient. The primary follow-up model uses the 240 post-baseline rows; baseline pain is carried as a covariate on each of those rows.

Before fitting, check that each patient has no more than one record at each visit, that the treatment assignment is constant within patient, and that all observed pain values belong to the specified ordered levels. For instance, two P01 records both labeled week 2 would be a duplicate to investigate, not two independent measurements to retain automatically.

In R, an ordered outcome can be created with levels none, mild, moderate, severe, and the follow-up data can be fit using ordinal::clmm. A suitable call is clmm(pain ~ treatment + visit + baseline_pain + (1 | patient), data = followup, link = "logit", Hess = TRUE), after setting the reference levels of treatment and visit. The random intercept groups the repeated ratings by patient.

The model assumes separate patients are independent, conditional on the predictors and patient random effects. Within a patient, the shared random intercept induces dependence across visits. It also assumes normally distributed random intercepts independent of the included predictors, a logistic cumulative link, and proportional odds: each predictor has the same coefficient at all three cumulative cutpoints. If missing follow-up ratings are present, likelihood-based inference additionally relies on missingness being ignorable given the observed data and model, rather than depending on unobserved pain after conditioning on that information.

Interpret the Treatment Coefficient

In this model, the treatment coefficient is not a change in mean pain score. It is a shift in the cumulative log odds of being above versus at or below each pain threshold. For a given threshold $k$, the model implies

\frac{\operatorname{odds}(Y_{ij}>k\mid \text{active},b_i,\text{covariates})} {\operatorname{odds}(Y_{ij}>k\mid \text{usual care},b_i,\text{covariates})} =\exp(\beta_{\text{treatment}}).

This is a patient-conditional comparison: it holds the random intercept and other predictors fixed. The proportional-odds assumption makes this same odds ratio apply at the none-versus-higher, none-or-mild-versus-higher, and none-or-mild-or-moderate-versus-severe divisions. It does not say that every category’s probability changes by the same amount.

Worked Example: Translating a Treatment Estimate Into an Odds Ratio

Suppose a hypothetical fitted model reports a treatment coefficient of $\hat{\beta}_{\text{treatment}}=-0.70$ with $SE=0.28$, where active treatment is compared with usual care. The odds ratio for being above any given pain threshold is

\exp(-0.70)=0.4966.

Thus, conditional on the same patient random intercept, visit, and baseline pain, the estimated odds of being in a higher rather than a lower-or-equal pain category under active treatment are about 0.497 times the corresponding odds under usual care. This is a 50.3% lower estimated cumulative odds, not a 50.3% reduction in pain ratings or in the probability of severe pain.

A Wald 95% interval on the log-odds scale is $-0.70\pm1.96(0.28)$, or $(-1.2488,-0.1512)$. Exponentiating both endpoints gives an odds-ratio interval of approximately $(0.287,0.860)$. Under the model and the Wald approximation, this interval is compatible with cumulative odds between 0.287 and 0.860 times those under usual care. The Wald statistic is $z=-0.70/0.28=-2.50$, giving a two-sided $p$-value of about 0.0124. This is evidence against a zero treatment coefficient under the model; it is not the probability that the null hypothesis is true.

If treatment was randomized and the model and missing-data assumptions are credible, the contrast can support a treatment-effect interpretation. The estimate itself is still conditional on the model specification, and it should not be described as proving that treatment caused an individual patient’s pain to decrease.

From Cumulative Probabilities to Category Probabilities

The fitted model returns cumulative probabilities at the thresholds. Individual category probabilities are obtained by differencing adjacent cumulative probabilities. This calculation makes the clinical category context concrete, but it should use unrounded cumulative probabilities until the final reporting step.

Worked Example: Calculating Four Pain-Category Probabilities

Consider illustrative thresholds $\theta_1=-1.2$, $\theta_2=0.2$, and $\theta_3=1.5$, and set the linear predictor $\eta$ to zero for a reference covariate profile with random intercept $b_i=0$. These are hypothetical values chosen to demonstrate the calculation, not results from the study above. With the logistic function $\operatorname{logistic}(u)=1/(1+\exp(-u))$, the cumulative probabilities are

P(Y\leq\text{none})=\operatorname{logistic}(-1.2)=0.231475, \quad P(Y\leq\text{mild})=\operatorname{logistic}(0.2)=0.549834, \quad P(Y\leq\text{moderate})=\operatorname{logistic}(1.5)=0.817574.

Therefore, the category probabilities are:

None: $0.231475$, reported as 0.2315.
Mild: $0.549834-0.231475=0.318359$, reported as 0.3184.
Moderate: $0.817574-0.549834=0.267741$, reported as 0.2677.
Severe: $1-0.817574=0.182426$, reported as 0.1824.

The unrounded category probabilities sum to 1.000000. In a treatment profile with $\eta=-0.70$, the cumulative logits become $\theta_k-\eta=\theta_k+0.70$, giving cumulative probabilities 0.377541, 0.710950, and 0.900250. Differencing these unrounded values gives probabilities 0.3775 for none, $0.710950-0.377541=0.3334$ for mild, $0.900250-0.710950=0.1893$ for moderate, and $1-0.900250=0.09975$ for severe. These illustrative probabilities apply only to the specified covariate and random-effect profile; they are not population-averaged probabilities.

Why the Odds-Ratio Interpretation Holds

The treatment interpretation follows directly from the cumulative model. For a fixed threshold $k$, write $F_k=P(Y\leq k\mid b_i,\mathbf{x})$. The model states $\operatorname{logit}(F_k)=\theta_k-\eta$. Since $\operatorname{odds}(Y>k)=\{1-F_k\}/F_k$, its log is $-\operatorname{logit}(F_k)=\eta-\theta_k$. Comparing active treatment with usual care at the same visit, baseline pain, and random intercept changes $\eta$ by $\beta_{\text{treatment}}$, while $\theta_k$ stays fixed. Subtracting these log odds gives $\beta_{\text{treatment}}$; exponentiation gives the cumulative odds ratio $\exp(\beta_{\text{treatment}})$. The same argument holds for every threshold, which is precisely the proportional-odds restriction.

Diagnostics, Alternatives, and Limitations

A numerical model fit is not sufficient evidence that the analysis is appropriate. Start with the records: tabulate pain categories by visit and treatment, check the intended category order, and inspect the number of missing ratings and the patients represented at each visit. Sparse or absent categories can make threshold estimates unstable or prevent estimation. Verify that visit and treatment reference levels match the stated contrasts.

Then check the fit. Inspect convergence messages, the Hessian, threshold estimates, and the estimated random-intercept variance. Nonconvergence, very large standard errors, or an estimated variance near a boundary can indicate weak information, separation-like patterns, or an unnecessarily complex model. A patient-level random-intercept model also assumes that this single source of subject heterogeneity adequately represents repeated dependence. A proportional-odds model may fit poorly if treatment or visit effects differ substantially across pain thresholds; compare observed category patterns with fitted patterns and consider whether a less restrictive ordinal specification is justified.

The mixed model is not the only approach. A population-averaged ordinal model targets marginal rather than patient-conditional associations and handles within-patient correlation differently. Separate ordinal models at each visit avoid a common treatment effect but discard the joint repeated structure and can yield less coherent comparisons. A linear mixed model treats category scores as approximately equally spaced numeric measurements, an assumption that may not be clinically defensible. The ordinal mixed model is useful when ordered categories and patient-specific heterogeneity are central, but its conditional odds ratio is not a marginal probability difference and should not be reported as one.

Other limitations deserve attention. A random intercept does not allow the treatment effect to vary among patients; adding random slopes requires enough repeated information and may complicate estimation. The common treatment coefficient here imposes no change in the treatment association across visits on the cumulative-logit scale. If that scientific assumption is implausible, a treatment-by-visit interaction is a different model and requires a different interpretation. Finally, fitted probabilities for a patient with $b_i=0$ are not automatically averages over the patient population; obtaining population-level probabilities requires integrating over the random-effects distribution.

Common pitfall: A negative treatment coefficient does not mean that the treatment lowers every patient’s pain by that number of categories. It means the model estimates lower conditional cumulative odds of being above each pain threshold. Interpret $\exp(\beta)$ as an odds ratio, and state the category threshold and conditioning context.

Check Your Understanding

Answer each question, then compare with the solutions below.

Why should pain be stored as an ordered factor rather than as arbitrary category numbers in this analysis?
If $\hat{\beta}_{\text{treatment}}=0.40$, what is the estimated odds ratio for being above each pain threshold, and what direction does it indicate?
Given $P(Y\leq\text{none})=0.20$ and $P(Y\leq\text{mild})=0.55$, what is the probability of mild pain?
What does the proportional-odds assumption say about the treatment coefficient?

Solutions. (1) The factor preserves the clinical order without assuming that adjacent categories have equal numerical spacing. (2) $\exp(0.40)=1.492$. The estimated conditional odds of being above any threshold are about 1.49 times as large under active treatment, indicating a shift toward higher pain categories. (3) $P(Y=\text{mild})=P(Y\leq\text{mild})-P(Y\leq\text{none})=0.55-0.20=0.35$, or 35%. (4) It says that the treatment coefficient, and hence its cumulative odds ratio, is the same at each division between adjacent pain categories.