2.1. Direct construction with dimwise()

When every coordinate of \(\theta_i\) shares the same additive and multiplicative bases (the most common case), wrap the common formulas with dimwise() and pass the result as dims:

library(gdpar)

spec_uniform <- amm_spec(
  p    = 2L,
  dims = dimwise(a = ~ x1 + x2, b = ~ x1),
  W    = W_basis(type = "polynomial", degree = 2)
)
print(spec_uniform)
#> <amm_spec> AMM Level 2
#>   p (dim theta_i)    : 2
#>   dims (per-k a, b)  :
#>     k = 1 : a = ~x1 + x2 ; b = ~x1
#>     k = 2 : a = ~x1 + x2 ; b = ~x1
#>   W (modulating)     : W_basis(type = 'polynomial')

The printer reports the AMM level, p, the per-coordinate a and b formulas (resolved from the uniform template), and the modulating basis. The internal slot spec_uniform$dims is a length-p list whose entries each carry an a and a b formula; bare formula values passed to dims are rejected to prevent silent recycling.

2.2. Per-dimension overrides with override()

When one or more coordinates need a different additive or multiplicative basis, layer override() on top of dimwise():

spec_mixed <- amm_spec(
  p    = 3L,
  dims = dimwise(a = ~ x1 + x2, b = ~ x1) |>
           override(k = 2L, a = ~ x1) |>
           override(k = 3L, b = NULL),
  W    = W_basis(type = "polynomial", degree = 2)
)
print(spec_mixed)
#> <amm_spec> AMM Level 2
#>   p (dim theta_i)    : 3
#>   dims (per-k a, b)  :
#>     k = 1 : a = ~x1 + x2 ; b = ~x1
#>     k = 2 : a = ~x1 ; b = ~x1
#>     k = 3 : a = ~x1 + x2 ; b = NULL
#>   W (modulating)     : W_basis(type = 'polynomial')

The semantics of “unchanged” versus “disabled” is encoded by missing() inside override(): omitting the argument inherits from the base, while passing NULL explicitly disables the component for that coordinate only. Multiple override() calls compose; calling override() twice with the same k replaces the previous override for that index.

2.3. Chainable construction with amm_build()

For programmatic construction (loops, generators, serialisation pipelines), the chainable amm_build() helper exposes the same semantics as dimwise() + override():

spec_chain <- amm_build(p = 3L) |>
  amm_set_a_uniform(~ x1 + x2) |>
  amm_set_a(k = 2L, ~ x1) |>
  amm_set_b_uniform(~ x1) |>
  amm_set_b(k = 3L, NULL) |>
  amm_set_W(W_basis(type = "polynomial", degree = 2)) |>
  amm_set_x_vars(c("x1", "x2")) |>
  as_amm_spec()
print(spec_chain)
#> <amm_spec> AMM Level 2
#>   p (dim theta_i)    : 3
#>   dims (per-k a, b)  :
#>     k = 1 : a = ~x1 + x2 ; b = ~x1
#>     k = 2 : a = ~x1 ; b = ~x1
#>     k = 3 : a = ~x1 + x2 ; b = NULL
#>   W (modulating)     : W_basis(type = 'polynomial')
#>   x_vars             : x1, x2

The chain is functionally equivalent to the direct construction of Section 2.2; as_amm_spec() finalises the builder by validating coherence and forwarding to amm_spec() on the appropriate path. Overrides accumulated via amm_set_a() / amm_set_b() survive subsequent calls to amm_set_a_uniform() / amm_set_b_uniform(), mirroring the composition behaviour of dimwise() |> override().

2.4. Canonical serialisation: amm_save_spec() and amm_load_spec()

The canonical file format is a small key: value grammar with a mandatory version header. Round-trip is bit-exact for polynomial and B-spline bases; user-defined W bases are rejected (their evaluator is an arbitrary R function that cannot be canonised). Parsing is purely lexical; no source() or eval() runs on the file contents, so loading from untrusted locations is safe.

tmp <- tempfile(fileext = ".gdpar")
amm_save_spec(spec_uniform, tmp)
cat(readLines(tmp), sep = "\n")
#> # gdpar_spec_version: 0.1.0
#> p: 2
#> a: NULL
#> b: NULL
#> x_vars: NULL
#> W.type: polynomial
#> W.degree: 2
#> dims.1.a: ~x1 + x2
#> dims.1.b: ~x1
#> dims.2.a: ~x1 + x2
#> dims.2.b: ~x1

