Checking Prior Sensitivity

Why the scale is the reader’s problem

An analysis with bgm() states its assumptions as prior distributions, so the natural question is whether the conclusions depend on them. The interaction (pairwise) prior distribution is the one to worry about first, because its scale encodes how large you expect edge weights to be, and that expectation feeds directly into every inclusion verdict.

The scale is interpretable. A pairwise interaction in bgm() is half the log odds ratio between adjacent response categories, so the default normal_prior(scale = 1) places about two thirds of its prior mass on odds ratios between roughly 1/7 and 7. A reader who thinks edges are typically weaker wants a narrower scale; a reader who expects strong dependencies wants a wider one. If the inclusion verdicts survive that disagreement, they belong to the data. If they move, they belong to the scale, and the analysis should say so.

library(bgms)
data = Wenchuan

The check

prior_sensitivity_check() is that comparison, automated. It recovers each edge’s continuous inclusion-Bayes-factor curve across a range of slab scales and classifies how each verdict behaves along it. No special fit is required: it works on any bgm() fit made with edge selection — in the simplest case prior_sensitivity_check(bgm(data)).

The run below is deliberately small: the first six Wenchuan variables, two chains, and refits long enough for the per-anchor convergence gate to clear on so short a fit. Everything about how to read the report carries over to a full analysis. It was run once when this vignette was written rather than on every build — the answer does not change, and an installing user should not pay for it — so the call is shown as you would type it and the report beneath it is the object that call returned.

fit = bgm(data[, 1:6], seed = 1234, chains = 2, cores = 2,
          display_progress = "none", verbose = FALSE)
ps = prior_sensitivity_check(fit, seed = 1234, iter = 3000, warmup = 1000,
                             cores = 2)
ps
#> Prior sensitivity check: are the edge verdicts robust to the slab scale?
#> Bayes-factor curve from 0.4x to 2.5x the chosen scale (anchors at 0.4x, 0.63x,
#> 1x, 1.6x, 2.5x; the 1x anchor is the original fit); 15 edges.
#> 
#> 13 of 15 verdicts hold across the whole 0.4x-2.5x range; the exceptions are
#> named below.
#> 
#>   robust (same verdict at every scale)     13
#>   changed, within run-to-run noise          0
#>   changed, beyond run-to-run noise          2
#> 
#> 2 edges' verdicts genuinely depend on the scale:
#>   edge            0.4x       0.63x      1x         1.6x       2.5x   
#>   dreams-physior  undecided  undecided  absence    absence    absence
#>   dreams-avoidth  undecided  undecided  undecided  undecided  absence
#> 
#> Verdict counts by scale (at the anchors):
#>             0.4x 0.63x 1x 1.6x 2.5x
#>   presence     9     9  9    9    9
#>   undecided    4     4  3    3    2
#>   absence      2     2  3    3    4
#> More absence at wider scales is expected: a wider slab strengthens evidence
#> against borderline edges.
#> 
#> Note: the chosen scale (1) is much wider than the estimated interactions (about
#> 0.201 [0.131, 0.309]); absence verdicts in particular depend on this choice.
#> 
#> Method:  43-point curve from 5 anchor fits (0.4x to 2.5x the chosen scale),
#>          joined by importance reweighting; the 1x anchor is the original fit.
#>          Points with reweighting effective sample size below 400 are not shown.
#> Refits:  5 nuts refits, warm-started from the original fit, 33 s total.
#> Noise:   two identical refits at 1.6x differed by up to 0.79 log BF across
#>          threshold-relevant edges; verdict moves smaller than that are reported
#>          as run-to-run noise, not prior sensitivity.
#>          1 edge saturated in one of the two and was left out of that spread.
#> See ?prior_sensitivity_check for the full construction.

The headline is the stability statement and the category counts: most verdicts are robust across the whole range, and the edges whose verdict genuinely depends on the scale are named, with their verdict at each anchor. The per-scale verdict table below them shows the expected drift: as the slab widens, evidence for inclusion weakens and edges move from undecided toward absence, since a wider slab penalizes the extra parameter an edge costs.

The anchored curve

