1. What this vignette covers

Path 1 of gdpar is a hierarchical Bayesian fit (Stan / HMC) of the AMM canonical decomposition

\[\theta_i = \theta_{\text{ref}} + a(x_i) + b(x_i) \odot \theta_{\text{ref}} + W(\theta_{\text{ref}})\,x_i,\]

where $\odot$ denotes the Hadamard (elementwise) product, coherent with the canonical notation of vignette("v00_framework_overview", package = "gdpar") §8.2 and vignette("v01_amm_identifiability", package = "gdpar") §3.3.

The additive component a and the modulating component W each carry a hierarchical scale parameter (sigma_a, sigma_W). For Stan to sample efficiently, each of these scales admits two parametrizations:

Non-centered (NCP). The hierarchical coefficient is written as raw * sigma, with raw ~ normal(0, 1). Recommended in the low-information regime, where the marginal posterior of the coefficient stays close to the prior and the funnel is broad.
Centered (CP). The hierarchical coefficient is sampled directly with coef ~ normal(0, sigma). Recommended in the high-information regime, where the marginal posterior is tightly identified and CP avoids the funnel at the bottom of the NCP geometry.

The function gdpar() exposes a flag parametrization (and per-component overrides parametrization_a, parametrization_W) that controls which parametrization is used for the long fit. Default value "auto" runs a short pre-flight NCP fit and decides per-component automatically. This vignette documents:

The API and when to use each mode.
The three filters that drive the automatic decision.
How to introspect the decision via fit$parametrization$meta.
A runnable end-to-end example.
Known limitations.

For the underlying canonical form and identifiability theory, see vignette("v01_amm_identifiability", package = "gdpar") and vignette("v04_asymptotics_path1_bayesian", package = "gdpar").

2. API

2.1. Global flag

fit <- gdpar(..., parametrization = "auto")  # default
fit <- gdpar(..., parametrization = "ncp")   # force NCP everywhere
fit <- gdpar(..., parametrization = "cp")    # force CP everywhere

"auto" runs the pre-flight diagnostic (a 200 warmup + 200 sampling, 2-chain NCP fit with adapt_delta = 0.95 and max_treedepth = 10) and applies three filters per component. Approximately 30% extra wall-time per gdpar() call.
"ncp" and "cp" skip the pre-flight and apply the choice to both a and W.

2.2. Per-component overrides

fit <- gdpar(..., parametrization_a = "cp",
                  parametrization_W = "ncp")

parametrization_a and parametrization_W accept "ncp", "cp", or NULL.
When NULL, the component inherits from the global parametrization.
When both per-component overrides are explicit, the pre-flight is skipped regardless of parametrization.
A mixed case parametrization = "auto", parametrization_a = "cp", parametrization_W = NULL runs the pre-flight only for W, forcing CP on a.

2.3. When to use which mode

Situation	Recommended mode
Default exploratory work, no prior knowledge of geometry	`parametrization = "auto"`
Calibration runs where reproducibility across configurations matters	Explicit `"ncp"` or `"cp"` to remove the data-dependent decision
Production pipelines where the pre-flight wall-time is unacceptable	Explicit choice, informed by an earlier `"auto"` run
Strong prior belief that the posterior is far from the prior (high `n`, low residual noise)	`"cp"`
Strong prior belief that the posterior is close to the prior (small `n`, weak signal)	`"ncp"`
`a` is high-info but `W` is low-info (or vice versa)	`parametrization_a` and `parametrization_W` separately

2.4. Anchor configuration

The anchor argument of gdpar() fixes the reference point against which the AMM deviation is measured. It admits three forms:

anchor = "prior_mean" (default): the anchor is set to the prior mean of $\theta_{\text{ref}}$. This is the standard choice for AMM identifiability under (C1)-(C7) and matches the canonical centering convention of the framework.
anchor = "empirical_y": the anchor is estimated by inverting the link function at the sample mean of $y$. Useful when no defensible prior on $\theta_{\text{ref}}$ is available or when an empirical reference point is preferred for interpretability.
anchor = <numeric>: an explicit numeric anchor. Accepts a scalar when $p = 1$, a length-$p$ vector when $p > 1$, and a length-$K$ vector when $K > 1$ (per-slot anchors). Validated by resolve_anchor() against the dimensions of the model.

