Worked Workshop: Standardized Ordinal Probabilities

From an Ordinal GLMM to Standardized Probabilities

An ordinal GLMM describes ordered outcomes while allowing observations from the same cluster to be associated through random effects. Its regression coefficients describe the model conditional on those random effects. A population-averaged question instead asks: what is the probability of each outcome category after averaging over the random-effect distribution and a specified distribution of covariates?

This workshop answers that question for a three-category ordinal outcome and compares two exposure settings. The target covariate distribution is the empirical distribution of two observed profiles, given equal weight. The calculations illustrate the distinction between integrating over random effects and standardizing over covariates: both are required to obtain the stated population-averaged comparison.

Model and target: For observation $i$ in cluster $j$, let $Y_{ij}\in\{1,2,3\}$ be an ordered outcome, $A_{ij}$ a binary exposure, and $x_{ij}$ an observed covariate profile. Consider a cumulative-logit random-intercept model $$ \operatorname{logit}\{P(Y_{ij}\leq k\mid A_{ij},x_{ij},b_j)\} =\kappa_k-\{\beta_A A_{ij}+f(x_{ij})+b_j\},\qquad k=1,2, $$ where $\kappa_1<\kappa_2$, $b_j\sim N(0,\sigma_b^2)$, and $b_j$ is independent of the covariates under the fitted model. The two cutpoints and the slopes are shared across the cumulative logits (the proportional-odds assumption). The standardized probability for category $c$ at exposure setting $a$ is $$ \widehat P_c(a)=\sum_{m=1}^{M}w_m\int P(Y=c\mid A=a,x_m,b)\,dF_b(b), \qquad \sum_{m=1}^{M}w_m=1, $$ where $x_m$ are target covariate profiles and $w_m$ are their target-population weights. This averages response probabilities, not linear predictors.

The model assumes that the cumulative probabilities are ordered for every covariate and random-effect value, that the conditional outcome distribution is correctly represented by the cumulative-logit form, and that the random intercept has the stated normal distribution. The random-effect independence assumption specifies how the fitted model averages over $b_j$ at a covariate profile. Standardizing over observed profiles then defines the covariate target; it does not, by itself, make an observational exposure contrast causal.

Integrating the Random Intercept

For a given exposure and covariate profile, write the two conditional cumulative probabilities as $F_1(b)=P(Y\leq1\mid A,x,b)$ and $F_2(b)=P(Y\leq2\mid A,x,b)$. The conditional category probabilities are obtained by differencing adjacent cumulative probabilities:

P(Y=1\mid A,x,b)=F_1(b),\quad P(Y=2\mid A,x,b)=F_2(b)-F_1(b),\quad P(Y=3\mid A,x,b)=1-F_2(b).

The same identities hold after integration because integration is linear. If $\bar F_k(A,x)=\int F_k(b)\,dF_b(b)$, then the integrated category probabilities are $\bar F_1$, $\bar F_2-\bar F_1$, and $1-\bar F_2$. In particular, it is valid to integrate the cumulative probabilities and then difference them; equivalently, one can integrate each category probability directly.

For a normal random intercept, the integral generally has no elementary closed form. A common approximation is three-point Gauss–Hermite quadrature. For $b\sim N(0,\sigma_b^2)$, it evaluates a function $g(b)$ approximately as

\int g(b)\,dF_b(b)\ \approx\ \frac{1}{6}g(-\sqrt{3}\sigma_b)+\frac{2}{3}g(0)+\frac{1}{6}g(\sqrt{3}\sigma_b).

The values of $g$ are response-scale probabilities at the quadrature nodes. In practice, more nodes or adaptive integration can be used to check whether the approximation is sufficiently accurate. Once probabilities have been integrated for each profile, the covariate-standardization step is a weighted average across profiles.

Worked Example: Recovering Category Probabilities After Integration

Suppose a hypothetical fitted ordinal GLMM gives the following quadrature-integrated cumulative probabilities for exposure $A=1$ at one observed covariate profile: $\bar F_1=0.1551$ and $\bar F_2=0.5249$. These are probabilities after integrating over the random-intercept distribution, not probabilities conditional on a particular cluster effect.

For category 1, the probability is $\bar F_1=0.1551$. For category 2, subtract the first cumulative probability from the second:

\bar P_2=\bar F_2-\bar F_1=0.5249-0.1551=0.3698.

For category 3, take the remaining probability:

\bar P_3=1-\bar F_2=1-0.5249=0.4751.