The continuous curve is what Bartos et al. present, and it is the deliverable here. A single fit cannot produce it for a global scale — reweighting one fit across the whole range collapses — so the curve is anchored: the model is fit at a handful of fixed scales (the log-spaced anchors, default c(0.4, 0.63, 1, 1.6, 2.5) times the chosen scale) and the values between anchors are filled in by importance reweighting. Reweighting a fixed-scale fit to a nearby scale needs only the per-draw slab-density ratios over the included edges; the likelihood cancels, so no extra fit is needed between anchors, and there is no normalizing constant to integrate.

The chosen-scale (1x) anchor is the original fit itself — it is never refit — so the chosen-scale verdicts the check reports are exactly the analysis you already ran. Each other anchor is one short warm-started refit.

Reweighting is only trustworthy near an anchor; far from every anchor the importance weights concentrate on a few draws and the estimate is unreliable. The check measures this with the importance effective sample size at each display point and masks any point below ess_floor (default 400) to NA rather than reporting a number. Log-spaced anchors make neighbouring usable ranges overlap, so the masked-out regions close up into a continuous curve; if two anchors are too far apart to overlap, the check warns and suggests adding an anchor.

Run-to-run noise

Verdict-identity comparisons are flaky by construction: an inclusion Bayes factor that sits within Monte Carlo error of a threshold can change category between reruns of the same prior, with no real sensitivity present. To calibrate that, the check repeats one anchor refit and takes the spread of \(|\Delta \log \mathrm{BF}|\) between the two as the empirical noise band (its 95th percentile is the figure the report’s Noise: line quotes, computed over edges whose evidence is near the decision thresholds).

An edge is reported as “changed, beyond run-to-run noise” only if its verdict changes along the curve and its natural log Bayes-factor swing exceeds that band (and its own Monte Carlo error, and a small tolerance). An edge that flips inside the band is “changed, within run-to-run noise”: the check refuses to call it a move. This is what keeps a robustness tool from crying wolf on boundary edges. In the $edges$mover column these categories are stored as stable, indistinguishable-from-wobble, and moved-beyond-wobble.

plot(ps)

The plot’s title states the answer. Each edge’s evidence curve runs across the anchored range with dots at the anchor scales, the undecided band shaded and the verdict zones labeled; edges whose verdict depends on the scale are colored and labeled by name, everything else is muted background, and any low-ESS curve points leave visible gaps.

The per-edge table carries the chosen-scale verdict and natural log Bayes factor with its Monte Carlo standard error, the per-anchor verdict columns, the stability range, the mover category, and an insufficient flag for edges the fit cannot certify (its chains disagree, or the Bayes-factor band straddles a threshold); the report counts these as “not certifiable (too noisy to assess)” and names them, when there are any. The full display grid, the per-point importance ESS, and which anchor each point was read from are in ps$multipliers, ps$curve$ess, and ps$curve$anchor_used.

head(ps$edges)

The chosen scale versus the data

One line of the report costs no refit: how the chosen scale compares with the size of the estimated interactions (the root-mean-square of the included interactions in the original fit, with an approximate standard error). When the chosen scale is far wider or narrower than that — more than about a factor of two — the report says so, because inclusion verdicts near the boundary then depend on the choice, and the curve shows by how much.

Cost and samplers

The chosen-scale anchor is the original fit, so it costs nothing to reuse; only the other anchors are refit. With warm starts each refit begins from the original fit’s per-chain end state, and an ordinal fit sampled with NUTS also carries the adapted step size and mass matrix, so a short warmup suffices. The anchors are refit one after another, each running its chains in parallel across cores, and on Wenchuan the whole check costs about one original fit of wall-clock. A fit sampled with adaptive Metropolis or Gibbs refits with the same sampler; when that makes the check cost more than about a minute, it suggests refit_sampler = "nuts", which is sound because every comparison runs refit-against-refit — the original fit’s own draws are only ever reweighted, never compared against a differently-sampled refit. Continuous (GGM) and mixed-MRF fits refit cold.

Because warm starts place every chain near the chosen-scale posterior, between-chain dispersion is reduced by construction, so split-\(\hat R\) is a weakened diagnostic at the refits; the per-anchor convergence gate therefore also leans on the effective sample size and on per-chain verdict agreement. An anchor that fails its gate is dropped from the curve rather than folded in.