spec_roundtrip <- amm_load_spec(tmp)
identical(spec_uniform$level, spec_roundtrip$level)
#> [1] TRUE

The serialised form records constructor inputs only, not derived state. In particular, a W_basis that was materialised at a specific \(p\) via W_basis(..., p = ...) is serialised as the equivalent unmaterialised basis; reconstruction at gdpar() time re-applies materialize_W_basis() with the p of the loaded spec.

3.1. Construction

Three input forms are accepted:

fam_mv1 <- gdpar_family_multi("gaussian", p = 2L)
print(fam_mv1)
#> <gdpar_family_multi> coord-wise factorization
#>   p              : 2
#>   homogeneous    : TRUE
#>   name           : gaussian
#>   link           : identity
#>   has_dispersion : TRUE
#>   did_status     : holds
#>   did_reference  : Block 1, Section 6.4 (Lemma 1B)
#>   param_specs    : mu (per_observation), sigma (population) per coord (x2)

fam_mv2 <- gdpar_family_multi(gdpar_family("poisson"), p = 3L)
print(fam_mv2)
#> <gdpar_family_multi> coord-wise factorization
#>   p              : 3
#>   homogeneous    : TRUE
#>   name           : poisson
#>   link           : log
#>   has_dispersion : FALSE
#>   did_status     : holds
#>   did_reference  : Block 1, Section 6.4 (Lemma 1B)
#>   param_specs    : mu (per_observation) per coord (x3)

fam_mv3 <- gdpar_family_multi(
  list(
    gdpar_family("gaussian"),
    gdpar_family("gaussian")
  ),
  p = 2L
)
print(fam_mv3)
#> <gdpar_family_multi> coord-wise factorization
#>   p              : 2
#>   homogeneous    : TRUE
#>   name           : gaussian
#>   link           : identity
#>   has_dispersion : TRUE
#>   did_status     : holds
#>   did_reference  : Block 1, Section 6.4 (Lemma 1B)
#>   param_specs    : mu (per_observation), sigma (population) per coord (x2)

The first form (name string) is the most common; the second wraps an existing gdpar_family object; the third supplies an explicit list, useful when you want to vary the link function across coordinates in a future heterogeneous-family release while keeping today’s homogeneous-family code working.

3.2. Auto-promotion

When gdpar() is called with amm$p > 1L and a univariate gdpar_family object, the function auto-promotes the family to gdpar_family_multi(family, p = amm$p) and emits an informational message (gdpar_family_promotion_message). The convenience saves typing in interactive use; production scripts that need silence can pass an explicit gdpar_family_multi or wrap the call in suppressMessages().

fit <- gdpar(
  formula = y ~ x1 + x2,
  family  = gdpar_family("gaussian"),  # auto-promoted to gdpar_family_multi
  amm     = amm_spec(p = 2L,
                     dims = dimwise(a = ~ x1 + x2)),
  data    = df
)
# Informational message emitted (class gdpar_family_promotion_message):
#   "Auto-promoted univariate family 'gaussian' to gdpar_family_multi
#    via gdpar_family_multi(family, p = 2) because amm$p > 1.
#    Pass an explicit gdpar_family_multi to silence."

3.3. Constraints

The current release enforces two coherence invariants that abort gdpar() before sampling:

• family$p must equal amm$p when a gdpar_family_multi is supplied explicitly.
• All per-coordinate families must share the same stan_id and link function (homogeneous case).

A gdpar_family_multi supplied to a spec with amm$p NULL or 1 is rejected with gdpar_input_error; the reverse mismatch (univariate family with amm$p > 1L) triggers the auto-promotion described in Section 3.2.

4.1. Materialisation at a given p

The polynomial and B-spline bases are separable: each coordinate of \(\theta_{\text{ref}}\) contributes an independent block of basis functions, concatenated in coordinate order. The total basis dimension and per-coordinate block indices are computed by materialize_W_basis(). This can run either eagerly at construction time (when p is supplied) or lazily inside gdpar():

wb_poly_p2 <- W_basis(type = "polynomial", degree = 2, p = 2L)
print(wb_poly_p2)
#> <W_basis>
#>   type   : polynomial
#>   degree : 2
#>   dim    : 4
#>   p      : 2 (multivariate, separable)
#>   block_indices: 1-2, 3-4
wb_poly_p2$block_indices
#> [[1]]
#> [1] 1 2
#> 
#> [[2]]
#> [1] 3 4