The category probabilities sum to $0.1551+0.3698+0.4751=1.0000$. The ordering $0.1551\leq0.5249$ also confirms that the two integrated cumulative probabilities are ordered, as required.

Standardizing Over Observed Covariates

For a standardized contrast, select the target covariate profiles and their weights first. Then set the exposure to each comparison value for every profile, integrate over the random effects, and average across the same profiles and weights. Using the same target distribution at both exposure settings makes the contrast a comparison for a common covariate population.

The following hypothetical quadrature results summarize a fitted model. There are two observed profiles, $x_1$ and $x_2$, each assigned weight $1/2$. The entries are category probabilities after integration over the random intercept; each row therefore sums to one up to rounding.

Exposure $A$	Profile	Category 1	Category 2	Category 3
0	$x_1$	0.2925	0.4150	0.2925
0	$x_2$	0.2090	0.3994	0.3916
1	$x_1$	0.2243	0.4049	0.3708
1	$x_2$	0.1551	0.3698	0.4751

Worked Example: Standardized Exposure Contrast

Use the two profiles above with equal weights to compare the standardized probabilities at $A=1$ with those at $A=0$. For category 1 at $A=0$, for example, the average is

\widehat P_1(0)=\frac{0.2925+0.2090}{2}=0.25075\approx0.2508.

Applying the same calculation to each category and exposure setting gives

\begin{aligned} \widehat P_1(0)&=(0.2925+0.2090)/2=0.2508,\\ \widehat P_2(0)&=(0.4150+0.3994)/2=0.4072,\\ \widehat P_3(0)&=(0.2925+0.3916)/2=0.3421,\\ \widehat P_1(1)&=(0.2243+0.1551)/2=0.1897,\\ \widehat P_2(1)&=(0.4049+0.3698)/2=0.38735\approx0.3874,\\ \widehat P_3(1)&=(0.3708+0.4751)/2=0.42295\approx0.4230. \end{aligned}

Define the category-specific contrast as the probability at $A=1$ minus the probability at $A=0$. The results are

\begin{aligned} \Delta_1&=0.1897-0.2508=-0.0611,\\ \Delta_2&=0.38735-0.4072=-0.01985\approx-0.0199,\\ \Delta_3&=0.42295-0.34205=0.08090\approx0.0809. \end{aligned}

In this hypothetical target population, setting exposure to 1 rather than 0 is associated with about a 6.11 percentage-point lower probability of category 1, a 1.99-point lower probability of category 2, and an 8.09-point higher probability of category 3. These are differences in standardized probabilities; the calculation alone does not establish a causal effect.

Because the three category probabilities sum to one under each exposure setting, their contrasts must sum to zero. Here, using unrounded values, $-0.06105-0.01985+0.08090=0$. This is a useful arithmetic check: increasing probability in some categories must be balanced by decreases in others.

Why the Target Weights Matter

The equal weights in the example define a population with half its profiles like $x_1$ and half like $x_2$. Different weights answer a different population question, even with the same fitted model and the same integrated profile probabilities. If the observed target population instead gives $x_1$ weight $0.75$ and $x_2$ weight $0.25$, its standardized probabilities are not the equal-weight averages.

Worked Example: Changing the Target Covariate Distribution

Keep the same exposure-specific integrated probabilities, but use weights $0.75$ for $x_1$ and $0.25$ for $x_2$. At $A=0$, the category 1 probability is

\widehat P_1(0)=0.75(0.2925)+0.25(0.2090)=0.219375+0.05225=0.271625\approx0.2716.

The remaining $A=0$ probabilities are

\begin{aligned} \widehat P_2(0)&=0.75(0.4150)+0.25(0.3994)=0.4111,\\ \widehat P_3(0)&=0.75(0.2925)+0.25(0.3916)=0.317275\approx0.3173. \end{aligned}

At $A=1$, the corresponding probabilities are

\begin{aligned} \widehat P_1(1)&=0.75(0.2243)+0.25(0.1551)=0.2070,\\ \widehat P_2(1)&=0.75(0.4049)+0.25(0.3698)=0.396125\approx0.3961,\\ \widehat P_3(1)&=0.75(0.3708)+0.25(0.4751)=0.396875\approx0.3969. \end{aligned}

The corresponding contrasts are $0.2070-0.271625=-0.064625$, $0.396125-0.4111=-0.014975$, and $0.396875-0.317275=0.079600$. Rounded to four decimal places, these are $-0.0646$, $-0.0150$, and $0.0796$. They sum to zero before rounding. The difference from the equal-weight results is not a numerical error: the target population has changed.

