Diagnostics for edge-selected parameters

With edge or difference selection active, the effect parameters are governed by spike-and-slab priors. The corresponding parameter is set to exactly zero when the effect is excluded, rather than being removed from the model. Because the parameter has a well-defined value at every iteration, mean, sd and Rhat are computed from the full chain, zeros included. The reported n_eff is deliberately not a full-chain ESS: run through the standard autocorrelation machinery, the spike-and-slab chain inflates the ESS in cells where the inclusion indicator, rather than the weight itself, drives the variability. The package therefore reports a composite estimator built from the two Monte Carlo errors that actually limit the model-averaged weight.

• n_eff in the pairwise table is the composite ESS of the model-averaged weight: it combines the Monte Carlo error of the conditional (included-only) weight with the Monte Carlo error of the inclusion probability.
• share_incl is the share of that composite Monte Carlo variance contributed by the inclusion part, so it says which of the two limits the precision of the weight.

The inclusion table carries its own diagnostics:

summary(fit)$indicator
#>                        mean         mcse           sd     n_eff      Rhat n0->1
#> intrusion-dreams  1.0000000           NA 0.000000e+00        NA        NA     0
#> intrusion-flash   1.0000000           NA 2.267338e-12        NA 0.9999995     0
#> intrusion-upset   0.9616211 1.354604e-02 1.875890e-01  191.7741 1.0091198    15
#> intrusion-physior 0.9597258 1.660955e-02 1.920383e-01  133.6780 1.0126830    11
#> dreams-flash      1.0000000           NA 0.000000e+00        NA        NA     0
#> dreams-upset      0.9874166 6.815380e-03 1.085553e-01  253.7007 1.0021993     6
#> dreams-physior    0.0997961 1.385844e-02 2.627090e-01  359.3527 1.0022622   110
#> flash-upset       0.1271870 1.344047e-02 3.016580e-01  503.7332 1.0000229   100
#> flash-physior     0.9999999 1.656812e-07 7.409489e-06 2000.0000 1.0000000     0
#> upset-physior     1.0000000           NA 0.000000e+00        NA        NA     0
#>                   n1->0
#> intrusion-dreams      0
#> intrusion-flash       0
#> intrusion-upset      15
#> intrusion-physior    11
#> dreams-flash          0
#> dreams-upset          5
#> dreams-physior      110
#> flash-upset         101
#> flash-physior         0
#> upset-physior         0

• mean, mcse, sd, n_eff, and Rhat describe the Rao-Blackwellized inclusion probability. The sampler records a one-step inclusion draw at every indicator update, and that draw is continuous, so the standard Monte Carlo machinery applies to it. They are NA when that draw carries too little variability for the machinery to work on, which happens by two routes. The first is an indicator whose Rao-Blackwellized draws are constant to double precision, where the inclusion probability sits at its numerical bound; sd is then exactly 0. The second is a draw whose variance is non-zero but falls below the numerical floor of the autocovariance kernel: mcse and n_eff come back NA while sd and Rhat are still reported. In the table above, intrusion-flash takes the second route — an sd of about 2e-11 with NA mcse and n_eff and a finite Rhat. Either way the reading is the same: the inclusion probability is pinned, not poorly estimated.
• n0->1 and n1->0 are the directional flip counts of the binary indicator chain. They record how much the chain explored the two states, which the Rao-Blackwellized columns cannot see: those measure precision conditional on exploration, so a chain stuck in one structure can show a large n_eff. Few flips are expected on an edge whose inclusion probability sits near 0 or 1, and the asymmetry between the two counts is itself informative.

Traceplots

Users can inspect traceplots by extracting raw samples directly. Here is an example for the pairwise effect parameter.

param_index = 1
chains = fit$raw_samples$pairwise
nchains = length(chains)
cols = c("firebrick", "steelblue", "darkgreen", "goldenrod")

plot(chains[[1]][, param_index],
  type = "l", col = cols[1],
  xlab = "Iteration", ylab = "Value",
  main = "Traceplot of pairwise[1]",
  ylim = range(sapply(chains, function(ch) range(ch[, param_index])))
)
if(nchains > 1) {
  for(c in 2:nchains) {
    lines(chains[[c]][, param_index], col = cols[c])
  }
}

E-BFMI

