Marginal Category Probabilities In Ordinal GLMMs

From Conditional Cumulative Probabilities to Category Probabilities

An ordinal response records an outcome in ordered categories, such as low, moderate, and high. In an ordinal generalized linear mixed model, the probability of being at or below a category is often specified conditional on covariates and a cluster-specific random effect. A population-averaged category probability answers a different question: after allowing for the random-effects distribution, what is the probability of each response category at the specified covariate values?

The calculation has two operations. First integrate each conditional cumulative probability over the random effects. Then take adjacent differences of those marginal cumulative probabilities. Keeping these operations distinct is important: a cumulative probability refers to a set of categories, whereas a category probability refers to one category alone.

Definition: Let $Y$ take ordered values $1,\ldots,K$, let $x$ denote covariates, and let $b$ be a random-effects vector with specified distribution $G$. Write $C_k(x,b)=P(Y\leq k\mid x,b)$, with $C_0(x,b)=0$ and $C_K(x,b)=1$. The population-averaged cumulative probability is $$ \bar C_k(x)=\int C_k(x,b)\,dG(b). $$ The population-averaged probability of category $k$ is $$ \bar p_k(x)=\bar C_k(x)-\bar C_{k-1}(x),\qquad k=1,\ldots,K. $$

For example, a cumulative-logit random-intercept model may be written

P(Y\leq k\mid x,b)=\operatorname{logit}^{-1}\{\theta_k-x^\mathsf{T}\beta-b\},\qquad \theta_1<\cdots<\theta_{K-1},\quad b\sim N(0,\sigma_b^2).

Here $\operatorname{logit}^{-1}(u)=1/(1+\exp(-u))$, the $\theta_k$ are ordered thresholds, and the same covariate effect applies at each threshold. The random intercept is assumed to follow the stated normal distribution and to have mean zero; conditional on $x$ and $b$, the response follows the specified ordinal model. For random slopes, $b$ is a vector and the integral is over its joint distribution, including its covariance structure.

These assumptions make the conditional cumulative probabilities ordered in $k$. Integration preserves that ordering, so $\bar C_1(x)\leq\cdots\leq\bar C_{K-1}(x)$. The adjacent differences are consequently nonnegative. The model must also define the boundary values $C_0=0$ and $C_K=1$; they ensure the first and last category probabilities are included.

Why Adjacent Differences Work

The event $Y=k$ is the event $Y\leq k$ with the event $Y\leq k-1$ removed. Conditional on $x,b$, this gives $P(Y=k\mid x,b)=C_k(x,b)-C_{k-1}(x,b)$. Taking expectations over the random-effects distribution yields the marginal category probability:

\begin{aligned} P(Y=k\mid x) &=\int P(Y=k\mid x,b)\,dG(b)\\ &=\int \{C_k(x,b)-C_{k-1}(x,b)\}\,dG(b)\\ &=\bar C_k(x)-\bar C_{k-1}(x). \end{aligned}

The last equality follows from linearity of integration. Thus, integrating conditional category probabilities directly gives the same result as first integrating the cumulative probabilities and then differencing. The latter route is often convenient because ordinal models are commonly parameterized in terms of cumulative probabilities.

Worked Example: Differencing Four Marginal Cumulative Probabilities

Suppose a fitted ordinal model has four ordered response categories: 1 = low, 2 = lower-middle, 3 = upper-middle, and 4 = high. At a specified covariate value, numerical integration of the fitted conditional model gives $\bar C_1=0.18$, $\bar C_2=0.61$, and $\bar C_3=0.88$. These values are hypothetical and are assumed to come from a model with ordered thresholds and a valid random-effects distribution.

The boundary cumulative probabilities are $\bar C_0=0$ and $\bar C_4=1$. Therefore,

\begin{aligned} \bar p_1&=0.18-0=0.18,\\ \bar p_2&=0.61-0.18=0.43,\\ \bar p_3&=0.88-0.61=0.27,\\ \bar p_4&=1-0.88=0.12. \end{aligned}