wb_bs <- W_basis(type = "bspline", degree = 3, df = 4, p = 2L)
print(wb_bs)
#> <W_basis>
#>   type   : bspline
#>   degree : 3
#>   df     : 4
#>   dim    : 8
#>   p      : 2 (multivariate, separable)
#>   block_indices: 1-4, 5-8
wb_bs$block_indices
#> [[1]]
#> [1] 1 2 3 4
#> 
#> [[2]]
#> [1] 5 6 7 8

The block_indices slot is a length-p list whose k-th entry holds the row indices of the basis output that correspond to coordinate k of \(\theta_{\text{ref}}\). For polynomial bases of degree degree, every block has length degree; for B-spline bases, every block has length df (or length(knots) + degree if knots is supplied instead of df).

4.2. Splitting a separable basis with as_per_k()

For introspection and per-coordinate diagnostics, as_per_k() returns a length-p list of univariate W_basis objects describing the contribution of each coordinate:

subs <- as_per_k(wb_poly_p2)
length(subs)
#> [1] 2
subs[[1L]]$dim
#> [1] 2
subs[[2L]]$dim
#> [1] 2

For user-defined W bases (type = "user"), as_per_k() returns NULL and emits a warning: the package cannot infer separability from an arbitrary R evaluator. Callers that need a per-coordinate decomposition for a user basis must construct it explicitly.

4.3. Why W stays top-level

Declaring W per coordinate would silently restrict the model class to the separable sub-class:

\[W(\theta_{\text{ref}}) x = \sum_{k=1}^{p} W_k(\theta_{\text{ref}}[k])\,x \quad \text{(separable)},\]

losing the cross-dimensional coupling that gives the AMM hierarchy its expressive power. The polynomial and B-spline bases are separable by construction, but they materialise as a single \(W\) object whose internal block structure is recorded in block_indices. Replacing this single object by a length-p list of independent bases is what the package refuses at construction time. The non-separable extension (cross-coupling polynomials and Gaussian-process priors on \(W\)) is on the roadmap; the current API is forward-compatible.

5.1. Outcome layout

The outcome must be a matrix column of data with ncol(y) == p. Two equivalent constructions:

df$y <- cbind(y1, y2)              # matrix column added to existing df
# or
df <- data.frame(x1 = ..., x2 = ...)
df$y <- y_matrix                   # y_matrix is a numeric matrix n x p

A vector outcome aborts with an informative gdpar_input_error. A matrix with ncol(y) != amm$p aborts as well. The package never imputes missing values: any NA in the outcome matrix aborts before sampling.

5.2. Parametrization aggregation

For p > 1, the pre-flight diagnostic is run per coordinate (Path B’ applied to each \(\theta_{\text{ref}}[k]\) independently). The per-coordinate decisions are aggregated to a per-component decision via the new argument parametrization_aggregation:

The argument is ignored for the univariate path (amm$p == 1L). The full per-coordinate report is stored at fit$parametrization$report as an object of class gdpar_preflight_report.

5.3. Full call signature for the multivariate path

fit <- gdpar(
  formula                     = y ~ x1 + x2,
  family                      = gdpar_family_multi("gaussian", p = 2L),
  amm                         = spec_uniform,
  data                        = df,                # df$y is n x 2 matrix column
  parametrization             = "auto",            # runs per-coord preflight
  parametrization_aggregation = "any_ncp",         # default; conservative
  iter_warmup                 = 1000L,
  iter_sampling               = 1000L,
  chains                      = 4L,
  seed                        = 42L
)

The user-facing arguments are the same as the scalar path; the multivariate dispatch is internal and based on amm$p > 1L. All scalar-path arguments (anchor, skip_id_check, adapt_delta, max_treedepth, refresh, verbose, seed, parametrization_a, parametrization_W) are honored identically on the multivariate path.

7.1. predict() returns a 3-D array

For multivariate fits, predict(fit, ...) returns a three-dimensional array of shape (S, n, p) with dimnames list(NULL, row_names, paste0("dim_", seq_len(p))). Three type modes and three summary modes mirror the scalar path:

arr <- predict(fit)                                 # type = "theta_i", summary = "draws"
dim(arr)                                            # (S, n, p)