Checks, Comparisons, and Limitations

A standardized ordinal probability is a population summary built from two averages: first over random effects, then over covariate profiles. It differs from a prediction at a representative profile, which does not average over the target covariate distribution. It also differs from a conditional prediction at a selected value of $b_j$, which describes an outcome conditional on that cluster effect rather than averaging across the random-effect distribution. The distinction between conditional and marginal results is essential because averaging a nonlinear response function generally does not give the response function evaluated at an averaged linear predictor.

Several checks help detect mistakes or model problems:

Category probabilities: each integrated profile row and each standardized vector should have probabilities between zero and one and sum to one, subject to rounding.
Cumulative ordering: $\bar F_1\leq\bar F_2$ must hold. A violation indicates an implementation or model problem.
Quadrature accuracy: increase the number of quadrature nodes or compare with adaptive integration. Material changes in probabilities or contrasts indicate that the original approximation was inadequate.
Target distribution: check that the profile weights sum to one and represent the intended population. For clustered data, explicitly decide whether the target weights observations or clusters if their sizes differ.
Ordinal-model fit: assess whether the proportional-odds structure is reasonable. If exposure or covariate effects differ across cumulative cutpoints, the fitted model may not represent the data well.
Random-effects assumptions: assess sensitivity to the assumed random-effect distribution and its estimated variance. The sensitivity of marginal effects to random-effects variance is a separate issue from the arithmetic of standardization.

The main alternative summaries answer different questions. Conditional coefficients describe changes on the cumulative-logit scale at fixed random effects. Marginalized predictions average over random effects but may be reported at selected covariate profiles. Standardization additionally averages those response-scale predictions over a stated target covariate distribution. None is universally preferable; the appropriate result depends on the scientific population and estimand.

These calculations are plug-in estimates: they use fitted model parameters and treat them as fixed. Inference for the standardized contrasts must also account for parameter uncertainty, for example through a delta-method calculation or simulation-based uncertainty procedure. The numerical contrasts above therefore describe point estimates, not confidence intervals or hypothesis tests. Finally, standardization over observed covariates does not resolve confounding, lack of exposure overlap, or misspecification; causal interpretation requires additional design and identification assumptions.

Common pitfall: A tempting shortcut is to average the fitted linear predictors across profiles and random effects, then apply the ordinal inverse-link once. This is generally wrong because the inverse-link is nonlinear: in general, $E\{g^{-1}(\eta)\}\ne g^{-1}\{E(\eta)\}$. The correction is to calculate category probabilities at each random-effect integration point, average those probabilities over the random effects, and only then standardize them over the target covariate profiles.

Check Your Understanding

Use the definitions and hypothetical results above. Answers follow each question so that the numerical and interpretive checks are explicit.

If integrated cumulative probabilities are $\bar F_1=0.22$ and $\bar F_2=0.61$, what are the three integrated category probabilities?
Solution: They are $0.22$, $0.61-0.22=0.39$, and $1-0.61=0.39$. They sum to 1.
With equal weights on $x_1$ and $x_2$, what is the standardized category 3 probability at $A=0$ using the table?
Solution: $(0.2925+0.3916)/2=0.34205$, or approximately $0.3421$.
Why must the three category-specific exposure contrasts sum to zero when the same profiles and weights are used at both exposure settings?
Solution: The category probabilities sum to one at each setting, so $\sum_{c=1}^3\{\widehat P_c(1)-\widehat P_c(0)\}=1-1=0$.
Does a standardized contrast from an observational GLMM automatically have a causal interpretation?
Solution: No. It is a model-based, standardized association unless causal identification assumptions, such as adequate control of confounding and exposure overlap, are also justified.
What is wrong with applying the inverse-logit to the mean linear predictor to obtain a marginal ordinal probability?
Solution: The inverse-link is nonlinear, so the inverse-link of a mean predictor generally differs from the mean of inverse-link probabilities. Integrate the response-scale category probabilities first, then standardize over covariates.

Exposure \(A\)	Profile	Category 1	Category 2	Category 3
0	\(x_1\)	0.2925	0.4150	0.2925
0	\(x_2\)	0.2090	0.3994	0.3916
1	\(x_1\)	0.2243	0.4049	0.3708
1	\(x_2\)	0.1551	0.3698	0.4751