The anchor interacts with the parametrization toggle: under CP the prior on $\theta_{\text{ref}}$ is centred on the anchor; under NCP the same anchor enters as the shift parameter in the non-centred reparametrization. The toggle is anchor-agnostic — both modes admit any of the three forms.

For the multi-slot case ($K > 1$), the anchor receives one value per slot; see vignettes vop05_distributional_K_dharma and vop07_eb_workflow for slot-specific recipes.

3. The “auto” decision: three filters

The pre-flight diagnostic of gdpar() is split into two paths: the multivariate path (R/preflight_multi.R) handles $p > 1$ and numbers its checks filters 1-3; the scalar path (R/preflight.R) handles $p = 1$ and numbers its checks filters 4-6 as the scalar homologues of filters 1-3 respectively (the numbering is historical and preserved for cross-reference with vop02 §6 where filters 1-3 are documented in their multivariate form). This vignette covers the scalar path; the three filters described below are the scalar homologues of the multivariate triple.

Orthogonal to the parametrization filters, the upstream basis-restricted identifiability check runs before the pre-flight fires. Its strictness is controlled by the gdpar() argument id_check_rigor = c("full", "fast") (default "full"), which propagates verbatim to gdpar_check_identifiability(rigor = ...). Under "full" the check aborts on any violation of C1-C4 (scalar path) or C4-bis (multivariate path); under "fast" the cross-coordinate C4-bis check is downgraded to a single consolidated warning at the end of the per-coordinate loop, letting the fit proceed. See vignette("v01_amm_identifiability", package = "gdpar") §6.6.1 for the conditions C1-C4 and §6.6.1.4 for the rigor = "full" versus rigor = "fast" contract.

When parametrization = "auto" (and at least one component has no explicit override), the pre-flight diagnostic samples a short NCP fit and applies three filters in priority order. The first filter to fire fixes the decision; subsequent filters only act on still-undecided components.

3.1. Filter 4 — Divergence attribution

Active only when the pre-flight produces at least one divergent transition.

For each component, every transition is scored by

\[g_t = -z_{\log \sigma}(t)\cdot z_{\|raw\|}(t),\]

where the z-scores are computed over all pre-flight transitions. High values of $g_t$ flag transitions concentrated in the funnel region of the NCP geometry. The test statistic is the centered mean over divergent transitions, divided by the standard error estimated from all transitions under the null that divergences are random with respect to the funnel geometry:

\[S = \frac{\bar g_{\text{div}} - \bar g_{\text{all}}}{\sqrt{\widehat{\mathrm{Var}}(g) / n_{\text{div}}}}.\]

Rejected at $\alpha = 0.025$ against the upper t-distribution with $n_{\text{div}}-1$ degrees of freedom — small $n_{\text{div}}$ self-regularizes via the heavier-tailed critical value. A rejection forces CP on that component with reason code filter_attribution.

3.2. Filter 5 — E-BFMI

Active when the minimum E-BFMI across chains is below 0.3.

Energy mixing pathologies in the pre-flight imply that the marginal of sigma is poorly explored regardless of which hierarchical coefficient drives the issue. The filter forces CP on any component not yet decided, with reason code filter_ebfmi. This filter cannot attribute the pathology per component because E-BFMI is a global chain-level diagnostic; it acts as a conservative blanket.

3.3. Filter 6 — Info ratio (Path B’)

Active on any component still undecided after filters 4 and 5.

This filter measures how much the data contracted the marginal posterior of the component’s effective coefficient relative to its prior scale. The implementation is Path B’, which differs from a naive per-coordinate t-test because of the rank deficiency of W evaluated at a scalar $\theta_{\text{ref}}$ (see Section 7).

3.3.1. Effective coefficients and reference scale

For component a, the effective coefficient at coordinate $j$ per draw $t$ is

\[\mathrm{eff}_a[t, j] = a_{\text{coef}}[t, j], \quad j = 1, \ldots, J_a^{\text{free}},\]

and the per-draw reference scale is $\mathrm{ref}_a[t] = \sigma_a[t]$.

For component W with a polynomial basis of degree $K$ and $d$ predictors, the effective coefficient at coordinate $j$ per draw $t$ is the linear contribution of the basis to the $j$-th column of X in $\eta$:

\[\mathrm{eff}_W[t, j] = \sum_{k=1}^K \bigl(\theta_{\text{ref}}[t]^k - \theta_{\text{anchor}}^k\bigr)\,W_{\text{raw}}[t, k, j]\,\sigma_W[t],\]