arr_resp <- predict(fit, type = "response")         # per-coord inverse link applied
arr_se   <- predict(fit, summary = "mean_se")       # list of p data.frames
arr_q    <- predict(fit, summary = "quantiles")     # list of p data.frames

arr_new  <- predict(fit, newdata = df_new)          # reconstruct from posterior draws

The reconstruction on new data mirrors exactly the formula encoded in the multivariate Stan template:

\[\eta_{i,k} = \theta_{\text{ref}, k} + Z_{a,k}[i, \cdot] \cdot a_{\text{coef}, k} + Z_{b,k}[i, \cdot] \cdot c_{b, k} + \sum_{j=1}^{W_{\text{per\_k\_dim}}} \bigl(\theta_{\text{ref}, k}^j - \theta_{\text{anchor}, k}^j\bigr)\,W_{\text{raw}}[r_{k,j}, \cdot]\,\sigma_W\,X[i, \cdot]^{\top},\]

with \(r_{k,j} = (k - 1) W_{\text{per\_k\_dim}} + j\) and the \(\sigma_W\) multiplier present only when the modulating component was sampled in the non-centered parametrization (cp_W = FALSE). The internal helper performs the per-coordinate centering via the means recorded at fit time.

7.2. coef() returns a unified gdpar_coef object

The S3 class gdpar_coef provides a single representation for both scalar and multivariate fits. The structure is:

gdpar_coef
  $p               : integer, the dimension of theta_i
  $summary_stats   : character vector, c("mean", "q05", "q50", "q95")
  $theta_ref       : data.frame; cols (k, mean, q05, q50, q95) with p
                     rows when grouping is inactive (J_groups == 1);
                     cols (g, k, mean, q05, q50, q95) with J_groups * p
                     rows when grouping is active
  $a               : NULL (component absent) or list of length p; each
                     entry NULL (coord inactive) or data.frame
                     (term, mean, q05, q50, q95)
  $b               : same conventions as $a
  $W               : NULL or list of length p; each entry NULL or
                     data.frame (basis_idx, x_name, mean, q05, q50, q95)
  $mu_theta_ref    : NULL when grouping is inactive; data.frame with
                     cols (k, mean, q05, q50, q95) and p rows under
                     grouping (per-coord posterior summary of the
                     hyper-mean shared across groups)
  $sigma_theta_ref : NULL when grouping is inactive; data.frame with
                     cols (k, mean, q05, q50, q95) and p rows under
                     grouping (per-coord posterior summary of the
                     hyper-scale)
  $J_groups        : integer, the number of grouping levels (1 when
                     grouping is inactive)
  $group_levels    : NULL when grouping is inactive; character vector
                     of length J_groups with the original group labels
                     preserved from the data (the integer code in
                     $theta_ref$g maps to group_levels[g])

For scalar fits (p = 1), the per-component slots are length-1 lists, not bare data.frames. This unification is deliberate: downstream code does not need to bifurcate on p == 1 versus p > 1.

Under grouping (J_groups > 1; see vignette("vop03_grouped_anchors", package = "gdpar") for the public API and the canonical use cases), the theta_ref data.frame gains an integer column g running over the grouping levels, and the hyper-parameter slots mu_theta_ref / sigma_theta_ref are populated with per-coord posterior summaries (NULL otherwise). The original group labels are preserved verbatim in group_levels, so group_levels[theta_ref$g] recovers them.

In the multivariate path, the W slot honors cp_W: if cp_W = TRUE, the centered parametrisation already absorbs the scale and the slot reports W_raw directly; if cp_W = FALSE, the non-centered draws are multiplied by sigma_W per sample before computing quantiles, so the reported coefficients are on the natural modulating scale.

Three S3 methods aid inspection:

cf <- coef(fit)
print(cf)                                       # level = "global", default
print(cf, level = "coord")                      # per-coord means only
print(cf, level = "full")                       # per-coord with full quantiles

summary(cf)                                     # aggregated counts + theta_ref mean
format(cf)                                      # one-line representation

7.3. as.data.frame() for tidy pipelines

The flattener returns a long-tidy data.frame with columns (component, k, identifier, x_name, mean, q05, q50, q95), ready for dplyr and ggplot2:

df_coef <- as.data.frame(coef(fit))
head(df_coef)

# Using dplyr with explicit namespace to avoid library() in vignettes:
df_coef |>
  dplyr::filter(component == "a") |>
  dplyr::group_by(k) |>
  dplyr::summarise(mean_abs = mean(abs(mean)),
                   max_q95  = max(q95))

