Skip to contents

MCMC convergence diagnostics

Source: R/diagnostics.R
mcmc_diagnostics.Rd

One table per parameter, a multivariate summary, a location summary and a verdict for the chains of gibbs() (or for any set of draws). Each number answers a different question, and the verdict reads them together; no diagnostic can prove convergence, it can only fail to find a problem.

Usage

mcmc_diagnostics(
  x,
  derived = TRUE,
  locations = TRUE,
  probs = c(0.05, 0.95),
  thresholds = list()
)

Arguments

x

a gibbs() fit; or draws: a numeric matrix iterations x chains (one parameter), an array iterations x chains x parameters (the layout of posterior::draws_array), a list of matrices iterations x parameters, one per chain (matched by column name), a draws object of the posterior package in any format (read through posterior::as_draws_array()), or an mcmc (one chain) or mcmc.list of the coda package. Raw draws need at least 4 per chain (with 4 or 5 the verdict is at best "run longer"); chains shorter than the longest are left out, and the verdict says so

derived

add the heritability and the correlations computed draw by draw (fits only)

locations

add the location R-hat (fits with two or more chains)

probs

the quantiles reported, with their MCSE

thresholds

a named list overriding any of rhat = 1.01, ess_per_chain = 100, disagreement = 1.5 (location and scale), isolated = 1.05, alpha = 0.01, khat = 0.5, each a single finite number

Value

an object of class breeding_mcmc_diag: table (one row per parameter, the columns above plus type, component or derived, status and note), multivariate, locations (one row per location term: quantiles of the level R-hat, share_above 1.01, tau_reference, limit, active, flagged, ess_typical, ess_low, and mcse_sd_q50 and mcse_sd_q90, the median and 90% quantile over the levels of the Monte Carlo error of the pooled mean in posterior sd, from the spread of the chain means), worst_levels, boundary (on a gibbs() fit: one row per covariance group and chain on the boundary, with group, cor for a group of two terms, chain, share and min_eigen, the smallest eigenvalue of the correlation matrix over the kept draws; no rows when no chain is there), verdict, reasons, notes, thresholds and info (chains, kept draws per chain, stopped chains, n_iter, burnin, thin, seconds and the smallest bulk-ESS per second).

Details

