---
title: "Checking Prior Sensitivity"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Checking Prior Sensitivity}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: refs.bib
csl: apa.csl
link-citations: TRUE
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 6
)
```

# 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.

```{r, eval = FALSE}
library(bgms)
data = Wenchuan
```

```{r, include = FALSE}
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.

```{r, eval = FALSE}
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
```

```{r, echo = FALSE}
# The report below is the real object that call produced, run once with those
# exact arguments by data-raw/make-prior-sensitivity.R and stored
# beside this file. Nothing here is transcribed: the check is not re-run on
# every build, but what is printed is still `print()` on a genuine
# prior_sensitivity_check() result.
ps = readRDS("prior-sensitivity-ps.rds")
ps
```

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`.

```{r, eval = FALSE}
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`.

```{r, eval = FALSE}
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.