The four probabilities sum to $0.18+0.43+0.27+0.12=1.00$, and each is nonnegative. At these covariate values, the population-averaged probability of the upper-middle category is 0.27, or 27%. This is an average over the assumed random-effects distribution, not a prediction for a cluster known to have random effect zero.

Integrating the Random Effects

For a logistic cumulative link with a normal random intercept, the integral generally has no elementary closed form. A standard numerical approach is Gaussian quadrature: approximate the integral by evaluating the conditional probability at selected random-effect values and applying weights. The weights sum to one, and the quadrature nodes and weights are chosen to approximate integration under the normal distribution. More nodes are usually needed when the random-effects variance is large or the integrand changes rapidly.

Worked Example: Logistic-Normal Category Probabilities by Quadrature

Consider a hypothetical three-category cumulative-logit model with thresholds $\theta_1=-1$ and $\theta_2=1$, linear predictor $\eta=0.35$, and random intercept $b\sim N(0,0.8^2)$. Thus $C_k=\operatorname{logit}^{-1}(\theta_k-\eta-b)$. Use the three-point normal quadrature approximation with nodes $b=-1.3856,0,1.3856$ and weights $1/6,2/3,1/6$. The thresholds are ordered, the random effect has the stated centered normal distribution, and the model’s cumulative probabilities are consequently ordered for every node.

For category 1, evaluate the first cumulative probability at each node:

\begin{aligned} C_1(-1.3856)&=\operatorname{logit}^{-1}(-1-0.35+1.3856)=0.5089,\\ C_1(0)&=\operatorname{logit}^{-1}(-1.35)=0.2059,\\ C_1(1.3856)&=\operatorname{logit}^{-1}(-2.7356)=0.0609. \end{aligned}

Weighting these values gives $\bar C_1\approx(1/6)(0.5089)+(2/3)(0.2059)+(1/6)(0.0609)=0.2322$. For the second cumulative probability, the node-specific values are $0.8844$, $0.6570$, and $0.3239$, respectively. Hence

\bar C_2\approx(1/6)(0.8844)+(2/3)(0.6570)+(1/6)(0.3239)=0.6394.

Including the boundaries, the marginal category probabilities are

\begin{aligned} \bar p_1&=0.2322,\\ \bar p_2&=0.6394-0.2322=0.4072,\\ \bar p_3&=1-0.6394=0.3606. \end{aligned}

They are nonnegative and sum to $1.0000$, up to rounding. These are quadrature approximations to probabilities under the model, not exact values from a closed-form integral. The three-point rule is useful for showing the calculation; in an analysis, one should check whether the numerical results stabilize when using a more accurate integration rule.

A Closed-Form Case: Probit Random Intercepts

A probit cumulative link has a useful exact result when the random intercept is normal. Suppose $P(Y\leq k\mid x,b)=\Phi(\theta_k-\eta-b)$, where $\Phi$ is the standard normal cumulative distribution function, $\eta=x^\mathsf{T}\beta$, and $b\sim N(0,\sigma_b^2)$. Equivalently, the cumulative event can be represented using an independent standard normal latent error $\epsilon$: it occurs when $\epsilon+b\leq\theta_k-\eta$. Since $\epsilon+b\sim N(0,1+\sigma_b^2)$,

\bar C_k(x)=\Phi\left(\frac{\theta_k-\eta}{\sqrt{1+\sigma_b^2}}\right).

The independence of the latent error and random intercept, along with their normal distributions, is essential to this result. The category probabilities still require adjacent differences, and are not obtained by applying $\Phi$ directly to a category-specific threshold.

Worked Example: Exact Marginal Probabilities in a Probit Model