Per parameter (the sampled components and, with derived = TRUE, the heritability computed draw by draw as in summary() and the correlation of every cov() with its two variances, since a failure can show only in a nonlinear function of the components):

  • rhat: the rank-normalized split R-hat of Vehtari et al. (2021), the larger of the bulk and the folded (tail) versions; rhat_bulk and rhat_tail separately.

  • ess_bulk, ess_tail: the multi-chain effective sample sizes of the same paper. The chains share one autocorrelation estimate, so chains that disagree LOWER the ESS (a sum of per-chain ESS cannot see that). ess_bulk_chain_min is the smallest per-chain bulk-ESS. disagreement is the sum of the per-chain bulk-ESS divided by the multi-chain one: near 1 when the chains agree, large when the between-chain variance exceeds what their autocorrelation predicts. disagreement_tail is the same ratio on the folded draws |x - median|, which sees chains that agree in location and differ in scale.

  • mcse_mean, mcse_sd, mcse_q*: Monte Carlo standard errors of the posterior mean, sd and quantiles (Vehtari et al. 2021). They say whether the reported digits are supported, which is adequacy, not convergence. mcse_sd also assumes a finite fourth moment of the draws (Pareto k below 0.25), which is not checked.

  • acf1, acf10, acf50: autocorrelation at lags 1, 10 and 50 of the kept draws, averaged over the chains.

  • geweke_z, geweke_p, hw_p, hw_discard: tests of the START of the chains. R-hat cannot see a short, extreme initial transient (the ranks bound its influence and it inflates the within-chain variance); these two can. The Geweke (1992) z compares the mean of the first 10% with the last 50%, with the spectral density at zero estimated ONLY in the final window and pooled over the chains, and is combined across chains (the sum of the per-chain z over sqrt(M) and the sum of their squares, Bonferroni over the two). The Heidelberger and Welch (1983) Cramer-von Mises test sums the Brownian bridges of the chains; its first look is the test, and hw_discard is the smallest share of the kept draws (10% to 50%) whose removal passes. Both p-values carry the noise of the estimated spectral density through a t reference (a chi-square scale mixture for the Cramer-von Mises statistic) with a quarter of the effective size of the final window as degrees of freedom, a rule tuned by simulation and not derived: on a stationary single AR(1) chain whose first window held 4 to 21 effective draws (seven settings), the normal reference rejected 6% to 11% of the time at 5% (Geweke) and the plain first look 7% to 11% (Heidelberger-Welch); with it, 3% to 8% and 3.5% to 6% (400 replicates per setting, about 1 point of Monte Carlo error; validation/mcmc_diagnostics_calibration.R, section A). window_ess is the effective size of the first window, summed over the chains. These numbers differ from coda::geweke.diag(), which estimates the spectral density in each window separately and lets the transient inflate its own variance.

  • start_scale: the scale those two tests ran on. Their reference for the noise of the spectral density assumes draws with a finite fourth moment; on heavy-tailed stationary draws (inverse gamma with Pareto k 0.5 to 0.83, an AR(1) copula with such a margin, t with 1.5 df) the raw tests called the burn-in too short in 6.5% to 15% of the runs at the 1% level (400 replicates per setting, section J). So the draws stay "raw" while the Pareto k-hat of their second half (the reference window of both tests) is at most 0.25, a quantity whose draws are all positive goes to "log" above that (the scale Gelman and Rubin 1992 recommend for variances), and the normal scores of the ranks ("rank") are the last resort, noted in notes; on the same draws the verdict then said "burn-in too short" 0% to 2.2% of the time. The switch is limited to heavy tails because it costs power: a start 50 sd above a positive quantity with light tails was caught 86% (one chain) and 100% (4 chains) of the time on the raw scale, and would have been 36% and 98% in log and 0% and 14% on the ranks (50 replicates). On a heavy-tailed quantity a start far above the posterior is a handful of huge draws: the raw tests caught it every time, but they also flag the stationary runs above; in log it is a short transient, caught 6% of the time, which barely moves the median and quantiles that the heavy-tail note points to.

  • khat: the Pareto k-hat of the two tails (Vehtari et al. 2024), as posterior::pareto_khat(). Above 0.5 the draws may have infinite variance, the empirical mean and its MCSE are unreliable, and both MCSE columns are left blank (0.7, the threshold of posterior, is for Pareto-smoothed estimates, not for the raw mean). It is NA with 29 or fewer draws in all (each tail then has 5 draws or fewer, too few for the fit) and when both tails are constant; posterior returns -Inf in both cases.

Multivariate (multivariate): the multivariate ESS of Vats, Flegal and Jones (2019) by batch means replicated over the chains (Gupta and Vats 2020) with the lugsail correction (Vats and Flegal 2022), the minimum ESS they require (min_ess), the multivariate PSRF of Brooks and Gelman (1998) with their factor (1 + 1/M) for M chains (coda::gelman.diag() uses (1 + 1/p) for p parameters; the two agree only when M = p), the R-hat of every pair of chains (the larger over the parameters) and the chains whose smallest pair R-hat exceeds isolated. The multivariate numbers are information: the multivariate ESS can exceed the smallest univariate one (a fast direction compensates a slow one), so it does not enter the verdict. An isolated chain does, with the per-parameter rule below; it mostly fires together with a disagreement, and alone (above the ESS floor) in at most 3% of the runs of the simulated settings where a chain was shifted or rescaled.

