Rank-Normalized R-Hat for Multiple Chains

Why Compare Multiple Chains?

Several chains can produce plausible-looking traces while exploring different parts of a target distribution. Comparing only their running means or autocorrelations may not reveal persistent disagreement in location or scale. Rank-normalized R-hat is a diagnostic that compares variation within chains with variation among chains, after splitting each chain and transforming pooled ranks to normal scores. The split helps reveal changes between the early and late portions of a run; the rank transformation makes the diagnostic less sensitive to the scale and tails of the original parameter.

Suppose \(M\) chains each provide \(N\) scalar draws \(X_{m,1},\ldots,X_{m,N}\), where \(m=1,\ldots,M\). Regard these draws as the output of \(M\) MCMC runs on a common probability space \((\Omega,\mathcal{F},\mathbb{P})\). The formulas below are applied to the recorded draws; they do not assume those draws are independent. For simplicity, take \(N\) even and divide each chain into its first and last \(N/2\) draws. This produces \(m=2M\) split chains, each of length \(n=N/2\). The letter \(m\) below denotes the number of split chains, not the number of original chains.

To rank-normalize, pool the \(S=mn\) draws from the split chains and assign each draw its pooled rank \(r\). For ties, use the average of the ranks the tied observations would occupy. Convert the rank to the score \(z=\Phi^{-1}((r-3/8)/(S+1/4))\), where \(\Phi\) is the standard normal distribution function. Then compute \(\hat R\) from the scores, not the original draws. This is the rank-normalized split R-hat.

A second calculation targets scale differences. Find the median of all pooled draws and replace each draw by its absolute deviation from that median. Rank-normalize these deviations using the same plotting-position formula, then calculate split R-hat again. This gives the folded rank-normalized R-hat. The reported rank-normalized R-hat is the larger of the ordinary rank-normalized and folded values. The folding is useful because chains can have similar centers but noticeably different spreads.

What the R-hat Formula Compares

Proof. Starting from the definition, divide the expression for \(\widehat{\operatorname{Var}}^{+}\) by \(W\), which is permitted because \(W>0\). This gives \(\hat R^2=(n-1)/n+B/(nW)\). Since \((n-1)/n=1-1/n\), rearranging gives \(\hat R^2=1+(B/W-1)/n\). Because \(\hat R\) is the nonnegative square root, comparing it with \(1\) is equivalent to comparing \(\hat R^2\) with \(1\). The displayed identity then gives the stated conditions. \(\square\)

The identity describes the calculation, not a convergence theorem. A value near one says that the estimated between-chain and within-chain variation are similar for these transformed draws. It does not establish that the chains have reached stationarity, explored every relevant region, or accurately estimated a desired quantity. A numerical cutoff is a screening convention, not proof.

Proof. Let \(g\) be strictly increasing. For any two pooled draws \(x\) and \(y\), \(x<y\) implies \(g(x)<g(y)\), and \(x>y\) implies \(g(x)>g(y)\). Therefore every pairwise ordering, and thus every pooled rank, is unchanged by applying \(g\). The score assigned to a rank depends only on that rank and on \(S\); it is consequently unchanged as well. The split-chain means, variances, and R-hat computed from those scores are then unchanged. \(\square\)

This rank invariance applies to the ordinary rank-normalized calculation. Folding is a separate operation on absolute deviations from the pooled median; its purpose is to compare spread as well as location, so it should not be interpreted as merely applying an increasing transformation to the original draws.

Rank Normalization in a Small Example

Folding to Detect Scale Disagreement

Interpretation and Common Pitfalls

For a scalar parameter, compute the ordinary rank-normalized statistic and the folded statistic, then report their maximum. For a vector parameter, apply the diagnostics to each scalar component of interest; a reassuring value for one component says nothing by itself about another. Rank normalization can reduce sensitivity to extreme magnitudes, while folding asks whether chains differ in their distance from a common center. Neither calculation repairs poor exploration or substitutes for inspecting traces and other diagnostics.

Common mistake: a small numerical difference in an intermediate calculation can materially change the final statistic. In the folded example, incorrectly using \(0.7025^2=0.493\) and then rounding too early can produce an inaccurate \(B\). Recompute from the displayed value: \(0.7025^2=0.49350625\), so \(2[4(0.7025^2)/3]=2(1.974025/3)=1.3160167\), not \(1.3174\). Using \(W=0.3112605\), the pooled estimate is \((0.3112605+1.3160167)/2\approx0.8136386\), and the folded statistic is approximately \(1.6168\). The correction is not to infer convergence from the revised value, but to calculate accurately and interpret the diagnostic cautiously.

Another pitfall is to treat a conventional cutoff, such as a value close to one, as a mathematical guarantee. The statistic is computed from finite output and summarizes particular comparisons of within-chain and between-chain variation. It cannot prove that the target distribution was reached or that all relevant regions were explored. Conversely, a high value is a reason to investigate the chains, not by itself an explanation of what went wrong. When all within-chain variances are zero, \(W=0\), so the ratio in the definition is undefined; do not report a finite R-hat by dividing by zero. Ties are handled through average ranks, and the same pooled ranks must be used for all split chains.