2.1. Layer A — Structural (class signatures, shapes)

The structural snapshot records:

• fit_class, coef_class: the class chain of the fit and of coef(fit).
• p, theta_ref_nrow, theta_ref_cols: shape of the theta_ref summary.
• components, parametrization_keys: which top-level slots are present in coef() and which CP/NCP slots are present in fit$parametrization. Both are treated as sets, not as ordered vectors, so a refactor that reorders the slot names but preserves the set passes.
• a_class, a_length, b_class, b_length, W_class, W_length: the shape of each per-component slot.
• a_per_k_cols: the column names of each per-coordinate a data.frame, in coordinate order.

The comparator distinguishes ordered comparisons (identical()) from unordered ones (setequal()) per field. Ordered fields detect changes in semantically meaningful order (class inheritance chain, fixed-position summary columns); unordered fields tolerate refactors that permute slot names without changing the set. The structural layer therefore detects silent API changes (new slots, renamed columns, removed components) before any numerical comparison runs.

The keys of fit$parametrization are sorted with sort(method = "radix") to make the snapshot deterministic across locales. Without radix, testthat::test() (which sets LC_COLLATE = C) and a manual Rscript run with the system locale would produce different orderings for the same set, generating spurious diffs.

2.2. Layer B — Discrete (bit-exact integer diagnostics)

The discrete snapshot records integer-valued sampler diagnostics:

• n_divergent: total number of divergent transitions across chains.
• treedepth_max_n: number of transitions that hit max_treedepth.
• treedepth_max_value: maximum tree depth observed across draws.
• n_leapfrog_total_per_chain: per-chain total leapfrog evaluations.
• ebfmi_min: minimum E-BFMI across chains (numeric, but compared bit-exact via the layer because the diagnostic surface is discrete in the same logical sense).

The comparator enforces bit-exact agreement on every integer field. The motivation is that any change in these diagnostics indicates a change in the sampler trajectory: a new divergence is a regression, a different leapfrog count means the NUTS retuned differently. The bit-exact contract is tied to a specific cmdstan version (currently 2.38.0). When upgrading cmdstan and NUTS retuning shifts the leapfrog count for the same seed and data, the recommended migration is to move n_leapfrog_total_per_chain from bit-exact to relative tolerance (e.g., within 5 %), not to weaken the rule across the board.

2.3. Layer C — Continuous (Monte Carlo standard error tolerance)

The continuous snapshot records, per parameter element matching the prefix of theta_ref, a_coef, sigma_a, sigma_y, and lp__:

• mean, sd, ess_bulk, ess_tail, rhat: standard posterior summaries from posterior::summarise_draws().
• mc_se: Monte Carlo standard error of the mean, computed as sd / sqrt(ess_bulk).

The comparator tolerance for the mean is

\[|\bar\mu_{\text{obs}} - \bar\mu_{\text{exp}}| \le k_\sigma \cdot \mathrm{MC\_SE}_{\text{exp}},\]

with k_sigma = 3 as the default. The contract is principled: MC_SE is the asymptotic standard error of the posterior mean estimator under MCMC sampling noise, so \(k_\sigma \cdot \mathrm{MC\_SE}\) is the band within which two repeated fits of the same model are expected to agree under independent RNG seeds. The contract auto-adapts to the precision of the golden simulation: a golden with high ESS produces a narrow tolerance band; a golden with low ESS produces a wider one.

Compared to a fixed relative tolerance (e.g., “within 5 %”), the k_sigma * MC_SE contract has two structural advantages:

• It does not over-flag near-zero parameters where small absolute drifts blow up relative differences.
• It does not under-flag parameters with high uncertainty where a 5 % drift is well below the sampling noise floor.

Compared to bit-exact comparison, it does not over-flag floating-point noise from cmdstan upgrades or basic-arithmetic refactors that change accumulation order.

When k_sigma should be tightened (more sensitive detection): if the project does not tolerate posterior mean drift larger than one MC SE for some parameter, set k_sigma = 1.5 or k_sigma = 2. When k_sigma should be loosened: if the golden was produced with too-low ESS and the band is unrealistically narrow, increase k_sigma to 4 or 5, but consider regenerating the golden with longer chains as the first option.

2.4. Layer D — Sanity (absolute floors)

The sanity snapshot records aggregates that must hold regardless of the golden:

• rhat_max: maximum R-hat across all parameters.
• ess_bulk_min, ess_tail_min: minimum bulk and tail ESS.
• divergent_pct: divergent transitions as a fraction of total sampling iterations.
• ebfmi_min: minimum E-BFMI.

The comparator applies absolute thresholds:

sanity_floor = list(
  rhat_max      = 1.05,
  ess_bulk_min  = 100,
  ess_tail_min  = 100,
  divergent_pct = 0.01,
  ebfmi_min     = 0.3
)

