Why Compare Chains?
MCMC convergence is not established merely because an algorithm has produced many draws. A useful diagnostic asks whether independently initialized chains appear to be sampling the same posterior distribution. Rank-normalized split Rhat is one such comparison: it measures whether variation between chains is large relative to variation within chains, after splitting chains and transforming draws using their pooled ranks.
In Bayesian NLME models, the diagnostic can be computed separately for each scalar quantity of interest: for example, a population clearance parameter, a residual-error parameter, or a random-effect standard deviation. A reassuring value is close to one, but that is evidence about agreement among the chains, not proof that they have explored the entire posterior or generated enough effectively independent information.
Splitting Chains and Comparing Variances
Suppose \(m\) chains each contain \(N\) retained draws for one scalar parameter. For simplicity, take \(N\) to be even. Split each chain into its first and second halves. This gives \(2m\) split chains, each containing \(n=N/2\) draws. The split is useful because a chain can have an early portion that differs from its later portion even when summaries of the whole chain look similar to those of other chains.
For a set of \(m\) equal-length chains \(x_{ij}\), where \(i=1,\ldots,m\) indexes chains and \(j=1,\ldots,n\) indexes draws, define the mean of chain \(i\) as \(\bar{x}_i\), and the grand mean as \(\bar{x}_{\cdot}\). Let \(s_i^2\) be the sample variance within chain \(i\). The within-chain variance \(W\), between-chain variance \(B\), variance estimate \(\widehat{V}\), and conventional Rhat are:
These equations describe the basic variance comparison used in Rhat. Rank-normalized Rhat applies it to transformed draws rather than directly to the parameter values. In every calculation, the chains must have equal length; splitting supplies equal-length chains when the retained lengths are even. Implementations may specify how to handle an odd number of draws, so the reported diagnostic should be understood according to the software’s documented procedure.
Rank and Folded Transformations
A direct variance comparison can be unreliable when the sampled parameter distribution is strongly skewed or heavy-tailed. Rank normalization reduces sensitivity to the magnitude of extreme values. Pool the draws from all split chains and assign each draw its rank \(r\) among the total of \(S\) draws. If values are tied, use their average rank. Convert each rank to a normal score using the Blom plotting-position transformation:
Here \(\Phi^{-1}\) is the quantile function of the standard normal distribution. The normal scores preserve the ordering of the pooled draws while replacing their original spacings with a more stable scale. Compute the variance comparison on these scores to obtain rank-normalized split Rhat.
A difference in chain spread can be important even when chain centers agree. To help detect it, compute the pooled median of the original draws, take each draw’s absolute deviation from that median, rank-normalize those deviations by the same procedure, and compute Rhat again. This is the folded rank-normalized comparison. The final diagnostic is the larger of the ordinary rank-normalized and folded values. A low ordinary value alone is therefore not enough if the folded value is appreciably higher.
Worked Calculations
Worked Example: A Rank-Normalized Rhat Near One
Suppose four hypothetical MCMC chains for a population clearance parameter are split into eight chains of \(n=100\) draws each. After rank normalization, the within-chain variance is \(W=1.000\), and the between-chain variance is \(B=1.010\). Assume the folded calculation is no larger.
First calculate the estimated variance:
Then:
The result is close to one, indicating little discrepancy in this variance comparison. It does not show that the chains explored every relevant region of the posterior. It also does not establish that the available draws are sufficient for precise posterior summaries.
Worked Example: A Noticeable Between-Chain Discrepancy
Consider a second hypothetical parameter, a random-effect standard deviation. For \(n=100\) draws in each split chain, suppose its rank-normalized summaries give \(W=1.000\) and \(B=2.000\). The folded Rhat is \(1.030\).
The variance estimate is:
Thus the ordinary rank-normalized comparison is:
This is approximately \(1.005\) to three decimal places. But the reported diagnostic is the larger of the ordinary and folded values: \(\max(1.00499,1.030)=1.030\). The higher folded value warns that the chains’ dispersion differs more than their rank-normalized centers suggest. The comparison calls for investigation rather than a claim that the parameter’s chains agree adequately.
Worked Example: Why Rhat Near One Is Not a Guarantee
Imagine a hypothetical posterior for a model parameter with two separated regions of high posterior density. Four chains are all initialized near the first region and remain there. Their rank-normalized split-chain summaries happen to give \(W=1.000\), \(B=1.000\), with \(n=80\). Suppose the folded calculation is also 1.000.
The calculation is:
The diagnostic shows agreement among these chains, but all four can agree while missing the second region. Their common initialization and failure to move between regions make this a plausible failure mode. Trace plots or other evidence of poor exploration may reveal the problem; Rhat by itself cannot. The same caution applies if chains are highly autocorrelated: agreement between chains does not mean that the draws provide a large amount of independent information.
Justification and Practical Interpretation
The variance comparison has a simple motivation. If all chains are sampling the same stationary distribution, their means should not be systematically separated, and within-chain variation should be compatible with variation across chains. In that case \(B\) and \(W\) should be of similar scale, so \(\widehat{V}/W\) should be near one. If chain means differ, the squared deviations in \(B\) increase, which raises \(\widehat{V}\) relative to \(W\). Taking the square root expresses the discrepancy on a scale where one represents approximate agreement.
This is a diagnostic rationale, not a proof of convergence. It assumes that the chains are being compared for the same scalar quantity, have equal lengths after splitting, and are suitable samples from the posterior target. If the chains have not reached the target distribution, the variance comparison may be misleading. Rank normalization makes the diagnostic less dependent on extreme numerical values, and folding adds sensitivity to differences in spread, but neither transformation can make unvisited posterior regions appear in the data.
Diagnostics, Alternatives, and Limitations
Interpret Rhat together with the sampling context and other diagnostics. A large value is a warning to inspect the chains; it is not a measure of how far a parameter estimate is from its true value. For NLME models, inspect relevant population parameters, variance components, correlations, and derived quantities—not just one convenient parameter. The trace-plot material gives a visual way to look for drifting, separated, or persistently trapped chains.
| Approach | What it contributes | Important limitation |
|---|---|---|
| Rank-normalized split Rhat | Compares split-chain variation and is less sensitive to extreme value spacing; folding adds a check on spread. | Agreement can occur even if all chains miss a posterior region or sampling is inadequate. |
| Unsplit, unranked Rhat | Provides a basic between- versus within-chain variance comparison on the original values. | Can be affected by heavy tails, and using whole chains can conceal differences between their early and late portions. |
| Trace plots | Show chain behavior over iteration, including visible drift or separation. | Visual agreement can be difficult to assess and cannot by itself quantify the amount of useful sampling. |
| Other sampling diagnostics | Can assess aspects of sampling that chain agreement alone does not capture. | Each diagnostic addresses a particular concern; none replaces consideration of the model and sampling behavior. |
A practical failure mode is to treat a small Rhat as a pass/fail certificate. Instead, investigate elevated values, compare the rank and folded results, inspect trace plots, and consider whether the chains’ initial values and posterior exploration make their agreement credible. Also check whether the diagnostic is defined: if transformed draws have zero within-chain variance, the variance ratio has a zero denominator and is not a reassuring value of one. Such a result signals a degenerate chain summary requiring investigation.
Rhat is also a parameter-by-parameter summary. A reassuring diagnostic for a population clearance parameter does not automatically establish adequate sampling for a random-effect correlation or another derived quantity. Finally, even a well-behaved Rhat does not determine the precision of posterior estimates. It should be reported as one part of a broader assessment, with the number of chains, the monitored quantities, and any concerning behavior made clear.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- Four original chains each have 200 retained draws. How many split chains are formed, and how many draws are in each split chain?
- For \(n=50\), \(W=1.000\), and \(B=1.400\), calculate the variance estimate and Rhat.
- Why might rank normalization be preferable to calculating Rhat directly from the original parameter values?
- Suppose ordinary rank-normalized Rhat is 1.004 and folded rank-normalized Rhat is 1.025. What value is reported, and what is the main diagnostic concern?
- Give one reason that Rhat equal to one might coexist with inadequate posterior sampling.
Solutions. (1) Splitting each of four chains in half gives \(2(4)=8\) split chains, each with \(200/2=100\) draws. (2) \(\widehat{V}=(49/50)(1.000)+(1/50)(1.400)=0.9800+0.02800=1.0080\). Therefore \(\widehat{R}=\sqrt{1.0080/1.000}\approx1.00399\), or about 1.004 to three decimal places. (3) Rank normalization replaces the original value spacings with normal scores based on pooled ranks, reducing sensitivity to extreme values and heavy tails. (4) Report \(\max(1.004,1.025)=1.025\); the larger folded value suggests a difference in chain spread that deserves investigation. (5) All chains could remain in the same posterior region and miss another region, or they could agree while being highly autocorrelated. In either case, chain agreement alone does not establish adequate sampling.