Locations (locations, worst_levels): the classic R-hat of every breeding value and fixed effect from the per-chain means and variances gibbs() keeps. It does not split the chains, so it cannot see a trend within a chain, and its Monte Carlo error has only M - 1 degrees of freedom, so with 2 chains the limit is high: for 860 levels and tau_reference 40 it is 1.07 at 5000 kept draws per chain and 1.19 at 1800, against 1.03 and 1.09 with 4 chains (section E5), and 4 chains are the usual choice. With Q levels over all the location terms the largest R-hat grows with Q even at convergence, so a level is flagged only above max(rhat, R_Q), where R_Q^2 is the 1 - alpha/Q quantile of (n-1)/n + tau chi2_{M-1}/((M-1) n), the null distribution of a level whose autocorrelation time is at most tau = tau_reference. The levels do not share one autocorrelation time (that of a breeding value grows with its distance from zero in posterior sd), so tau_reference is the larger of the typical one (from the median over the levels of the term, ess_typical = M n / tau) and a ceiling. In the Gaussian sampler, which draws all locations given the components and then the components given the locations, the autocorrelation time of a location is at most that of a function of the components plus 2 (Liu, Wong and Kong 1994), and the slowest row of the table stands in for that function (with every component held, the locations are independent draws). The stand-in is a heuristic: a nonlinear function of several components can mix slower than each of them. The limit was calibrated on simulated chain means with a valid ceiling given (the true largest time plus 2); the ceiling as a fit computes it was checked on simulated data of 860 animals in three models, where the largest level time over that of the slowest row was 0.17 (animal model), 0.33 (direct-maternal, chains far from mixed) and 0.20 (a two-coefficient reaction norm with 5 records per animal; sections G and G2): no genomic model, no probit, nothing at scale. In the probit model the latent liabilities make a third block, the ceiling does not hold, and only the typical time is used, without calibration. A flagged level mixes slower than the ceiling allows or its chains disagree, and one case no run length fixes is a fixed effect whose posterior is improper, such as a covariate that separates the two categories in the probit model (a contemporary group with every record in one category is refused by gibbs() before the chain). ess_low says that the typical effective size of a term is below ess_per_chain per chain by a one-sided test at alpha (Bonferroni over the terms) on the median over its levels: with few levels that median is nearly one chi-square with M - 1 degrees of freedom, and the point estimate alone flagged 7% to 25% of the runs with independent levels in terms of 1 to 3 levels, where the test flagged none in 2000 (its power there is low too: 5% to 41% at a quarter to half of the floor; section E4). active is FALSE when the limit sits above sqrt((n-1)/n + 1/M), the R-hat of a level with one chain one posterior sd away from the others: the check could not flag even that level, which happens when a component mixes slowly, and a note says so instead of letting flagged = 0 read as a pass. worst_levels holds the ten largest location R-hat of each term; in a group with more than one effect per level (direct and maternal, indirect(), the coefficients of rn()) each name carries the effect in parentheses, "F016 (maternal)", as the term column of solutions().

The status of a row is the most severe condition it triggers, and every condition it triggers is written in reasons:

  • "not mixed": a chain that does not move while the others do (or a sampled column that moves in no chain); or R-hat at or above rhat with the chains apart beyond their autocorrelation (disagreement or disagreement_tail above its threshold, or a chain whose R-hat with every other exceeds isolated, from 3 chains) and at least ess_per_chain bulk draws in the weakest chain. Below that floor a disagreement cannot be told from noise: on stationary chains with 10 to 26 bulk draws per chain, without the floor, 12.5% to 61% of the runs would be "not mixed" from the disagreements alone and 37% to 97% counting the isolated chain (4 and 8 chains, section C); the row is "run longer" instead.

  • "burn-in": the Geweke or Heidelberger-Welch test rejects at family-wise level alpha (Bonferroni over the parameters and the two tests).

  • "run longer": bulk- or tail-ESS below ess_per_chain per chain (M chains need M times as many), fewer than 6 kept draws per chain (the halves are then too short to estimate an ESS, and the verdict would rest on the R-hat alone), the split R-hat of a single chain at or above rhat (its two halves disagree), or an R-hat at or above rhat across chains with no chain set apart, which the ESS explains: R-hat^2 is about 1 + 2M/ESS even at convergence, and at 100 bulk draws per chain more than a third of converged runs pass 1.01 (37% to 42% with 2, 4 and 8 chains, section H).

  • "ok", or "fixed" for a column the model holds: a kernel(fixed =) component, the probit residual, every component under theta_fixed, and a derived row when every component is held. For raw draws, "fixed" is a column constant on the scale of its values.