Layer D is defence-in-depth: a fit that passes layers A-C against an obsolete golden but violates the absolute floors is still flagged. This catches the failure mode where a buggy refactor produces a fit that matches a buggy golden numerically but is itself diagnostically broken. The sanity_floor argument of gdpar_golden_compare() accepts a partial override that is merged onto the defaults via utils::modifyList().

2.5. The parametrization_resolved audit field

The fifth field of the snapshot records the resolved CP/NCP decisions of the fit:

• cp_a, cp_W: aggregated per-component flags.
• cp_a_per_k, cp_W_per_k: per-coordinate decisions for multivariate fits.
• aggregation: the strategy used ("any_ncp", "majority", "per_k").

This field is not directly compared by gdpar_golden_compare(); differences usually surface indirectly via the discrete or continuous layers (a different parametrisation typically produces a different leapfrog count, different MC SE, etc.). The field is recorded so that failure forensics can immediately distinguish “the model was sampled differently” from “the sampler regressed in the same parametrisation”.

3.1. Step 1 — Snapshot a reference fit

Fit the model under the configuration you want to lock (seed, iter counts, parametrisation), then call gdpar_snapshot_fit():

library(gdpar)

set.seed(42L)
n <- 200L
df <- data.frame(x1 = rnorm(n), x2 = rnorm(n))
df$y <- 1 + 0.6 * df$x1 - 0.4 * df$x2 + rnorm(n, sd = 0.3)

spec <- amm_spec(a = ~ x1 + x2)

fit <- gdpar(
  formula       = y ~ x1 + x2,
  family        = gdpar_family("gaussian"),
  amm           = spec,
  data          = df,
  iter_warmup   = 500L,
  iter_sampling = 500L,
  chains        = 2L,
  refresh       = 0L,
  verbose       = FALSE,
  seed          = 42L
)

snap_baseline <- gdpar_snapshot_fit(fit)
names(snap_baseline)
#> [1] "structural"               "discrete"                
#> [3] "continuous"               "sanity"                  
#> [5] "parametrization_resolved"

The snapshot is a plain list; printing it directly is verbose. Inspect specific layers:

snap_baseline$structural$p
#> [1] 1
snap_baseline$structural$components
#> [1] "theta_ref"       "a"               "b"               "W"              
#> [5] "mu_theta_ref"    "sigma_theta_ref" "J_groups"        "group_levels"
snap_baseline$discrete
#> $n_divergent
#> [1] 0
#> 
#> $treedepth_max_n
#> [1] 0
#> 
#> $n_leapfrog_total_per_chain
#> [1] 3382 3324
#> 
#> $treedepth_max_value
#> [1] 3
#> 
#> $ebfmi_min
#> [1] 0.9241451
snap_baseline$sanity
#> $rhat_max
#> [1] 1.016322
#> 
#> $ess_bulk_min
#> [1] 351.7967
#> 
#> $ess_tail_min
#> [1] 478.6572
#> 
#> $divergent_pct
#> [1] 0
#> 
#> $ebfmi_min
#> [1] 0.9241451
head(names(snap_baseline$continuous$theta_ref))
#> [1] "theta_ref[1]"
snap_baseline$continuous$theta_ref[[1L]]
#> $mean
#> [1] 0.9617158
#> 
#> $sd
#> [1] 0.02378407
#> 
#> $ess_bulk
#> [1] 1144.412
#> 
#> $ess_tail
#> [1] 600.6642
#> 
#> $rhat
#> [1] 1.003517
#> 
#> $mc_se
#> [1] 0.000703064

3.2. Step 2 — Persist the snapshot

Store the snapshot wrapped together with provenance metadata (cmdstan version, seed, iter counts, package version):

golden <- c(
  list(
    schema_version  = 1L,
    scenario_id     = "my_baseline",
    gdpar_version   = as.character(utils::packageVersion("gdpar")),
    cmdstan_version = cmdstanr::cmdstan_version(),
    R_version       = R.version.string,
    seed            = 42L,
    n               = nrow(df),
    p               = 1L,
    n_warmup        = 500L,
    n_sampling      = 500L,
    n_chains        = 2L,
    generated_at    = format(Sys.time(), "%Y-%m-%d %H:%M:%S %Z")
  ),
  snap_baseline
)
saveRDS(golden, file = "tests/regression/golden_my_baseline.rds")

The schema versioning field allows a future schema bump to detect obsolete goldens and refuse the comparison until they are regenerated. The provenance fields are not consumed by the comparator but make the failure forensics tractable when the golden is years old.

3.3. Step 3 — Compare a fresh fit against the golden

Refit with the same configuration (or under a new package version after a refactor) and compare:

golden <- readRDS("tests/regression/golden_my_baseline.rds")

fit_new <- gdpar(
  formula       = y ~ x1 + x2,
  family        = gdpar_family("gaussian"),
  amm           = spec,
  data          = df,
  iter_warmup   = golden$n_warmup,
  iter_sampling = golden$n_sampling,
  chains        = golden$n_chains,
  refresh       = 0L,
  verbose       = FALSE,
  seed          = golden$seed
)