E-BFMI (Energy Bayesian Fraction of Missing Information) measures how efficiently the sampler explores the posterior. It compares the typical size of energy changes between successive samples to the overall spread of energies. Values close to 1 indicate that the sampler moves freely across the energy landscape; values below 0.3 suggest the sampler may be getting stuck or that the chain has not yet settled into its stationary distribution.

A low E-BFMI does not necessarily mean your results are wrong, but it does warrant further investigation. In models with edge selection, the most common cause is that the warmup period was too short for the discrete graph structure to equilibrate. Increasing warmup often resolves this.

Divergent transitions

Divergent transitions occur when the numerical integrator encounters regions of the posterior where the curvature changes too rapidly for the current step size. A small number of divergences (say, fewer than 0.1% of samples) is generally acceptable. However, many divergences indicate that the sampler may be missing important parts of the posterior.

If you see a large number of divergences, consider increasing target_accept (which makes the sampler use a smaller step size) and, if this does not fix it, switching to update_method = "adaptive-metropolis".

Tree depth

NUTS builds trajectories by repeatedly doubling their length until a “U-turn” criterion is satisfied. If the trajectory frequently reaches the maximum allowed depth (nuts_max_depth, default 10), it suggests the sampler may benefit from longer trajectories to explore the posterior efficiently. Hitting the maximum depth occasionally is normal; hitting it on most iterations may indicate challenging posterior geometry. If this happens, consider increasing nuts_max_depth.

Non-reversible steps

For MRFs with continuous variables, the leapfrog integrator enforces equality constraints through a projection step. After each forward step, the integrator checks whether reversing the step returns to the starting point. When the round-trip error exceeds a tolerance scaled by the square of the step size, the step is flagged as non-reversible.

A small number of non-reversible steps is not a concern. A large number indicates that the step size is too large for the constraint geometry. Because the step size is tuned during warmup, the most effective remedy is to increase warmup so the adapter has more time to find an appropriate step size. If non-reversible steps persist after increasing warmup, switch to update_method = "adaptive-metropolis".

Warmup and equilibration

Standard HMC/NUTS warmup is designed to tune the step size and mass matrix for the continuous parameters. In models with edge selection, the discrete graph structure may take longer to reach its stationary distribution than the continuous parameters. As a result, even after warmup completes, the first portion of the sampling phase may still show transient behavior (i.e., non-stationarity).

The warmup_check component provides simple diagnostics that compare the first and second halves of the post-warmup samples:

fit$nuts_diag$warmup_check
#> $warmup_incomplete
#> [1] FALSE FALSE
#> 
#> $energy_slope
#>      time_idx      time_idx 
#> -3.621392e-04 -8.963804e-05 
#> 
#> $slope_significant
#> time_idx time_idx 
#>    FALSE    FALSE 
#> 
#> $ebfmi_first_half
#> [1] 1.0201273 0.9853728
#> 
#> $ebfmi_second_half
#> [1] 0.8946329 1.0427978
#> 
#> $var_ratio
#> [1] 1.0088005 0.9924032

The returned list contains the following fields (one value per chain):

• warmup_incomplete: A logical flag that is TRUE when any of the indicators below suggest the chain may not have reached stationarity.
• energy_slope: The slope of a linear regression of energy against iteration number. A slope near zero indicates stable energy; a significant negative slope suggests the chain is still drifting toward higher-probability regions.
• slope_significant: TRUE if the energy slope is statistically significant (p < 0.01).
• ebfmi_first_half and ebfmi_second_half: E-BFMI computed separately for the first and second halves of the post-warmup samples. If the first-half value is much lower (for example, below 0.3) while the second-half value is healthy, the early samples were likely still settling.
• var_ratio: The ratio of energy variance in the first half to that in the second half. A ratio much greater than 1 (for example, above 2) indicates higher variability early on, consistent with transient behavior.

If these diagnostics suggest the chain was still settling, increase warmup and re-run the model. If diagnostics remain problematic after a substantial increase (for example, doubling or tripling warmup), consider re-fitting with update_method = "adaptive-metropolis" and comparing the posterior summaries. If the two samplers produce similar results, the estimates are likely trustworthy despite the warnings; if they differ substantially, that warrants further investigation of the model or data.

The two alarms

