Per-group hierarchical anchors: Operational Guide

API, pre-flight (C7), coefficient summaries, and cross-group prediction (Block 6.5)

José Mauricio Gómez Julián

2026-07-06

1. What this vignette covers

Block 6.5 promotes the population-level reference \(\theta_{\text{ref}}\) of the AMM canonical decomposition

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

from a single global scalar (Block 6 default) to a per-group vector \(\theta_{\text{ref}}[g]\), \(g = 1, \ldots, J_{\text{groups}}\), sampled hierarchically from \(\mathrm{Normal}(\mu_{\theta_{\text{ref}}}, \sigma_{\theta_{\text{ref}}})\) with both hyperparameters estimated. The model under grouping is

\[\theta_i = \theta_{\text{ref}}[g_i] + a(x_i) + b(x_i) \odot \theta_{\text{ref}}[g_i] + W(\theta_{\text{ref}}[g_i])\,x_i, \qquad \theta_{\text{ref}}[g] \sim \mathrm{Normal}(\mu_{\theta_{\text{ref}}}, \sigma_{\theta_{\text{ref}}}),\]

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 user activates the grouping via the one-sided formula argument group of gdpar(). Default group = NULL reduces bit-exactly to the Block 6 scalar regime, so existing fits and tests do not change.

This vignette documents:

  1. The group API and what changes inside the model.
  2. The pre-flight condition (C7) that prevents structural aliasing of the per-group anchor with a or b.
  3. The coef() schema under grouping: per-group anchors plus the population-level hyperparameters.
  4. The predict() semantics for newdata that contains levels unseen at fit time.
  5. The relationship with the parametrization toggle of vop01_parametrization_toggle.

2. API

2.1. The group argument

group is a one-sided formula identifying a single variable in data:

fit <- gdpar(
  formula = y ~ x1 + x2,
  family  = gdpar_family("gaussian"),
  amm     = amm_spec(a = ~ x1 + x2),
  data    = df,
  group   = ~ school    # promotes theta_ref to per-school anchor
)

Constraints accepted by the resolver .resolve_group_argument():

2.2. What changes inside the model

3. Pre-flight condition (C7)

When group is supplied, the pre-flight check .check_group_aliasing_c7() runs in addition to the standard Block 1 identifiability check. It enforces that neither a nor b contains columns that are linearly dependent with the group indicator \(G\). The check has two layers:

  1. Constant-per-group detection. Any column of \(Z_a\) or \(Z_b\) that is constant within every level of group is rank-deficient with \(G\) and creates a perfect alias with \(\theta_{\text{ref}}[g]\). This catches the most frequent mistake: including factor(group) (or any deterministic function of group) in a or b.
  2. Joint QR rank check. For the centered design block \(Z\), form \([G \mid Z]\), normalize column-wise, and verify \(\mathrm{rank}([G \mid Z]) = \mathrm{ncol}([G \mid Z])\). This catches indirect aliases that the per-column check misses (a non-trivial linear combination of \(Z\) columns matching a column of \(G\)).

On violation, gdpar() aborts with gdpar_input_error naming the offending columns. The user either removes them from the AMM spec or drops the group argument.

Defense in depth: the standard Block 1 identifiability check (gdpar_check_identifiability()) frequently catches the direct factor(group) case before C7 even runs, because \(Z_a\) centered is rank-deficient with the implicit intercept that the model carries via \(\theta_{\text{ref}}\). C7 covers the residual subtler aliases.

See vignette("v01_amm_identifiability", package = "gdpar") Section 6.6.2 for the formal statement and proofs.

4. Coefficient summary under grouping

coef(fit) returns an object of class gdpar_coef with the following slots when use_groups = 1:

Under use_groups = 0 (default), the theta_ref data.frame keeps the Block 6 columns (k, mean, q05, q50, q95) without g, preserving backward compatibility bit-exactly.

as.data.frame(coef(fit)) returns a long-tidy table that always includes the g column (NA when grouping is inactive) plus rows tagged component = "mu_theta_ref" and component = "sigma_theta_ref" when present.

Example (untested):

cf <- coef(fit)
cf$theta_ref          # one row per (g, k)
cf$mu_theta_ref       # population-level mu
cf$sigma_theta_ref    # population-level sigma
as.data.frame(cf)     # long-tidy: component, g, k, identifier, x_name, mean, q05, q50, q95

The print method shows the population-level hyperparameter blocks first, followed by the per-group anchors and the per-component summaries.

5. Cross-group prediction

predict(fit, newdata = new_df, ...) under a grouped fit resolves new_df[[var_name]] against fit$group_info$levels to obtain the per-row group index:

The newdata data.frame must contain the same grouping variable. If absent, predict() aborts with gdpar_input_error.

The b-block of the predictor uses c_b directly under grouping (consistent with the Stan model that does not report b_coef per group). Mathematically the formula is the same as in Block 6 with Z_b %*% c_b replacing Z_b %*% b_coef * theta_ref, since c_b = theta_ref * b_coef by construction.

Example (untested):

new_df <- df[1:10, ]
new_df$school <- factor(c(rep("school_A", 8), rep("new_school_X", 2)),
                        levels = c(levels(df$school), "new_school_X"))
# predict() emits a single warning naming "new_school_X" and uses the
# prior predictive marginal for those two rows:
pred <- predict(fit, newdata = new_df, summary = "draws")

6. Interaction with the parametrization toggle

The parametrization toggle (parametrization, parametrization_a, parametrization_W) discussed in vignette("vop01_parametrization_toggle", package = "gdpar") is orthogonal to grouping: it controls the CP/NCP geometry of the hierarchical scales sigma_a and sigma_W on the basis coefficients, not the hierarchical scale sigma_theta_ref on the per-group anchors. The per-group anchor is always sampled from a standard random-intercept hierarchy without an additional CP/NCP toggle.

In practice, when a fit shows divergences under grouping, the diagnostic order is:

  1. Confirm pre-flight (C7) is satisfied (otherwise the divergences are structural, not geometric).
  2. Inspect fit$parametrization$meta for the per-component CP/NCP decision; switch to the alternative parametrization with the standard toggle if sigma_a or sigma_W is the culprit.
  3. If the divergences persist and concentrate on theta_ref[g], consider tightening gdpar_prior(sigma_theta_ref = ...), since with very few observations per group the random-intercept funnel can be severe regardless of the per-component parametrization.

7. Reproducibility checklist

8. See also

mirror server hosted at Truenetwork, Russian Federation.