Why a Point-Estimate Ranking Is Not the Whole Answer
A table of empirical Bayes predictions can make groups appear to have a clear order: sort the predicted random effects from largest to smallest. But each prediction is uncertain. If two groups have similar predictions relative to that uncertainty, their order may change often across plausible values of the effects. A ranked list based only on posterior means can therefore give a sharper impression of separation than the model supports.
The relevant question is not only which predicted effect is larger on average, but also how often one group’s effect exceeds another’s under the model’s uncertainty. This tutorial focuses on pairwise comparisons and rank probabilities for predicted group effects. It does not turn a ranking into a claim that group membership caused an outcome difference.
In an empirical Bayes analysis, predictions are commonly formed using estimated model parameters. A probability calculated from a posterior distribution conditional on those estimates reflects uncertainty about the latent effects conditional on the fitted parameters. It does not automatically include uncertainty in the estimated fixed effects or variance components. The earlier material on empirical Bayes prediction and latent-effect versus parameter uncertainty distinguishes these targets; the ranking calculation must be explicit about which uncertainty it uses.
Pairwise Comparisons and Rank Probabilities
For each pair of groups, define the contrast \(D_{jk}=b_j-b_k\). Then \(P(b_j>b_k\mid y)=P(D_{jk}>0\mid y)\). This comparison uses the joint uncertainty in the two effects. It is not generally valid to calculate it from two separate standard errors while ignoring their dependence.
If the posterior distribution of the effects is jointly normal, with posterior means \(m_j,m_k\), variances \(v_j,v_k\), and covariance \(c_{jk}\), then the difference is normal. Its mean is \(m_j-m_k\), and its variance is \(v_j+v_k-2c_{jk}\). Consequently:
Here \(\Phi\) is the standard normal cumulative distribution function. The formula is a special case; binary GLMMs and other models need not yield normal conditional distributions for their random effects. In those settings, the probability is still defined, but it is usually calculated from joint simulation draws rather than by imposing a normal approximation.
If there are no ties, the expected rank of group \(j\), counting rank 1 as highest, follows from counting how many groups exceed it:
This expected rank summarizes all pairwise comparisons involving a group. It is not the same as ranking groups by their posterior means, and it does not say that the group has that rank with certainty. The full rank distribution, including the probability of ranking first, may be more informative when the practical question concerns the top group.
Worked Calculations
Worked Example: A Normal Pairwise Comparison
Suppose two hypothetical groups have conditional posterior distributions \(b_A\mid y\sim N(0.60,0.25^2)\) and \(b_B\mid y\sim N(0.40,0.30^2)\). Assume their posterior covariance is zero. The effects are on a log-odds scale, and the goal is to assess whether group A’s random intercept exceeds group B’s.
The posterior mean difference is \(0.60-0.40=0.20\). Under the stated independence assumption, the variance of the difference is \(0.25^2+0.30^2=0.0625+0.0900=0.1525\), so its standard deviation is \(\sqrt{0.1525}=0.3905\). Thus:
Conditional on the fitted model and the stated uncertainty distribution, A’s random intercept exceeds B’s in about 69.6% of the posterior distribution. A is more likely than not to have the larger effect, but the probability is far from certainty. A normal 95% posterior interval for \(b_A-b_B\) is \(0.20\pm1.96(0.3905)=(-0.565,0.965)\) log-odds units, which includes zero. Neither the probability nor this interval establishes a difference in causal effects.
Worked Example: Ranking from Joint Posterior Draws
Suppose ten hypothetical joint draws of three groups’ random intercepts, in log-odds units, are summarized as follows. The first five rows represent five draws with the same values; the last five represent five other draws. Ties do not occur.
| Draws | \(b_A\) | \(b_B\) | \(b_C\) | Highest group |
|---|---|---|---|---|
| 1–5 (each draw) | 0.90 | 0.20 | 0.10 | A |
| 6–10 (each draw) | 0.12 | 0.58 | 0.66 | C |
The posterior means from these draws are \(E(b_A\mid y)\approx[5(0.90)+5(0.12)]/10=5.10/10=0.510\), \(E(b_B\mid y)\approx[5(0.20)+5(0.58)]/10=3.90/10=0.390\), and \(E(b_C\mid y)\approx[5(0.10)+5(0.66)]/10=3.80/10=0.380\). Ranking these means gives A, then B, then C.
The draws give a different picture of certainty. A is highest in 5 of 10 draws, so its estimated top-rank probability is \(5/10=0.50\); C is highest in the other 5, also giving \(0.50\). B is never highest in this small illustration. A exceeds B in 5 of 10 draws, and B exceeds C in 5 of 10 draws, so both corresponding pairwise probabilities are estimated as 0.50. The posterior means order B above C, yet the draws show each of those two effects larger than the other in half the draws.
These ten draws illustrate the calculation, not a precise estimate of a probability. For an estimated probability of 0.50 based on ten independent draws, the Monte Carlo standard error is approximately \(\sqrt{0.50(1-0.50)/10}=0.158\). In actual computation, enough well-mixed draws are needed to make simulation error acceptably small; correlated simulation draws require attention to effective sample size.
Worked Example: Why the Joint Distribution Matters
For two hypothetical groups, suppose the conditional posterior means are \(m_A=0.60\) and \(m_B=0.40\), each posterior variance is \(0.25\), and their posterior covariance is \(0.20\). These values imply posterior correlation \(0.20/\sqrt{0.25(0.25)}=0.80\), so the two predictions are strongly positively associated.
The difference has mean \(0.60-0.40=0.20\), but its variance uses the covariance:
Therefore \(P(b_A>b_B\mid y)=\Phi(0.20/0.3162)=\Phi(0.6325)\approx0.736\). If the covariance were incorrectly set to zero, the calculation would instead use variance \(0.25+0.25=0.50\) and give \(\Phi(0.20/\sqrt{0.50})=\Phi(0.2828)\approx0.611\). Both calculations have the same marginal means and variances; only the first uses the stated joint posterior distribution. This example shows why pairwise uncertainty cannot in general be recovered from marginal summaries alone.
Justification and Practical Computation
The normal contrast formula follows directly from properties of a multivariate normal distribution. If \((b_j,b_k)\mid y\) is jointly normal, any linear combination is normal; \(D_{jk}\) is the linear combination \(b_j-b_k\). Its expectation is \(E(b_j\mid y)-E(b_k\mid y)=m_j-m_k\). Expanding its variance gives:
Standardizing this normal difference gives the displayed probability formula. This argument requires a jointly normal posterior for the pair and the correct covariance; it does not require the posterior means to be unbiased estimates of fixed group-specific constants.
For simulation-based inference, obtain \(M\) joint draws \(b^{(1)},\ldots,b^{(M)}\) from the uncertainty distribution being used. Then calculate:
Here \(I(\cdot)\) equals 1 when its condition is true and 0 otherwise. Each draw must preserve dependence among groups. For a continuous distribution, exact ties have probability zero; if a model or reported values permit ties, a tie convention must be stated and applied consistently. If uncertainty in model parameters is part of the intended target, the draws must also reflect that uncertainty rather than treating estimated parameters as fixed.
Specify which random effects are compared, the direction of “higher,” and whether probabilities are conditional on fitted parameters or include parameter uncertainty.
Calculate contrasts from the joint distribution or from joint draws. Do not infer pairwise probabilities from separate marginal intervals alone.
Give the posterior-mean order if useful, but accompany it with pairwise or rank probabilities that show how uncertain that order is.
Diagnostics, Alternatives, and Limitations
Ranking is most useful when its uncertainty is reported, not hidden. Before interpreting probabilities, check that the fitted model and computation support them. Inspect convergence and mixing of posterior or simulation draws, especially for group effects with sparse observations. Poor mixing can make the estimated frequency of a pairwise ordering unreliable. Compare conclusions across reasonable modeling choices when the result is consequential, and assess whether the fitted model represents the observed data adequately. A ranking can be numerically precise under a model that is a poor description of the data.
Several shortcuts answer different questions or can mislead:
- Sorting posterior means. This provides a deterministic summary order, not the probability that each ordering holds. Small differences between means do not imply stable rankings.
- Comparing separate intervals. Overlap or non-overlap of marginal intervals is not itself a test of a pairwise contrast. The uncertainty of \(b_j-b_k\) depends on the joint distribution.
- Using raw group outcomes. Raw rates or averages are not the same as model-based random-effect predictions. They may be especially unstable for groups with little data; the model predictions also depend on the specified GLMM.
- Reporting only a highest-ranked group. A top-ranked label omits whether its top-rank probability is 0.51 or 0.99. Those cases support very different levels of confidence in the ranking.
Pairwise probabilities and rank probabilities also have limitations. They describe uncertainty under a specified model and target; they do not establish that one group is meaningfully different from another in a practical sense. A high probability of being higher can accompany a negligible effect-size difference, while a low probability can reflect limited information rather than equality. Ranking many groups creates many comparisons, so a few extreme-looking orderings may arise even when most groups are similar. For a random-slope model, a group’s rank may depend on the covariate value at which the slope contributes; compare predictions at a clearly specified value rather than treating a slope as a universal group ranking.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- Two groups have independent normal posterior distributions with means 0.30 and 0.10 and standard deviations 0.20 and 0.20. Calculate the probability that the first effect is larger.
- In 200 joint draws, group R ranks first 124 times. Estimate its top-rank probability. What does the estimate mean?
- Why can separate marginal posterior intervals fail to determine a pairwise ranking probability?
- With rank 1 defined as highest and no ties, write the formula for a group’s expected rank in terms of pairwise probabilities.
Solutions. (1) The difference has mean \(0.30-0.10=0.20\) and standard deviation \(\sqrt{0.20^2+0.20^2}=\sqrt{0.08}=0.2828\). Thus \(P(b_1>b_2\mid y)=\Phi(0.20/0.2828)=\Phi(0.7071)\approx0.760\). (2) The estimate is \(124/200=0.620\). Under the uncertainty distribution represented by the draws, group R ranks first in an estimated 62.0% of draws; it is not a 62.0% probability that a null hypothesis is true. (3) The probability depends on the distribution of the difference, including the covariance between effects. Marginal intervals alone do not provide that joint information. (4) \(E[\operatorname{rank}(b_j)\mid y]=1+\sum_{k\ne j}P(b_k>b_j\mid y)\), because the rank is one plus the number of effects exceeding \(b_j\).