The gauge runs by default. It audits the chain’s decisions in a short set of assessment sweeps after sampling: on a subset of the edge moves the chain actually made, it redoes the add/remove decision the exact (slow) way and checks whether the outcome would have been different. Auditing the chain’s own moves is deliberate: it measures the approximation where the sampler is actually spending its time, rather than at arbitrary test configurations. The audit is produced under the hierarchical prior, and its summary is stored in fit$zratio_diag, with one row per chain. Turn it off with options(bgms.zratio_gauge_sweeps = 0L) before fitting, which leaves fit$zratio_diag as NULL.

Only edge moves whose mediating block holds two or more variables are audited. On a sparse posterior no block is that large, the normalizer ratio is exact, and the gauge reports nothing because there is nothing to audit: silence there means the approximation was never exercised, not that it went unchecked. Its cost follows the same rule. Two sweeps audit a capped number of moves regardless of how long the chain ran, so the gauge adds nothing to a sparse fit and about 5 to 10 seconds per chain to a dense fit at 100 to 200 variables — negligible against a production-length run, and worth switching off for short exploratory ones.

The gauge reports two complementary alarms. The first is flip_rate, the fraction of audited decisions that would have come out differently under the exact calculation. A chain is flagged on this channel above a 1% tolerance, after allowing for the small amount of Monte Carlo noise in the exact reference itself. Because the flip rate counts decisions rather than measuring a distance in abstract units, the 1% tolerance means the same thing regardless of the number of variables or the choice of prior.

The flip rate only registers error that lands near an accept/reject boundary, so a small but consistent error can pass it silently while still shifting the recovered inclusion probabilities: applied to every move a chain makes, a consistent lean relocates where the chain settles even when no single decision visibly changes. The second channel, harm_pred, projects the measured error onto that scale directly: it averages each audited edge’s signed error weighted by that edge’s own inclusion sensitivity (an error on an edge the data have pinned near 0 or 1 cannot move the result and does not count), then multiplies by the feedback amplification of the edge prior (learned inclusion probabilities amplify a consistent lean; a fixed inclusion probability does not). The result is the first-order predicted distortion of the mean posterior inclusion probability, flagged when it exceeds 0.01 and the weighted error is resolved above its own Monte Carlo uncertainty. It is a targeted alarm for a coherent lean, not a bound on all approximation error: edge-specific distortions that cancel in the mean are outside its scope. This channel is computed for the Bernoulli and Beta-Bernoulli edge priors; under other priors it is reported as NA.

The per_chain data frame in fit$zratio_diag records the following fields (one row per chain):

• flip_rate: The fraction of audited decisions that would have come out differently under the exact calculation. Together with harm_pred, one of the two fields to watch.
• flag: TRUE when the chain’s flip rate exceeds the 1% tolerance after allowing for reference noise.
• noise_floor: How much of the flip rate could be explained by Monte Carlo noise in the exact reference alone. The flag compares flip_rate against the tolerance plus this floor, so a chain is only flagged when the effect is larger than the reference’s own uncertainty.
• se_mean: The average difference between the approximate and exact edge calculations, measured on the log scale the sampler works on. A value near zero means the approximation is unbiased on this chain’s moves; a value consistently away from zero means it leans in one direction. An approximation can lean this way without changing any decisions when the data strongly favor one outcome; the harm channel exists to translate exactly that situation into consequences for the recovered graph.
• se_sd: The spread of that difference across the audited moves.
• se_mcse, se_se: The Monte Carlo uncertainty of se_mean: se_mcse is the contribution of the exact reference’s own noise, and se_se combines it with the spread across audited moves. The harm channel only flags a lean that is resolved above twice its own uncertainty (accounting for repeated audits of the same edge and for the reference noise), so a flag reflects a measured effect rather than noise.
• amplification, kappa, harm_pred, harm_flag: The harm channel. amplification is the feedback factor estimated from the chain’s own inclusion probabilities and the edge prior; kappa is the predicted shift of the mean inclusion probability per nat of coherent error (a regime summary: how exposed this fit is to approximation error, independent of how large the error turned out to be); harm_pred is the projected distortion itself; harm_flag is TRUE when that projection exceeds 0.01 and is resolved. A fit can be flagged here with a quiet flip rate: that combination means the error is not changing individual decisions but is consistently leaning the recovered graph in one direction.
• n_ent, n_ref, n_capped: How much the audit rests on. n_ent is the number of non-trivial edge moves the gauge saw, n_ref the number it actually checked against the exact calculation, and n_capped the number left unchecked because a per-sweep budget was reached. A small n_ref means both channels are based on few comparisons and should be read as a rough indication. On a dense large model the gap is wide by design — a few tens of checked moves against tens of thousands of non-trivial ones — so read the audit as a sample. A sample of that size resolves coherent error, error that leans the same way across edges, which is what harm_pred targets and what moves a recovered graph. It does not resolve rare edge-specific failures, and flip_rate is the flip rate among the audited decisions rather than a guarantee about every edge.
• block_lo, block_hi: The smallest and largest mediating block the audit covered, in variables. Read a flag against these: blocks beyond the approximation’s trained size range are where it is expected to drift.

