# Arbitrary `p`: Operational Cookbook for Multivariate Fits

### Building, fitting, predicting, and reporting AMM models with `p > 1` (Path 1)

#### José Mauricio Gómez Julián

#### 2026-07-06

# 1. What this vignette covers

The canonical AMM decomposition fitted by Path 1 of gdpar generalises from a scalar individual parameter \(\theta_i \in \mathbb{R}\) to a vector \(\theta_i \in \mathbb{R}^p\):

\[\theta_i[k] = \theta_{\text{ref}}[k] + a_k(x_i) + b_k(x_i)\,\theta_{\text{ref}}[k] + \bigl(W_k(\theta_{\text{ref}}) - W_k(\theta_{\text{anchor}})\bigr)\,x_i, \qquad k = 1, \ldots, p,\]

stacked across \(k\), this is the vector form

\[\theta_i = \theta_{\text{ref}} + a(x_i) + b(x_i) \odot \theta_{\text{ref}} + \bigl(W(\theta_{\text{ref}}) - W(\theta_{\text{anchor}})\bigr)\,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 package factorises the likelihood across coordinates (architectural Option B canonised in Phase F of Block 5.2):

\[p(y_i \mid \theta_i) = \prod_{k=1}^{p} D_k\bigl(y_{ik} \mid \theta_i[k]\bigr),\]

with cross-dimensional coupling carried exclusively by the modulating component \(W(\theta_{\text{ref}})\). The additive component \(a_k\) and the multiplicative component \(b_k\) depend only on the covariates \(x_i\) and therefore factorise per coordinate; declaring them per coordinate is the user’s responsibility and is what the multivariate API of `amm_spec()` enforces.

This vignette is the operational cookbook for \(p > 1\) fits. It documents:

1. The four construction paths for a multivariate `amm_spec`: direct, with `dimwise()` + `override()`, with the chainable `amm_build()` helper, and the canonical serialised form via `amm_save_spec()` / `amm_load_spec()`.
2. The multivariate family `gdpar_family_multi()`, its auto-promotion semantics, and the homogeneous-family restriction of the current release.
3. The modulating component on a multivariate spec: `W_basis(..., p = ...)` materialisation, the structural asymmetry between `a` / `b` and `W`, and the separable polynomial / B-spline default.
4. The full call cycle of `gdpar()` for `p > 1`: outcome layout as a matrix column, `parametrization = "auto"` together with `parametrization_aggregation`, and the per-coordinate pre-flight report.
5. Extracting predictions and coefficients: `predict()` returning a three-dimensional array, `coef()` returning the unified `gdpar_coef` class, and `as.data.frame()` for downstream `dplyr` / `ggplot2` pipelines.
6. A runnable end-to-end smoke fit with `p = 2`.
7. A short sub-section showing how to read the pre-flight calibration report produced by `inst/benchmarks/scripts/report_hit_rate_multi.R` so that users can build their own scenario-driven calibration workflows.
8. Known limitations and pointers to the planned multi-parametric extension (Block 8).

For the underlying canonical form, identifiability, and the cross-dimension condition C4-bis, see `vignette("v01_amm_identifiability", package = "gdpar")`. For the scalar parametrization toggle and the three-filter pre-flight diagnostic, see `vignette("vop01_parametrization_toggle", package = "gdpar")`. The asymptotic theory of Path 1 lives in `vignette("v04_asymptotics_path1_bayesian", package = "gdpar")`.

# 2. Building a multivariate spec

A multivariate spec differs from the scalar spec in one structural axis only: the per-dimension components `a` and `b` are supplied via the `dims` argument of `amm_spec()`, not via the top-level `a` and `b` arguments. The package rejects mixing the two paths in either direction (`amm_spec(p = 2L, a = ~ x1)` and `amm_spec(p = 1L, dims = ...)` both abort with an informative error). The modulating component `W` remains a single top-level argument regardless of `p`, because it depends on the full \(\theta_{\text{ref}}\) vector and therefore couples the dimensions.

## 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. Multivariate families: `gdpar_family_multi()`

The multivariate family declares one univariate family per coordinate of \(\theta_i\). The current release of Path 1 restricts the per-coordinate families to a homogeneous set (all entries must share the same `stan_id`); heterogeneous families per coordinate are deferred to a later sub-phase.

## 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. The modulating component for `p > 1`

The asymmetry between per-dimension components (`a`, `b`) and the cross-dimension component (`W`) is what makes the canonical AMM form expressive beyond the separable sub-class. The package therefore reflects this asymmetry at the API level: `dims` for `a` and `b`, a single top-level `W` argument for \(W\).

## 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. Calling `gdpar()` for `p > 1`

## 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`:

Strategy | Behaviour | When to use
--- | --- | ---
`"any_ncp"` (default) | Component is CP only when every coordinate’s decision is CP. A single NCP flips the component-wide decision to NCP. Conservative. | Default. NCP geometry is monotonically safer in the worst case (more warmup steps in high-info, fewer divergences in low-info).
`"majority"` | Component is CP if the strict majority votes CP. Ties break toward NCP. | When per-coordinate decisions are expected to be near-uniform and you want the dominant signal.
`"per_k"` | No aggregation; each coordinate keeps its own CP/NCP at the Stan-template level via `segment()`-based priors (Phase H.2 of Block 5.2). | When coordinates have markedly different information regimes; needs the per-k Stan template wiring active.

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.

# 6. Inspecting the pre-flight report

The S3 class `gdpar_preflight_report` bundles per-coordinate, per-component CP/NCP decisions together with an aggregated per-component summary. Access it via the exported accessors `preflight_per_dim()` and `preflight_global_decision()`, or via the S3 methods `print()`, `summary()`, `as.data.frame()`, and `format()`.

    # After the fit:
    rep <- fit$parametrization$report   # gdpar_preflight_report
    print(rep)                           # default: level = "global"
    print(rep, level = "dim")            # per-coordinate detail
    print(rep, level = "both")           # global + per-coord

    preflight_per_dim(rep)               # tidy data.frame, one row per (component, dim)
    preflight_global_decision(rep)       # one row per component
    as.data.frame(rep)                   # same as preflight_per_dim()
    summary(rep)                         # named list with aggregates

The per-coordinate table has columns `(component, dim, decision, decision_reason, n_divergent, div_pct, ebfmi_min, t_attribution, t_info_cp, t_info_ncp)`. The `decision_reason` column carries the active-filter code (`filter_attribution`, `filter_ebfmi`, `filter_info_high`, `filter_info_low`, `filter_info_ambiguous_ncp`, `filter_info_undefined_ncp`, `absent_or_degenerate`). The reason codes are consistent with the univariate guide; see `vignette("vop01_parametrization_toggle", package = "gdpar")` Section 4 for the canonical table.

The global table has columns `(component, global_decision, agreement, method)`. The `agreement` column reports the share of effective per-coordinate decisions matching the aggregated `global_decision` (or the modal-decision frequency when `method = "per_k"`). An agreement of 1.0 flags `uniform`; values strictly above 0.5 flag `mixed`; values at or below 0.5 flag `split`.

When per-coordinate decisions for \(W\) are heterogeneous (mixed CP / NCP), an informational message of class `gdpar_W_per_k_heterogeneous_message` is emitted at fit time documenting that the current model uses a single global `sigma_W[1]` shared across blocks. The per-coordinate decisions are recorded in the report for auditability; the sampler honors only the aggregated `cp_W`. Block 8 (multi-parametric extension) may promote `sigma_W` to `array[p]` and honor per-coordinate \(W\) decisions.

# 7. Extracting predictions and coefficients

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

# 8. End-to-end worked example

A minimal Gaussian smoke fit with `p = 2`, two coordinates of \(\theta_i\), two covariates entering both the additive and the modulating components.

    library(gdpar)

    set.seed(42L)
    n <- 300L
    df <- data.frame(
      x1 = rnorm(n),
      x2 = rnorm(n)
    )
    # True theta_ref = c(0.5, -0.5), beta_a 2x2, sigma_y = 0.3 per coord
    true_theta_ref <- c(0.5, -0.5)
    true_beta_a <- matrix(c(0.8, -0.6,
                            0.4,  0.7), nrow = 2L, byrow = TRUE)
    y_mat <- matrix(NA_real_, nrow = n, ncol = 2L)
    for (k in seq_len(2L)) {
      y_mat[, k] <- true_theta_ref[k] +
                    true_beta_a[k, 1L] * df$x1 +
                    true_beta_a[k, 2L] * df$x2 +
                    rnorm(n, sd = 0.3)
    }
    df$y <- y_mat
    str(df)
    #> 'data.frame':    300 obs. of  3 variables:
    #>  $ x1: num  1.371 -0.565 0.363 0.633 0.404 ...
    #>  $ x2: num  -0.00462 0.76024 0.03899 0.73507 -0.14647 ...
    #>  $ y : num [1:300, 1:2] 1.525 -0.281 1.063 0.816 0.713 ...

The spec uses a uniform additive basis across both coordinates (each \(\theta_i[k]\) depends on `x1` and `x2`):

    spec <- amm_spec(
      p    = 2L,
      dims = dimwise(a = ~ x1 + x2)
    )
    print(spec)
    #> <amm_spec> AMM Level 1
    #>   p (dim theta_i)    : 2
    #>   dims (per-k a, b)  :
    #>     k = 1 : a = ~x1 + x2 ; b = NULL
    #>     k = 2 : a = ~x1 + x2 ; b = NULL
    #>   W (modulating)     : NULL

The fit with `parametrization = "auto"` runs the per-coordinate pre-flight, aggregates with the conservative default `"any_ncp"`, and proceeds to the long fit. Short iteration counts are used here for vignette responsiveness; production calibration should match the defaults of the calibration scripts (`iter_warmup = 1000`, `iter_sampling = 1000`, `chains = 2-4`).

    fit <- gdpar(
      formula                     = y ~ x1 + x2,
      family                      = gdpar_family("gaussian"),  # auto-promoted
      amm                         = spec,
      data                        = df,
      parametrization             = "auto",
      parametrization_aggregation = "any_ncp",
      iter_warmup                 = 300L,
      iter_sampling               = 300L,
      chains                      = 2L,
      refresh                     = 0L,
      verbose                     = FALSE,
      seed                        = 42L
    )
    print(fit)
    #> <gdpar_fit>
    #>   path                 : bayes
    #>   family               : gaussian (link = identity)
    #>   AMM Level            : 1
    #>   p (theta_ref dim)    : 2
    #>   anchor               : [0, 0]
    #>   observations         : 300
    #>   identifiability_pass : TRUE
    #>   converged            : FALSE
    #>   rhat_max             : 1.011
    #>   ess_bulk_min         : 326.2
    #>   divergent_count      : 0

The print method reports the AMM level, `p`, the anchor vector, observation count, identifiability pass, and the convergence verdict together with R-hat / ESS / divergent summaries.

    fit$parametrization$cp_a
    #> [1] TRUE
    fit$parametrization$cp_W
    #> [1] FALSE
    fit$parametrization$cp_a_per_k
    #> [1] TRUE TRUE
    fit$parametrization$cp_W_per_k
    #> [1] FALSE FALSE

The pre-flight report exposes the per-coordinate decisions and the active filter for each:

    rep <- fit$parametrization$report
    print(rep, level = "both")
    #> <gdpar_preflight_report>
    #>   p (theta_ref dim) : 2
    #>   aggregation       : any_ncp
    #>   components        : a, W
    #>
    #>   global decisions:
    #>     component  global_decision  agreement  flag
    #>     a          CP               1.00       uniform
    #>     W          absent            NA
    #>
    #>
    #>   per-coordinate decisions:
    #>     component  dim    decision  decision_reason       t_info_cp  t_info_ncp
    #>     a              1  CP        filter_info_high      41.5       57.2
    #>     a              2  CP        filter_info_high      30.8         44
    #>     W              1  absent    absent_or_degenerate    NA         NA
    #>     W              2  absent    absent_or_degenerate    NA         NA

Coefficient extraction and tidy flattening:

    cf <- coef(fit)
    print(cf, level = "coord")
    #> <gdpar_coef>
    #>   p                 : 2
    #>   summary_stats     : mean, q05, q50, q95
    #>   components active : a(2/2) b(0/2) W(0/2)
    #>
    #>   theta_ref:
    #>  k    mean     q05     q50     q95
    #>  1  0.4799  0.4544   0.48   0.5073
    #>  2 -0.5278 -0.5578 -0.5275 -0.4973
    #>
    #>   coord k = 1:
    #>     a:
    #>  term    mean
    #>    x1  0.7795
    #>    x2 -0.6181
    #>
    #>   coord k = 2:
    #>     a:
    #>  term   mean
    #>    x1 0.4029
    #>    x2  0.698

    df_coef <- as.data.frame(cf)
    head(df_coef, 8L)
    #>   component  g k identifier x_name       mean        q05        q50        q95
    #> 1 theta_ref NA 1       <NA>   <NA>  0.4799498  0.4543930  0.4799642  0.5073288
    #> 2 theta_ref NA 2       <NA>   <NA> -0.5277693 -0.5577847 -0.5275476 -0.4972924
    #> 3         a NA 1         x1   <NA>  0.7794671  0.7522783  0.7793699  0.8054150
    #> 4         a NA 1         x2   <NA> -0.6181115 -0.6449905 -0.6179875 -0.5899020
    #> 5         a NA 2         x1   <NA>  0.4029126  0.3701883  0.4030815  0.4350416
    #> 6         a NA 2         x2   <NA>  0.6980163  0.6689458  0.6979309  0.7261308

Prediction with `summary = "quantiles"` returns a list of `p` data.frames, one per coordinate:

    q_list <- predict(fit, type = "response", summary = "quantiles")
    length(q_list)
    #> [1] 2
    head(q_list$dim_1, 4L)
    #>          q05        q50        q95
    #> 1  1.5039983  1.5532166  1.5952615
    #> 2 -0.4689861 -0.4295067 -0.3930356
    #> 3  0.7110963  0.7387112  0.7666663
    #> 4  0.4807776  0.5195424  0.5549394
    head(q_list$dim_2, 4L)
    #>            q05         q50        q95
    #> 1 -0.001862122  0.04888494  0.1043911
    #> 2 -0.239240930 -0.19585643 -0.1544547
    #> 3 -0.358338547 -0.32636237 -0.2939195
    #> 4  0.223606238  0.26874689  0.3107521

Diagnostics access via the convenience function:

    diagnostics(fit)
    #> <gdpar_diagnostics>
    #>   converged           : FALSE
    #>   rhat_max            : 1.011
    #>   ess_bulk_min        : 326.2
    #>   ess_tail_min        : 221.9
    #>   divergent_count     : 0
    #>   treedepth_saturated : 0
    #>   efmi_min            : 0.8772

The fit can be archived together with the spec via `amm_save_spec()`; the spec is the only object needed to re-build the design at load time (the data and the seed reproduce the rest).

# 9. Reporting parametrization decisions across scenarios

For workflows that calibrate the pre-flight against custom scenarios, the package ships two scripts in `inst/benchmarks/`:

- `calibrate_cp_ncp_multi.R` runs `preflight_parametrization_multi()` over a battery of multivariate scenarios (the package ships eight canonical ones covering homogeneous high-info, low-info, mixed per-coordinate, W-only, both-active with overlap, borderline, and `p = 3` heterogeneous). For each scenario, the script also runs a contrastive fit with the alternative parametrisation to measure the divergence delta. Wall-time roughly 15-25 minutes with `cmdstanr` cached.
- `report_hit_rate_multi.R` reads the resulting CSV (`inst/benchmarks/results/cp_ncp_hit_rate_multi.csv`) and produces a hit-rate summary, a borderline breakdown, a deduplicated table of `n_div_pred` versus `n_div_alt` per scenario, and two `ggplot2` faceted plots written next to the CSV.

The reporter is intended as a template; copy and adapt it for domain-specific scenario batteries.

## 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. PSIS-LOO via `gdpar_loo()`

The Stan template emits the per-observation log-likelihood as a generated quantity (`log_lik[i, k]` for `p > 1`, `log_lik[i]` for `p = 1`). The helper `gdpar_loo()` wraps these draws into the `loo::loo()` workflow and returns the standard `psis_loo` object with `elpd_loo`, its standard error, and the Pareto-\(k\) diagnostics.

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

## 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_id`s 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`.

# 12. References and cross-references

- `vignette("v01_amm_identifiability", package = "gdpar")` — canonical form, identifiability conditions (C1)-(C4), the cross-dimension condition C4-bis for \(p > 1\).

- `vignette("vop01_parametrization_toggle", package = "gdpar")` — scalar parametrization toggle, three-filter pre-flight, reason codes.

- `vignette("vop03_regression_testing", package = "gdpar")` — the four-layer comparator `gdpar_golden_compare()` and the `gdpar_snapshot_fit()` API, both currently flagged as experimental.

- `vignette("v04_asymptotics_path1_bayesian", package = "gdpar")` — asymptotic theory of Path 1.

- `?gdpar`, `?amm_spec`, `?dimwise`, `?override`, `?amm_build`, `?W_basis`, `?as_per_k`, `?gdpar_family_multi`, `?gdpar_check_identifiability`, `?preflight_per_dim`, `?preflight_global_decision`, `?diagnostics`, `?amm_save_spec`, `?amm_load_spec`, `?gdpar_loo`.

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

- Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2017). Stan: A probabilistic programming language. Journal of Statistical Software, 76(1).

- Vehtari, A., Gelman, A., and Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27(5), 1413-1432.