The identifier column carries the term for the a and b slots, the basis_idx (formatted as a string) for W, and NA for theta_ref. The x_name column is NA everywhere except for W rows.

9.1. Reading the CSV directly

csv_path <- system.file(
  "benchmarks", "results", "cp_ncp_hit_rate_multi.csv",
  package = "gdpar"
)
results <- utils::read.csv(csv_path, stringsAsFactors = FALSE)
str(results)

The CSV is long-tidy with one row per (scenario, p, k, component) combination. The column regime_truth is NA for borderline scenarios where no ground-truth parametrisation is defined a priori; hit is NA for those rows. The columns n_div_pred and n_div_alt are the divergent counts of the predicted and contrastive (alternative) fits.

9.2. Hit-rate aggregation

The reporter aggregates hits per component and per scenario, restricting to rows with a defined truth:

sub <- results[!is.na(results$regime_truth) & !is.na(results$hit), ]
hits_comp <- aggregate(hit ~ component, data = sub,
                        FUN = function(x) mean(as.logical(x)))
hits_sc <- aggregate(hit ~ scenario, data = sub,
                      FUN = function(x) mean(as.logical(x)))
print(hits_comp)
print(hits_sc)

The package ships a baseline run where the hit rate is 18 / 18 = 1.00 over the truth-defined rows of the canonical eight scenarios. In the contrastive comparison, n_div_alt >= n_div_pred holds for every scenario, confirming empirically that the pre-flight chose the parametrisation with fewer divergences.

9.3. Faceted plot

When ggplot2 is available, the reporter writes cp_ncp_hit_rate_multi.png (facets scenario ~ component, fill by hit / miss / borderline) and cp_ncp_div_pred_vs_alt.png (dodge per scenario, predicted vs alternative divergences). The relevant chunk is:

# Using ggplot2 with explicit namespace to avoid library() in vignettes;
# the actual reporter at inst/benchmarks/scripts/report_hit_rate_multi.R
# follows the same convention.
results$hit_label <- ifelse(
  is.na(results$hit), "borderline",
  ifelse(as.logical(results$hit), "hit", "miss")
)
ggplot2::ggplot(
  results,
  ggplot2::aes(x = factor(k), fill = hit_label)
) +
  ggplot2::geom_bar(width = 0.7) +
  ggplot2::facet_grid(scenario ~ component, scales = "free_x",
                       space = "free_x") +
  ggplot2::scale_fill_manual(values = c(
    "hit"        = "#2c7bb6",
    "miss"       = "#d7191c",
    "borderline" = "#fdae61"
  )) +
  ggplot2::labs(x = "k (coordinate)", y = "count", fill = "verdict") +
  ggplot2::theme_minimal()

Users adapting the reporter to their own scenarios should preserve the long-tidy CSV schema (scenario, p, k, component, regime_truth, regime_pred, hit, decision_reason, n_div_pred, n_div_alt, ebfmi_min, t_attr, t_info_cp, t_info_ncp) so that the aggregation code keeps working unchanged.

The companion vignette vop03_regression_testing documents a complementary tool: the four-layer comparator gdpar_golden_compare() that locks the posterior of a reference fit so that future runs detect any regression in the sampling-side output. The reporter of this section focuses on the decision of the pre-flight; the comparator focuses on the realised draws after the long fit.

10.1. Default aggregation: per-subject

For p > 1, the observational unit is the row (subject) of the input data. Following the coord-wise factorisation \(p(y_i \mid \theta_i) = \prod_k D_k(y_{ik} \mid \theta_i[k])\), the per-subject log-likelihood is \(\log p(y_i \mid \theta_i) = \sum_k \log p(y_{ik} \mid \theta_i[k])\). This is the default aggregation = "subject" and it matches the convention used by brms multivariate fits with set_rescor(FALSE), so the resulting elpd_loo values are directly comparable to per-coordinate competitors aggregated identically.

fit <- gdpar(
  formula  = y ~ x1 + x2,
  family   = gdpar_family_multi("gaussian", p = 2L),
  amm      = amm_spec(p = 2L, dims = dimwise(a = ~ x1 + x2, b = NULL)),
  data     = train_df,
  refresh  = 0L, seed = 42L
)
lo <- gdpar_loo(fit)
print(lo)
# elpd_loo, se_elpd_loo, Pareto-k summary as in loo::loo