Suppose a hypothetical three-category probit random-intercept model has thresholds $-0.5$ and $0.8$, linear predictor $\eta=0.3$, and $b\sim N(0,0.5^2)$. The thresholds are ordered, and the latent standard normal error is assumed independent of $b$. Here $\sqrt{1+\sigma_b^2}=\sqrt{1.25}=1.1180$. Thus,

\bar C_1=\Phi\left(\frac{-0.5-0.3}{1.1180}\right)=\Phi(-0.7155)\approx0.2371,\qquad \bar C_2=\Phi\left(\frac{0.8-0.3}{1.1180}\right)=\Phi(0.4472)\approx0.6726.

Taking adjacent differences gives $\bar p_1=0.2371$, $\bar p_2=0.6726-0.2371=0.4355$, and $\bar p_3=1-0.6726=0.3274$. Rounding accounts for the final decimal in the sum, which is approximately 1.000. For comparison, setting $b=0$ gives conditional cumulative probabilities $\Phi(-0.8)\approx0.2119$ and $\Phi(0.5)\approx0.6915$, and category probabilities approximately $0.2119,0.4796,0.3085$. The difference illustrates why a prediction at the mean random effect is not generally the population average.

Checks, Alternatives, and Common Pitfalls

Several checks help identify mistakes in either the model calculation or its numerical implementation:

Check cumulative ordering. The integrated values must satisfy $0\leq\bar C_1\leq\cdots\leq\bar C_{K-1}\leq1$. A violation suggests a coding, parameter, or integration error.
Check category probabilities. Each adjacent difference must be nonnegative, and the probabilities must sum to one, apart from rounding error.
Check integration accuracy. For numerical quadrature, increase the number of nodes or compare a second reliable integration method. Results that change materially need more accurate integration.
Check the random-effects distribution. The answer is averaged over the distribution specified by the model. A misspecified variance or distribution can change the marginal category probabilities.

The main alternative is to report a conditional prediction at a chosen random effect, often $b=0$. This describes the model conditional on that value; it is not generally population-averaged. For probit models under the normal assumptions above, the closed form avoids numerical quadrature. For logistic-normal models, numerical integration is usually required. Direct integration of the conditional category probabilities is also valid, but it should produce the same result as integrating cumulative probabilities and differencing them.

A common mistake is to calculate $\operatorname{logit}^{-1}(\theta_k-\eta)$ for each threshold and call the resulting differences marginal probabilities. That sets $b=0$; it does not integrate over $b$. Another mistake is to average conditional cumulative probabilities but forget to subtract adjacent values. The correction is to integrate the cumulative probabilities under the specified random-effects distribution, then form all $K$ differences using the boundary values 0 and 1. These probabilities describe model-based averages at fixed covariates; they do not by themselves account for uncertainty in estimated model parameters.

Key takeaway: Integrate each conditional cumulative probability over the random-effects distribution, then difference adjacent marginal cumulative probabilities. Verify ordering, nonnegativity, and a total probability of one.

Check Your Understanding

Answer each question, then compare with the solutions below.

If a three-category model gives $\bar C_1=0.24$ and $\bar C_2=0.70$, what are its three marginal category probabilities?
Why does evaluating a cumulative-logit model at $b=0$ not generally give population-averaged category probabilities?
State two assumptions behind the closed-form probit marginal cumulative probability.
A numerical calculation returns $\bar C_1=0.62$ and $\bar C_2=0.51$ for a three-category model. What is wrong, and what diagnostic reveals it?

Solutions. 1. The probabilities are $0.24$, $0.70-0.24=0.46$, and $1-0.70=0.30$; they sum to 1. 2. Setting $b=0$ gives a conditional prediction at one random-effect value. Because the inverse link is nonlinear, this is generally not equal to integrating predictions over the random-effects distribution. 3. The latent probit error and random intercept are independent, and both are normal, with latent error variance 1 and random-intercept mean zero. 4. The cumulative probabilities decrease from the first threshold to the second, contrary to the required ordering. The ordering check fails, so investigate the model parameterization or integration implementation before differencing.