The verdict is the worst row, in the order "chain stopped" > "not mixed" > "burn-in too short" > "run longer" > "ok"; a flagged location, or a location term with ess_low, makes it at least "run longer". On a gibbs() fit, a covariance group (two or more terms joined by a cov(): direct and indirect, direct and maternal, the coefficients of rn()) with more than 1% of the kept draws of a chain within 1e-6 of a singular correlation matrix (its smallest eigenvalue; for two terms, |cor| > 0.999999) makes it at least "not mixed", with that reason first (after a stopped chain, if any): the chain sits on the boundary of the covariance matrices, where the default prior of gibbs() gives the posterior infinite mass, and a longer chain does not leave it (see prior in gibbs()). Under a posterior with no infinite mass there, and whose correlation is not itself at +-1, the share of draws that close to the boundary is of the order of 1e-6. The groups and chains are in boundary. "chain stopped" means a chain of gibbs() stopped early (an inadmissible theta or a failed Cholesky inside the chain): it is left out of the pooled results and of every diagnostic here. With one chain (the default of gibbs()) the R-hat is the split R-hat of its two halves, and the between-chain checks (disagreement, isolated chain, location R-hat) are inactive: run chains = 4 for them.

The thresholds are rules of thumb, not calibrated error rates, with three exceptions measured by simulation, on stationary Gaussian AR(1) chains and, for the burn-in tests, also on heavy-tailed draws: the burn-in tests at alpha, the location limit and the ess_low test. What the verdict did on simulated chains (validation/mcmc_diagnostics_calibration.R, sections B to D, J and K): with three independent parameters on stationary Gaussian chains of 200 to 600 bulk draws per chain (200 replicates per setting), "not mixed" never came out with 1, 2, 4 or 8 chains and "burn-in too short" at most 2% of the time; at 10 to 53 per chain, never "not mixed"; on the heavy-tailed draws of section J, "burn-in too short" at most 2.2% and never "not mixed"; with 4 to 7 draws per chain, never "ok". With one of 4 AR(1) chains shifted by 0.5 and by 1 posterior sd, "not mixed" 83% and 100% of the time, while a shift of 0.25 sd was not told from noise (2%); with one of 2 iid chains at 1.5 or 2 times the sd of the other, 100% (69% at 1.25 times; 100 replicates). A start 50 sd away with no burn-in gave "burn-in too short" every time with 4 chains and 86% of the time with one (50 replicates), where R-hat passed 1.01 in 2% of the 4-chain runs.

References

Vehtari, A., Gelman, A., Simpson, D., Carpenter, B. & Burkner, P.-C. (2021). Rank-normalization, folding, and localization: an improved R-hat for assessing convergence of MCMC. Bayesian Analysis 16:667-718.

Vats, D., Flegal, J.M. & Jones, G.L. (2019). Multivariate output analysis for Markov chain Monte Carlo. Biometrika 106:321-337.

Vats, D. & Flegal, J.M. (2022). Lugsail lag windows for estimating time-average covariance matrices. Biometrika 109:735-750.

Gupta, K. & Vats, D. (2020). Estimating Monte Carlo variance from multiple Markov chains. arXiv:2007.04229.

Brooks, S.P. & Gelman, A. (1998). General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics 7:434-455.

Gelman, A. & Rubin, D.B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science 7:457-472.

Vats, D. & Knudson, C. (2021). Revisiting the Gelman-Rubin diagnostic. Statistical Science 36:518-529.

Geweke, J. (1992). Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. In Bayesian Statistics 4. Oxford University Press.

Heidelberger, P. & Welch, P.D. (1983). Simulation run length control in the presence of an initial transient. Operations Research 31:1109-1144.

Anderson, T.W. & Darling, D.A. (1952). Asymptotic theory of certain "goodness of fit" criteria based on stochastic processes. Annals of Mathematical Statistics 23:193-212.

Vehtari, A., Simpson, D., Gelman, A., Yao, Y. & Gabry, J. (2024). Pareto smoothed importance sampling. Journal of Machine Learning Research 25(72):1-58.

Zhang, J. & Stephens, M.A. (2009). A new and efficient estimation method for the generalized Pareto distribution. Technometrics 51:316-325.

Liu, J.S., Wong, W.H. & Kong, A. (1994). Covariance structure of the Gibbs sampler with applications to the comparisons of estimators and augmentation schemes. Biometrika 81:27-40.

Gu, C.C., Lakatos, E. & Hamis, S. (2026). The MCMC convergence graph: a diagnostic for multimodal posteriors. arXiv:2609.37352 (preprint).

See also

gibbs(), rhat(), ess(), geweke_z()