10.2. Diagnostic aggregation: per-cell

aggregation = "cell" treats each pair \((i, k)\) as an independent observation, yielding PSIS-LOO over \(n \cdot p\) cells. It is useful when Pareto-\(k\) mass concentrates in a specific coordinate (a marginally identified component for that dimension), but it conflates subject-level and coordinate-level cross-validation: do not report cell-aggregated elpd_loo as comparable to subject-aggregated values from other methods.

lo_cell <- gdpar_loo(fit, aggregation = "cell")
sum(lo_cell$diagnostics$pareto_k > 0.7) # per-cell concentration

10.3. Pareto-\(k\) caveats

Pareto-\(k\) values above 0.7 signal that the PSIS approximation is unreliable for the affected observations. The standard refinements are loo::loo_moment_match() (cheap re-weighting) and loo::reloo() (per-observation re-fit, expensive). Both accept the gdpar_loo() output and the corresponding fit$fit cmdstanr object directly.

10.4. Experimental status

gdpar_loo() is flagged with @keywords experimental. The aggregation rule is stable; the signature may gain additional arguments in future versions (for example integrand for non-pointwise predictive quantities).

11.1. Separable W only

The package-provided \(W\) bases (polynomial and B-spline) are separable in the per-coordinate sense described in Section 4.3. Non-separable bases that cross-couple coordinates of \(\theta_{\text{ref}}\) require either W_basis(type = "user", basis_fn = ...) (with the caveat that as_per_k() returns NULL for user bases) or the future non-separable extension on the package roadmap.

11.2. Per-coordinate heterogeneous families are a future block (per-slot heterogeneity is implemented as of 8.3.7)

The release enforces that all per-coordinate families of gdpar_family_multi() share the same stan_id and link function. The package distinguishes two orthogonal kinds of family heterogeneity:

• Per-coordinate heterogeneity (p-side; e.g., one Gaussian coord and one Poisson coord within a single gdpar_family_multi()) is scoped for a later sub-phase and is not currently implemented. The list constructor of gdpar_family_multi(list(...), p = p) accepts heterogeneous inputs at the API level for forward compatibility, but rejects them at the homogeneity check until the heterogeneous multivariate Stan template lands.
• Per-slot heterogeneity (K-side; the \(K \geq 2\) distributional-regression slots of a single coord, e.g., a Gaussian fit with K = 2 slots where slot 1 carries the mean and slot 2 carries the log-scale, possibly with different links and even different stan_ids such as lognormal_loc_scale accessed via gdpar_family_custom_K()) is implemented as of Sub-phase 8.3.7 via the named-list public API on the family argument of gdpar(). See vignette("vop04_amm_intermediate", package = "gdpar") §3 for the per-slot heterogeneous recipe and the Stan helper apply_inv_link_by_id(link_id, eta) (canonised in inst/stan/amm_distrib_K.stan:228) that dispatches the per-slot inverse link.

11.3. Single global sigma_W shared across blocks

The multivariate Stan template uses one sigma_W[1] shared across all \(W\) blocks. When the pre-flight per-coordinate decisions for \(W\) are heterogeneous, the package emits a gdpar_W_per_k_heterogeneous_message documenting that the sampler honors only the aggregated cp_W. The per-coordinate decisions remain in the report for auditability. Promotion of sigma_W to array[p] is scoped for Block 8 (multi-parametric extension) together with the related per-coordinate prior policy.

11.4. Multi-parametric extension is a future block

The factorisation \(p(y_i \mid \theta_i) = \prod_k D_k(y_{ik} \mid \theta_i[k])\) canonised here is the coord-wise option (Option B of Phase F). The alternative multi-parametric option (Option A: a single univariate outcome parametrised by the entire vector \(\theta_i\), e.g., Gaussian with \(\theta_i = (\mu_i, \log\sigma_i)\) in the distributional regression sense) is scoped for Block 8, after the coord-wise validation against TOP3 competitors (Blocks 6-7).

11.5. Pre-flight wall-time

The per-coordinate pre-flight adds roughly 30 % wall-time per gdpar() call (one compilation, one short fit with iter_warmup = iter_sampling = 200, two chains, adapt_delta = 0.95, max_treedepth = 10). In production pipelines where this cost is unacceptable, run parametrization = "auto" once during prototyping, read the resolved decisions from fit$parametrization, and pass them explicitly in subsequent calls via parametrization_a / parametrization_W.