Fits that skip the approximation instead

The approximation described above corrects each edge weight for the mediating structure around that edge: the other variables that connect its two endpoints. How much that correction matters depends on the precision diagonal’s prior. A gamma_prior() shape above 10 concentrates the diagonal enough that the precision matrix becomes diagonally dominant and the mediating correction all but disappears, and past that shape the approximation has not been scored.

bgms handles that case by not approximating. Above shape 10 it drops the mediating correction and uses the isolated-edge value, which is exact for an edge with no mediating structure at all, so the entire error is the mediation left out. Measured against the slow exact calculation at shapes 12, 15 and 20, that is at most 0.00028 nats — roughly a tenth of the 0.003 nats the approximation is held to inside its scored range. A fit that takes this route says so in a note and records the count as n_isolated in its per-chain counters.

The note carries one caveat worth reading. The mediating correction grows with the standardized diagonal rate (eta), and the bound above was measured only at rates up to 2. Above that rate the same route is used, and the note says plainly that the measured bound does not cover your fit. The trust gauge still runs, and on such a fit it measures that residual directly rather than relying on the bound, so it remains the thing to read.

What to do when a chain is flagged

A flagged chain tells you the approximation is affecting that chain’s recovered structure — a flag on the flip channel means individual edge decisions are changing; a harm_flag means a consistent lean is shifting the inclusion probabilities — so the structure should be treated with caution. The approximation is a smooth surface fitted once per analysis over a range of component sizes and densities; it drifts when a fit visits graphs outside that trained range, most often in large or unusually dense models. Work through the following in order.

First, read what was measured. The flagged bullet reports the error share and the mediating-block sizes the audit covered (block_lo and block_hi in per_chain), and n_ref against n_ent says how much of the chain the audit saw. A flag driven by blocks well beyond the anchored range is a different situation from one driven by ordinary-sized blocks. The surface is anchored to components of up to 80 variables, and beyond that it is continued along its own boundary slope; a fit that relies on the continuation says so in a separate note, with the retained share and the largest block involved.

Second, resolve whether the signal is real. Raise options(bgms.zratio_gauge_sweeps) above its default of 2 and re-fit. More sweeps audit more moves, which narrows the Monte Carlo uncertainty on both channels; a flag that survives a larger audit is a measured effect, and one that disappears was noise. Narrowing the prior scale toward its trained range is worth trying at this point too, since a very wide scale is the other common cause of drift.

Third, and only if a flag persists, consider the joint prior specification, precision_graph_prior = "joint", which evaluates the edge weights without this approximation. Note what that change does: the joint specification is a different model, not a more accurate version of this one. Its graph marginal is the edge prior reweighted by the per-graph normalizer rather than the edge prior itself, so the inferential target moves. bgms will never make that switch for you.

You can also keep the recovered structure and treat the flagged edges as provisional, confirming them another way, for example by re-fitting the neighborhood of those edges with a narrower prior and checking that they survive.

For evidence-free fits (prior-only sampling, or data too weak to pin the edges) under a Beta-Bernoulli edge prior, a harm_pred flag deserves extra weight: the learned inclusion probability amplifies a small per-edge lean, and a narrower scale may not remove the residual. A direct check is available: under either specification the average posterior inclusion probability of a prior-only fit must equal the edge prior’s mean, which for beta_bernoulli_prior(a, b) is a / (a + b).