snap_new <- gdpar_snapshot_fit(fit_new)
cmp <- gdpar_golden_compare(snap_new, golden, k_sigma = 3)
cmp$passed
cmp$by_layer
if (!cmp$passed) print(cmp$failures, row.names = FALSE)

The passed field is TRUE only when every layer is clean. The by_layer field reports the count of failures per layer (structural, discrete, continuous, sanity). The failures data.frame has one row per failed item with columns:

The compact string formatter (format_compact(), internal) is the same used by all four layer comparators, so failure rows are visually consistent across layers.

4.1. k_sigma: sensitivity of layer C

The default k_sigma = 3 corresponds to a one-sided \(z\) of 3, roughly the 99.87 % upper quantile of a standard normal. Tighten when the project requires sub-MC-SE drift detection:

cmp <- gdpar_golden_compare(snap_new, golden, k_sigma = 2)

Loosen when the golden was produced at too-low ESS:

cmp <- gdpar_golden_compare(snap_new, golden, k_sigma = 5)

The recommended default lower bound is k_sigma = 1.5; below that, layer C effectively requires bit-exact agreement on continuous summaries, which is brittle under floating-point arithmetic. The recommended default upper bound is k_sigma = 5; above that, layer C ceases to detect regressions that matter. If k_sigma = 5 is not enough, regenerate the golden with longer chains.

4.2. sanity_floor: domain-specific absolute thresholds

The default sanity floor is intentionally permissive: rhat_max <= 1.05 (rather than the stricter 1.01), ess_bulk_min >= 100 (rather than 400), etc. The intent is to catch broken fits, not to enforce production-quality convergence.

Tighten the floor for production goldens of long, well-converged fits:

cmp <- gdpar_golden_compare(
  snap_new, golden,
  sanity_floor = list(
    rhat_max = 1.01, ess_bulk_min = 400, ess_tail_min = 400,
    divergent_pct = 0.001, ebfmi_min = 0.5
  )
)

The override is merged onto the defaults via utils::modifyList(), so partial specifications work:

# Only tighten rhat, leave the rest at defaults.
cmp <- gdpar_golden_compare(
  snap_new, golden,
  sanity_floor = list(rhat_max = 1.02)
)

4.3. Combining k_sigma and sanity_floor

The two knobs are orthogonal. A common production configuration:

cmp <- gdpar_golden_compare(
  snap_new, golden,
  k_sigma = 2,
  sanity_floor = list(
    rhat_max = 1.01, ess_bulk_min = 400, ess_tail_min = 400,
    divergent_pct = 0, ebfmi_min = 0.5
  )
)

This refuses any divergence in the observed fit, demands \(R\!-\!\hat{}\) within 1 %, and tolerates only two MC SE of drift in posterior means.

7.1. Experimental status

Both gdpar_snapshot_fit() and gdpar_golden_compare() are flagged @keywords experimental in their Roxygen. The schema version is at 1L. Future Blocks 6-9 may add fields (e.g., per-coordinate posterior summaries for \(W\), leave-one-out information criteria, PSIS-LOO diagnostics) and bump the schema with a documented migration. Until the first stable release, treat the contract as version-tied and re-generate goldens after every minor package bump.

7.2. Stan implementation drift

The bit-exact comparison on n_leapfrog_total_per_chain (layer B) is tied to a specific cmdstan version. NUTS retuning across cmdstan versions can change the leapfrog count for the same seed and data without affecting the posterior. When upgrading cmdstan, regenerate the goldens; if regeneration is not viable, the recommended migration is to move n_leapfrog_total_per_chain from bit-exact (layer B) to relative tolerance (e.g., 5 %), not to loosen the rule globally.

7.3. RNG seed semantics across hardware

The package does not currently snapshot the RNG implementation. In principle, the same seed could produce different draws on platforms whose RNG implementations differ at the floating-point level. In practice, R 4.6 + cmdstan 2.38 are reproducible across the supported platforms (Linux x86_64, macOS arm64) for the seeds tested in the package’s own golden suite. Treat this as an empirical observation, not a contract.

7.4. Layer C operates on posterior means only

The continuous layer compares mean per parameter element. Posterior sd, ess_bulk, ess_tail, and rhat are recorded in the snapshot for forensic use but are not directly compared (the sanity layer aggregates rhat and ess into floor checks). Adding per-element comparisons on sd or quantiles is straightforward but currently not done; if your application needs it, fork the comparator and add a layer-C extension.

7.5. The structural layer is not exhaustive

Layer A compares a fixed set of structural fields chosen to detect the most common refactors (slot rename, column rename, class chain change, parametrisation slot rename). It does not introspect the contents of the parametrization_resolved field beyond recording the keys. New structural invariants that should be enforced after future API changes need to be added to build_structural_snapshot_golden() (internal helper) and to golden_compare_structural() together with a schema_version bump.