with the per-draw reference scale being the conditional prior standard deviation of $\mathrm{eff}_W[t, j]$ given the hyperparameters at draw $t$:

\[\mathrm{ref}_W[t] = \sigma_W[t]\,\sqrt{\sum_{k=1}^K \bigl(\theta_{\text{ref}}[t]^k - \theta_{\text{anchor}}^k\bigr)^2}.\]

3.3.2. Log info ratio and block bootstrap

Per coordinate $j$, the log info ratio on the full pre-flight sample is

\[\widehat{\ell}_j = \log\!\left(\frac{\overline{\mathrm{ref}}}{\widehat{\mathrm{sd}}(\mathrm{eff}_{\cdot, j})}\right),\]

and the statistic averaged across coordinates is $\widehat{m} = \overline{\widehat\ell}$. The standard error of $\widehat{m}$ is estimated by a chain-aware block bootstrap: within each chain, blocks of 10 contiguous draws are resampled with replacement, and the statistic is recomputed; the bootstrap standard deviation across $B = 200$ replicates gives $\widehat{\mathrm{SE}}(\widehat{m})$. The chain-aware design preserves the local MCMC autocorrelation that an iid bootstrap would erase, yielding a more honest (typically larger) standard error.

3.3.3. Z-tests and thresholds

Two asymptotic z-tests are applied against null hypotheses on the log info ratio:

\[z_{\text{cp}} = \frac{\widehat{m} - \log \tau_{\text{cp}}}{\widehat{\mathrm{SE}}(\widehat{m})}, \qquad z_{\text{ncp}} = \frac{\widehat{m} - \log \tau_{\text{ncp}}}{\widehat{\mathrm{SE}}(\widehat{m})}.\]

Decision rule (one-sided, $\alpha = 0.025$, critical value $z^* \approx 1.96$):

$z_{\text{cp}} > z^*$: declare CP, reason filter_info_high.
$z_{\text{ncp}} < -z^*$: declare NCP, reason filter_info_low.
Otherwise: declare NCP, reason filter_info_ambiguous_ncp (conservative).

Default thresholds (as of 2026-05-10) are $\tau_{\text{cp}} = 5$ and $\tau_{\text{ncp}} = 2$, calibrated against eight canonical scenarios in inst/benchmarks/calibrate_cp_ncp.R. The benchmark achieves overall hit rate 0.92 against the well-defined scenarios; the single residual miss is attributable to structural confounding (Section 7), not to a threshold defect.

4. Introspecting the decision

After fitting, the resolved parametrization and the pre-flight metadata are stored at fit$parametrization:

fit$parametrization$cp_a          # logical: was CP chosen for a?
fit$parametrization$cp_W          # logical: was CP chosen for W?
fit$parametrization$meta          # list with diagnostic statistics

The meta list contains:

Field	Meaning
`used_preflight`	Logical. `TRUE` when the pre-flight ran (i.e., at least one component needed a data-driven decision).
`n_divergent`	Integer. Divergent transitions in the pre-flight.
`div_pct`	Numeric. Proportion of divergent transitions.
`ebfmi_min`	Numeric. Minimum E-BFMI across pre-flight chains.
`t_attribution_a`, `t_attribution_W`	Numeric. Filter 4 t-statistic per component (`NA` if filter 4 did not run for that component).
`t_info_cp_a`, `t_info_ncp_a`	Numeric. Filter 6 z-statistics for component `a` against $\tau_{\text{cp}}$ and $\tau_{\text{ncp}}$.
`t_info_cp_W`, `t_info_ncp_W`	Numeric. Same for component `W`.
`decision_reason_a`, `decision_reason_W`	Character. Reason code per component (see below).

Reason codes:

Code	Meaning
`user_global`	The user passed `parametrization = "ncp"` or `"cp"`; pre-flight skipped.
`user_explicit_a` / `user_explicit_W`	The user passed an explicit per-component override.
`absent_or_degenerate`	The component is not declared in `amm_spec` or its dimension is degenerate.
`filter_attribution`	Filter 4 rejected; CP forced.
`filter_ebfmi`	Filter 5 fired; CP forced.
`filter_info_high`	Filter 6 rejected against $\tau_{\text{cp}}$; CP chosen.
`filter_info_low`	Filter 6 rejected against $\tau_{\text{ncp}}$; NCP chosen.
`filter_info_ambiguous_ncp`	Filter 6 did not reach either threshold; NCP chosen conservatively.
`filter_info_undefined_ncp`	Filter 6 returned `NA` (degenerate inputs); NCP chosen conservatively.

