Checks Internals
Two functions turn a fit into a judgment. verdicts() works on the edge indicators, the binary variables that record whether each edge is in the model. It reads each indicator’s inclusion Bayes factor, translates it into one of three verdicts, and says whether a rerun of the sampler could change that verdict. calibration_check() asks whether the model’s conditional predictions land where they claim to. A third piece of machinery sits behind the prior sensitivity check. This is the anchored curve, which takes a handful of fits at fixed prior scales and fills in a continuous curve between them, without refitting.
The three live in R/verdicts.R, R/calibration_check.R, and R/anchor_curve.R. For argument lists and return formats, see Checking Tools in the Reference. This page documents the estimators.
The three-way rule
The decision rule itself is four lines. Write \(\ell\) for an indicator’s inclusion Bayes factor on the natural log scale: the evidence for including that edge over excluding it. Write \(t\) for the evidence threshold. With \(\lambda = \log t\):
\[ \text{verdict} = \begin{cases} \text{presence} & \ell \ge \lambda \\ \text{absence} & \ell \le -\lambda \\ \text{undecided} & \text{otherwise.} \end{cases} \]
A missing Bayes factor gives a missing verdict, not “undecided”. The distinction matters for indicators that were never updated: for those there is no Bayes factor to read, and “undecided” would wrongly present that gap as a weighed verdict. The threshold is validated to be a single finite number greater than one. That guarantees \(\lambda > 0\), so the two boundaries \(\lambda\) and \(-\lambda\) are always on opposite sides of zero. The same rule serves verdicts() and the per-anchor verdicts of the sensitivity check.
The Bayes factor comes from extract_inclusion_bf(), read on the natural log scale. Working in logs is not cosmetic here. The log_bf value stays finite where bf saturates at zero or infinity, and the fragility machinery below needs a finite number, because it measures distances to the verdict boundaries.
verdicts() requires edge selection to have been on. Without selection every edge is in the model, so there is no inclusion Bayes factor to read. The function then stops with that explanation rather than returning an empty table.
Reading the indicators in the right order
The verdict table has one row per indicator, in the order the raw indicator draws lay them out. That order is not the same for every model. A GGM (Gaussian graphical model) or ordinal fit lays its indicators out as the row-major upper triangle of the variable order. A mixed fit lays them out by block: discrete by discrete first, then continuous by continuous, then the cross block. That block layout is a different permutation whenever the discrete and continuous columns interleave in the original data. indicator_pair_index() recovers the right permutation by filling the mixed layout with draw positions and reading the positions back off. The Bayes factors and inclusion probabilities are then looked up out of their symmetric matrices at the positions the draws actually used.
bgmCompare() has its own layout again. For each variable it places that variable’s main-effect difference on the diagonal, followed by its pairwise differences. compare_indicator_index() builds that layout.
The fragility flag
A verdict near a boundary can flip on a rerun purely through Monte Carlo noise. The fragility flag turns that fact into a per-edge statement, so the warning attaches to the specific edges at risk rather than living as a general caveat in the prose.
The operating point. In a calibration study the data-generating truth was known, so every verdict could be graded. The study covered 37,010 graded edge-fits across ordinal, binary, and Gaussian graphical models. All 66 verdict errors sat within 0.58 of a threshold on the log Bayes factor scale, and no edge further out was ever misclassified. Verdict error is therefore a boundary phenomenon, and the flag marks the boundary region: an edge is fragile when a verdict boundary lies within two standard errors of its estimated evidence.
Two standard errors, because neither is enough alone. The flag uses two different standard errors. Both are standard errors of the logit inclusion probability. That scale is the log inclusion Bayes factor scale shifted by the constant prior inclusion odds, so both standard errors are directly comparable with \(\ell\).
The two-state standard error, two_state_se_logit(), models each binary indicator chain as a first-order two-state Markov chain: a process that sits in one of the indicator’s two states, 0 or 1, and flips between them. The routine counts transitions within each chain only, never across the join between chains. It applies Jeffreys smoothing to the transition rates:
\[ \hat a = \frac{n_{01} + \tfrac12}{n_0 + 1}, \qquad \hat b = \frac{n_{10} + \tfrac12}{n_1 + 1}. \]
Here \(n_{01}\) counts the transitions from state 0 to state 1, \(n_{10}\) counts the transitions in the other direction, and \(n_0\) and \(n_1\) count the draws from which each kind of transition could start. The half counts are the smoothing, and the smoothing is what keeps the estimate defined at zero flips, where the unsmoothed rates would be \(0/0\). From \(\hat a\) and \(\hat b\) the routine forms the stationary probability \(p = \hat a / (\hat a + \hat b)\) and the two-state effective sample size. It returns \(1/\sqrt{\text{ESS}\,p(1-p)}\). On its own this standard error caught 74% of the study’s verdict errors.
The Rao-Blackwellized standard error, rb_se_logit(), starts from the Monte Carlo standard error of the one-step inclusion draws. Every model class that runs edge selection produces those draws, the ordinal, mixed, and comparison models as much as the GGM; the construction is described once, under GGM Internals, the page where the machinery is introduced. It moves that standard error to the logit scale by the delta method, which here means multiplying by \(1/(p(1-p))\). It is NA wherever the incoming Monte Carlo standard error is NA, and that happens where the draws are constant to double precision. On its own it caught 94%.
The union rule. boundary_distance() measures the gap from the evidence to the nearer verdict boundary, \(\min(|\ell - \lambda|, |\ell + \lambda|)\), and divides that gap by each standard error in turn. An edge is fragile when either distance falls below 2. A missing distance abstains rather than votes: it neither flags the edge itself nor prevents the other standard error from flagging it. In the study the union caught all 66 errors, at the cost of flagging 3.0% of correct verdicts. In absolute terms that cost is roughly seventeen flags on correct verdicts for every error caught, which is the deliberate price of catching them all: most flagged verdicts are correct ones standing near a boundary. The rule also transfers across model types, so the same operating point serves Gaussian graphical models as it serves ordinal and binary ones.
A fragile verdict is not a wrong verdict. It is a verdict the run was too short to settle, and the remedy is more iterations.
Where the flag is not calibrated. Every arm of the study fitted a single graphical model, not a group comparison. On bgmCompare() difference indicators the same arithmetic runs, and it still marks verdicts near a boundary. But no study has measured what share of difference-verdict errors it catches, or how many correct verdicts it rejects. The returned object therefore carries flag_validated = FALSE for those fits, and the print method says so in plain language. The print method adds a second caveat: difference verdicts are contingent on difference_scale, whose default calibration is under study.
verdicts() draws no random numbers. Given a fit and a threshold it is deterministic, so there is no seed argument: there is nothing for a seed to control.
The reporting cap
format_log_bf() caps what it prints at \(10^4\), and the cap applies to the log Bayes factor \(\ell\), not to the Bayes factor itself. Inside the cap it prints a rounded value; outside the cap it prints an inequality instead. The reason is honesty about precision. Past the cap the digits carry no information a longer run would reproduce, and a saturated Bayes factor has no value to print at all. A log Bayes factor that rounds to zero prints as 0, never -0.0. The sign of a rounded-away quantity is not something the run established, so printing it would report more than is known.
Calibration
calibration_check() asks a different question: not whether an edge is there, but whether the model’s predictions are honest. The unit of the check is the conditional distribution \(P(x_{ij} \mid \text{the rest of case } i)\), the model’s prediction for variable \(j\) of case \(i\) given that case’s other variables. This conditional distribution is the model’s own regression unit, so the check runs one variable at a time. It does not pool across variables, because pooling would let miscalibrations in opposite directions cancel.
Discrete variables
For every case and every category threshold, the model issues a cumulative probability: the predicted chance that the case falls at or below that category. The data then record whether the case actually did. Isotonic regression, run through stats::isoreg(), estimates the observed frequency as a monotone function of the predicted probability; this is the construction of the CORP reliability diagram (Dimitriadis et al., 2021). For a calibrated variable that estimated curve tracks the diagonal, where observed frequency equals predicted probability. The final cumulative probability is 1 for every case, so it says nothing about calibration and is dropped.
The predictions are the posterior-mean ones from predict().
The consistency band. The band answers one question: how far from the diagonal does a model that is calibrated by construction wander on data this size? It is built by simulation. For each of nrep replicates, pav_panels() draws one category per case from that case’s own predicted distribution, using the inverse of the cumulative distribution function. It then recomputes the isotonic curve for that replicate. Pointwise quantiles across the replicates form the band.
Resampling the category, rather than the threshold events, is the load-bearing choice. One case’s cumulative threshold events are nested by construction: if the case falls at or below one category, it falls at or below every later one. Resampling those events independently would ignore that nesting and produce a band that is too narrow, and a too-narrow band makes ordinary wander look like miscalibration.
Continuous variables
A continuous variable’s prediction is a density, so there is no category to have fallen in. The same question is therefore asked through the probability integral transform, following the visual predictive checking recommendations of Säilynoja et al. (2025). Write \(y_{ij}\) for the observed value of variable \(j\) in case \(i\), \(y_{i,-j}\) for that case’s other variables, and \(F_j\) for the model’s conditional predictive distribution function of variable \(j\). The transformed value \(u_{ij} = F_j(y_{ij} \mid y_{i,-j})\) is Uniform(0, 1) exactly when the conditional predictive distribution is right. The panel plots the empirical distribution function of the \(u\) values against the uniform diagonal.
conditional_pit() builds \(F_j\) as the predictive mixture over ndraws posterior draws,
\[ u_{ij} = \frac{1}{D} \sum_{d=1}^{D} \Phi\!\left( \frac{y_{ij} - m_{ij}^{(d)}}{s_{ij}^{(d)}} \right), \]
where \(D\) is the number of draws, \(\Phi\) is the standard normal distribution function, and \(m_{ij}^{(d)}\) and \(s_{ij}^{(d)}\) are the conditional mean and standard deviation for variable \(j\) of case \(i\) under posterior draw \(d\). To receive those per-draw conditional parameters, rather than their average, the function uses the internal return_draws path of the prediction helpers. Parameter uncertainty therefore sits inside the distribution the observation is transformed by. That placement matters. A single draw’s Gaussian would be too narrow, by exactly the amount the posterior is uncertain, and the uniform reference would then no longer be the right one.
Its band is not resampled from the model. Under the transform the null distribution is Uniform(0, 1) no matter what the conditional density was. uniform_ecdf_band() can therefore simulate the empirical distribution functions of \(n\) independent uniforms, with \(n\) the number of \(u\) values, and take pointwise quantiles of those. The band depends on \(n\) alone, and that is why one band serves every continuous variable in a fit.
The two panel kinds therefore condition on different things. The isotonic panel conditions on posterior-mean probabilities; the transform panel conditions on the full predictive distribution. Both live on the unit square against the diagonal, so a mixed fit produces one figure and one summary table. A kind column records which construction produced each row.
What is seeded
Both bands consume random numbers from R’s own generator, and seed is the only control over them. calibration_check() calls set.seed() when seed is supplied, and it does nothing otherwise. The selection of posterior draws for the predictive mixture also goes through R’s generator, so one seed fixes the whole check. Leave seed unset and the summary numbers move a little between calls, by the width of the band’s own Monte Carlo error.
Per group
On a bgmCompare() fit the check runs per group. A group’s cases are the ones its own parameters predict. Pooling the groups would therefore let a variable that is predicted too high in one group cancel against another group, exactly as pooling variables would. Every variable in a comparison fit is discrete, so every panel is isotonic. The ndraws argument has no effect here, because predict.bgmCompare() issues posterior-mean predictions only.
The row alignment is worth stating, because it is a real trap. The fit stores its cases sorted by group, so the group membership vector the check uses is in that internal, sorted order. Supplied newdata, by contrast, arrives in the order the data were originally given. The check permutes it by the same stable sort before use. That is why the method insists on one row per fitted case, in the original order, rather than accepting an arbitrary matrix.
What the check does not say
Evaluated on the fitted data the check is in-sample. Each observation helped shape the parameters it is now judged against, so the transform is mildly under-dispersed and the check is conservative. Held-out newdata removes that advantage, which makes the check strictly harder to pass.
More important, calibrated conditional predictions do not imply that the model reproduces the joint distribution. A model can predict each variable well from the others and still understate how strongly they depend on one another. That question needs a display built on simulate().
The anchored curve
R/anchor_curve.R is the estimator behind the prior sensitivity check’s continuous curve. The check refits the model at a handful of fixed slab scales: fixed values for the scale of the slab prior, the prior placed on the included edge parameters. These fits are the anchors. This file fills in everything between the anchors without refitting. Sensitivity and Refits documents the grid, the refits, and the gate that decides which anchors are allowed to contribute.
Reweighting
Take a fit at anchor scale \(s_a\) and ask what its posterior would have been at a nearby scale \(s\). The importance weight of a draw is the ratio of the two unnormalized posterior densities evaluated at that draw. In that ratio the likelihood appears as the same factor in numerator and denominator, since the data and the draw are the same, so it cancels; the spike terms, the prior factors of the excluded edges, do not depend on the slab scale and cancel the same way. And because an anchor is fitted at a fixed scale, there is no posterior over the scale itself to integrate. What is left is a per-draw ratio of slab densities. anchor_log_weights() forms the logarithm of that ratio over the included edges. Write \(w_t(s)\) for draw \(t\)’s weight at target scale \(s\), \(m_t\) for the number of included edges in draw \(t\), \(\gamma_{te}\) for edge \(e\)’s indicator in draw \(t\), and \(\theta_{te}\) for that edge’s parameter:
\[ \log w_t(s) = -m_t \log\frac{s}{s_a} - \frac{1}{2} \left( \frac{1}{s^2} - \frac{1}{s_a^2} \right) \sum_{e:\,\gamma_{te}=1} \theta_{te}^2 \]
for a Normal slab. A Cauchy slab keeps the same leading \(-m_t \log(s / s_a)\) term and replaces the second one with the matching sum of \(\log(1 + (\theta/s)^2)\) differences. Mapping a Bayes factor as a function of the prior without refitting for every candidate prior goes back at least to Sinharay & Stern (2002), who traced the point-mass-prior Bayes factor over a grid; the per-draw form used here is the posterior-density-ratio reweighting identity of Bartoš et al. (2026), applied locally around each anchor. The anchored construction that stitches the anchors together is specific to bgms.
anchor_draws() supplies \(\theta\) in the frame the slab prior applies to. For a GGM fit that means each precision draw is multiplied by \(-\tfrac12\) before it is used, since the slab prior sits on the partial associations rather than on the precision elements. For a bgmCompare() fit one indicator can gate several parameters, meaning one indicator switches several parameters on and off together. A pairwise difference indicator gates that pair’s difference in every contrast. A main-effect difference indicator gates the whole block of that variable’s threshold differences in every contrast. compare_anchor_draws() builds the map of which parameters each indicator owns. The weight then sums over the gated parameters, while the inclusion probability stays reported per indicator.
anchor_reweight() then does self-normalized importance sampling at each target scale: it reweights the anchor’s draws by \(w_t(s)\) and normalizes the weights to sum to one. It does this both pooled across chains and per chain. Alongside each estimate it records the importance effective sample size, \((\sum w)^2 / \sum w^2\), a measure of how many draws effectively carry the weighted estimate.
Stitching
assemble_curve() builds the displayed curve one grid point at a time. At each point it asks which anchors may contribute. An anchor qualifies when its importance ESS at that point clears ess_floor and when it passed its own convergence gate. Each qualifying anchor contributes its reweighted inclusion probability, weighted by inverse variance on that scale:
\[ w_a = \frac{\text{ESS}_a}{\max\!\left(p_a(1 - p_a),\, 10^{-6}\right)}. \]
Here \(\text{ESS}_a\) is anchor \(a\)’s importance effective sample size at the grid point and \(p_a\) is its reweighted inclusion probability there. The \(10^{-6}\) clamp in the denominator keeps a saturated edge, one whose \(p_a\) sits at exactly 0 or 1, from dividing by zero. Pooling happens on the inclusion-probability scale, and the pooled value is transformed to the log Bayes factor afterwards. This order is what removes the staircase seams and the infinities that winner-take-all anchor selection produced at hand-off points: a capped-edge anchor no longer hands off discontinuously to a finite one.
Grid points where no anchor clears the floor are masked to NA rather than extrapolated. Each surviving point is tagged with its highest-ESS anchor for reporting. At an anchor’s own grid position, that tag is the anchor itself. The per-chain curves pool with the same weights as the pooled curve, so the downstream per-point Monte Carlo standard error and the chain-unanimity checks carry over unchanged.
The pooled rows are an approximation, and they are treated as one. Exactness lives on the anchor fits’ own Rao-Blackwellized statistics. The per-anchor verdict columns and every chosen-scale quantity are read straight from each fit, so the 1x column is exactly the original fit’s reported analysis. This is the multistate-bridge direction; fully self-consistent weights are future work.
See also
Sensitivity and Refits, Extractor Internals, Simulation and Prediction, Convergence Diagnostics, GGM Internals, Checking Tools reference, MCMC Diagnostics.