--- title: "**Arbitrary `p`: Operational Cookbook for Multivariate Fits**" subtitle: "Building, fitting, predicting, and reporting AMM models with `p > 1` (Path 1)" author: "**José Mauricio Gómez Julián**" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 3 vignette: > %\VignetteIndexEntry{Arbitrary p: Operational Cookbook for Multivariate Fits} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( echo = TRUE, message = FALSE, warning = FALSE, collapse = TRUE, comment = "#>" ) have_cmdstan <- requireNamespace("cmdstanr", quietly = TRUE) ``` --- # **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`: ```{r build-uniform, eval=TRUE} 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) ``` 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()`: ```{r build-override, eval=TRUE} 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) ``` 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()`: ```{r build-chainable, eval=TRUE} 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) ``` 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. ```{r serialize, eval=TRUE} tmp <- tempfile(fileext = ".gdpar") amm_save_spec(spec_uniform, tmp) cat(readLines(tmp), sep = "\n") spec_roundtrip <- amm_load_spec(tmp) identical(spec_uniform$level, spec_roundtrip$level) ``` 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: ```{r families, eval=TRUE} fam_mv1 <- gdpar_family_multi("gaussian", p = 2L) print(fam_mv1) fam_mv2 <- gdpar_family_multi(gdpar_family("poisson"), p = 3L) print(fam_mv2) fam_mv3 <- gdpar_family_multi( list( gdpar_family("gaussian"), gdpar_family("gaussian") ), p = 2L ) print(fam_mv3) ``` 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()`. ```{r promotion-illustrative, eval=FALSE} 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()`: ```{r W-materialize, eval=TRUE} wb_poly_p2 <- W_basis(type = "polynomial", degree = 2, p = 2L) print(wb_poly_p2) wb_poly_p2$block_indices wb_bs <- W_basis(type = "bspline", degree = 3, df = 4, p = 2L) print(wb_bs) wb_bs$block_indices ``` 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: ```{r W-per-k, eval=TRUE} subs <- as_per_k(wb_poly_p2) length(subs) subs[[1L]]$dim subs[[2L]]$dim ``` 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: ```{r outcome-illustrative, eval=FALSE} 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** ```{r gdpar-illustrative, eval=FALSE} 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()`. ```{r preflight-shape, eval=FALSE} # 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: ```{r predict-illustrative, eval=FALSE} 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: ```{r coef-illustrative, eval=FALSE} 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`: ```{r asdf-illustrative, eval=FALSE} 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. ```{r smoke-data, eval=have_cmdstan} 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) ``` The spec uses a uniform additive basis across both coordinates (each $\theta_i[k]$ depends on `x1` and `x2`): ```{r smoke-spec, eval=have_cmdstan} spec <- amm_spec( p = 2L, dims = dimwise(a = ~ x1 + x2) ) print(spec) ``` 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`). ```{r smoke-fit, eval=have_cmdstan} 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) ``` 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. ```{r smoke-params, eval=have_cmdstan} fit$parametrization$cp_a fit$parametrization$cp_W fit$parametrization$cp_a_per_k fit$parametrization$cp_W_per_k ``` The pre-flight report exposes the per-coordinate decisions and the active filter for each: ```{r smoke-preflight, eval=have_cmdstan} rep <- fit$parametrization$report print(rep, level = "both") ``` Coefficient extraction and tidy flattening: ```{r smoke-coef, eval=have_cmdstan} cf <- coef(fit) print(cf, level = "coord") df_coef <- as.data.frame(cf) head(df_coef, 8L) ``` Prediction with `summary = "quantiles"` returns a list of `p` data.frames, one per coordinate: ```{r smoke-predict, eval=have_cmdstan} q_list <- predict(fit, type = "response", summary = "quantiles") length(q_list) head(q_list$dim_1, 4L) head(q_list$dim_2, 4L) ``` Diagnostics access via the convenience function: ```{r smoke-diag, eval=have_cmdstan} diagnostics(fit) ``` 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** ```{r reporter-csv, eval=FALSE} 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: ```{r reporter-agg, eval=FALSE} 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: ```{r reporter-plot, eval=FALSE} # 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. ```{r psis-loo-subject, eval = FALSE} 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. ```{r psis-loo-cell, eval = FALSE} 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. - 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.