5. Worked example

A minimal Gaussian fit with parametrization = "auto". The data are simulated under a regime that should land on filter_info_high for a: large beta coefficients with low residual noise, n = 80.

library(gdpar)

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

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

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

Inspect the resolved parametrization:

fit$parametrization$cp_a
fit$parametrization$cp_W

Inspect the pre-flight metadata:

str(fit$parametrization$meta)

The relevant fields for a high-information regime are decision_reason_a (expected "filter_info_high") and the z-statistic t_info_cp_a (expected to exceed 1.96). In this configuration W is absent from the spec, so decision_reason_W is "absent_or_degenerate".

fit$parametrization$meta$decision_reason_a
fit$parametrization$meta$t_info_cp_a
fit$parametrization$meta$n_divergent
fit$parametrization$meta$ebfmi_min

6. Known limitations

6.1. Confounding between `a` and `W` when covariates overlap

When a and W share predictors in x_vars, the marginal posterior of the effective W coefficient is genuinely wider because data variability is split between the two components. Filter 6 reads this as “low information” and recommends NCP — which is statistically the correct choice given the marginal geometry, even when the true generating coefficients are large. The user-facing implication is that under confounding, decision_reason_W may be filter_info_low or filter_info_ambiguous_ncp even in regimes with strong signal. This is documented behavior, not a defect.

If you have prior knowledge that W is informative despite the confounding, override with parametrization_W = "cp" explicitly.

6.2. Identifiability of individual basis components of `W`

The Stan model evaluates W at the scalar theta_ref, so the data linearly identify only one combination per coordinate: $\sum_k (\theta_{\text{ref}}^k - \theta_{\text{anchor}}^k)\,W_{\text{raw}}[k, j]$. Individual W_raw[k, j] for different $k$ at the same $j$ are not separately identified beyond what the prior implies. Path B’ addresses this by computing the info ratio on the effective coefficient (the sum across $k$) rather than on individual W_raw[k, j], which is the correct identifiable quantity to measure.

6.3. Block bootstrap versus naive bootstrap

The chain-aware block bootstrap (block size 10) is conservative: it preserves local MCMC autocorrelation that an iid bootstrap would erase, yielding a typically larger standard error and a more demanding rejection criterion. Users who want to experiment with different block sizes can call gdpar:::preflight_info_ratio_t() directly with custom block_size and n_boot arguments, but this is an internal API not covered by stability guarantees.

6.4. Pre-flight wall-time cost

The pre-flight adds approximately 30% wall-time per gdpar() call (one compilation, one short fit with 200 warmup + 200 sampling, 2 chains, adapt_delta = 0.95, max_treedepth = 10). In production pipelines where this cost is unacceptable, run parametrization = "auto" once during prototyping, record the resolved decision, and pass it explicitly via parametrization_a / parametrization_W in subsequent calls.

7. Summary recommendations

Start with parametrization = "auto". It costs ~30% extra wall-time but removes a manual decision.
Inspect fit$parametrization$meta after the fit to learn how the decision was made; the reason code points directly to the active filter.
Override with explicit parametrization_a / parametrization_W when you have prior knowledge (especially under confounding), when calibrating across configurations, or when the pre-flight cost is unacceptable.
Default thresholds $\tau_{\text{cp}} = 5$, $\tau_{\text{ncp}} = 2$ are empirically calibrated against the eight canonical scenarios; in domain-specific applications you can tune them via inst/benchmarks/calibrate_cp_ncp.R.

References

vignette("v01_amm_identifiability", package = "gdpar") — canonical form and identifiability conditions.
vignette("v04_asymptotics_path1_bayesian", package = "gdpar") — asymptotic theory of Path 1.
Stan Development Team (2024). Stan User’s Guide, version 2.35. https://mc-stan.org/docs/stan-users-guide/
Betancourt, M. and Girolami, M. (2015). Hamiltonian Monte Carlo for hierarchical models. Current Trends in Bayesian Methodology with Applications, CRC Press.
Künsch, H. R. (1989). The jackknife and the bootstrap for general stationary observations. The Annals of Statistics, 17(3), 1217-1241. (Block bootstrap.)

Parametrization Toggle: Operational Guide