1. What this vignette covers

This is an operational recipe. It walks you through gdpar_eb() end to end across the four path regimes canonized in Sub-phase 8.6 of gdpar:

It then shows how to compare an EB fit against a Fully-Bayes (FB) fit on the same data via gdpar_compare_eb_fb(), how to read the numerical diagnostics that surface when the marginal Hessian is poorly conditioned, and when to give up on EB and fall back to FB.

The theoretical canonization of EB-vs-FB — the asymptotic equivalence (Theorem 7A and its multivariate extension Theorem 7A*), the higher-order coverage discrepancy (Proposition 7B scalar / 7B* matricial / 7B* tensorial), the finite-sample compound decision bound (Theorem 7C / 7C* compound multi-slot), and the four discrepancy conditions of Proposition 7D — lives in the canonical vignettes vignette("v07_eb_vs_fb") (scalar) and vignette("v07b_eb_multivariate") (multivariate extension). Read those if you want to know why a given choice was made; read this one if you want to do it.

The chunks below default to eval = FALSE because they compile Stan models and take several minutes per fit; re-enable evaluation on a per-chunk basis or globally via knitr::opts_chunk$set(eval = TRUE) if you want to reproduce the runs locally.

2. Setup

library(gdpar)
# These two are Suggests; gdpar_eb() and gdpar_compare_eb_fb() require
# them at runtime.
library(cmdstanr)
library(posterior)

We will work with a synthetic dataset of size n = 150 on a single continuous outcome with one covariate, then enrich it to multivariate and multi-slot variants as we walk through the four regimes.

set.seed(20260526L)
n  <- 150L
df <- data.frame(x = stats::rnorm(n))
df$y_scalar <- 1.0 + 0.4 * df$x + stats::rnorm(n, sd = 0.3)
# Multivariate (p = 2) outcome for Path A.
df$y_p2 <- cbind(
  1.0 + 0.4 * df$x + stats::rnorm(n, sd = 0.3),
  -0.5 + 0.2 * df$x + stats::rnorm(n, sd = 0.4)
)
# Same dataset is fine for Path B (K > 1, p = 1) and Path C (K > 1,
# p > 1) by reusing y_scalar / y_p2 with a multi-slot family below.

3. Minimal gdpar_eb() call (K = 1 + p = 1)

The base regime mirrors the canonical gdpar() signature; the only new arguments are eb_correction = TRUE (default; applies the Proposition 7B scalar inflation to the conditional credible intervals) and laplace_control = list(...) (controls the multi-start Laplace maximizer of v07 Section 11.1 step (i)).

fit_eb <- gdpar_eb(
  formula        = y_scalar ~ x,
  family         = gdpar_family("gaussian"),
  amm            = amm_spec(a = ~ x),
  data           = df,
  iter_warmup    = 500L,
  iter_sampling  = 500L,
  chains         = 2L,
  refresh        = 0L,
  seed           = 1L,
  laplace_control = list(multi_start_M = 5L)
)
print(fit_eb)

The output reports the EB plug-in point estimate \(\widehat\theta_{\text{ref}}^{\text{EB}}\) (from the marginal Laplace), its marginal standard error from the Laplace covariance, the marginal Hessian condition number \(\kappa(H)\), the multi-start dispersion across the multi_start_M = 5 independent inits, and the conditional HMC convergence diagnostics (\(\widehat R\), ESS, divergences).

summary(fit_eb) returns a tidy table of EB credible intervals with the Proposition 7B scalar inflation applied:

summary(fit_eb)

4. Path A (K = 1, p > 1): multivariate outcome

For a \(p\)-dimensional outcome, the amm spec uses dimwise() (or a plain list of length \(p\)) to declare the per-coordinate components, and the family is promoted to a gdpar_family_multi of dimension \(p\). The Stan template pair amm_eb_marginal_multi.stan + amm_eb_conditional_multi.stan (canonized in Sub-phase 8.6.C under decision D34) is dispatched automatically:

fit_eb_A <- gdpar_eb(
  formula        = y_p2 ~ x,
  family         = gdpar_family_multi("gaussian", p = 2L),
  amm            = amm_spec(p = 2L, dims = dimwise(a = ~ x)),
  data           = df,
  iter_warmup    = 500L,
  iter_sampling  = 500L,
  chains         = 2L,
  refresh        = 0L,
  seed           = 2L
)
print(fit_eb_A)

The corresponding correction is matricial, \(C^*_{g,\alpha} \in \mathbb{R}^{p\times p}\) (Proposition 7B* of v07b Section 5.1). It reduces algebraically to the scalar Proposition 7B at \(p = 1\), so the upgrade from base to Path A is transparent.

5. Path B (K > 1, p = 1): multi-parametric family

For a multi-slot distributional regression (e.g. modelling both mu and sigma of a Gaussian K=2, or mu and phi of a Negative Binomial K=2), the amm input is either a named list of amm_spec (one per slot) or a gdpar_formula_set via gdpar_bf(...). The Stan template pair amm_eb_marginal_K.stan + amm_eb_conditional_K.stan is dispatched automatically:

fs <- gdpar_bf(y_scalar ~ a(x), sigma ~ a(x))
fit_eb_B <- gdpar_eb(
  formula        = fs,
  family         = gdpar_family("gaussian"),
  data           = df,
  iter_warmup    = 500L,
  iter_sampling  = 500L,
  chains         = 2L,
  refresh        = 0L,
  seed           = 3L,
  skip_id_check  = TRUE
)
print(fit_eb_B)

Coverage of stan_ids per Sub-phase 8.6.C decision D33 (relaxed): {1, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13} — Gaussian, Negative Binomial, Beta, Gamma, Lognormal-loc-scale, Student-t, Tweedie, ZIP, ZINB, Hurdle-Poisson, Hurdle-NB.

6. Path C (K > 1, p > 1): full K x p extension

Path C, canonized in Sub-phase 8.6.D under decision D36 = (alpha) + D37 = (i) + D38’’ = (h), composes Path A coordinate-wise factorization with Path B multi-parametric K-slot semantics: a single outcome matrix-column y[n, p] is shared across the K distributional slots, each carrying its own per-coordinate linear predictor.

The initial 8.6.D iteration restricts Path C to family$stan_id %in% c(1, 3) (Gaussian K=2, NB K=2) under decision D40’ to avoid the numerical caveat of Section 6.1 of the opening handoff (HMC condicional bajo plug-in EB cerca del borde de soporte logit/log links + warmup corto). Beta / Gamma / Lognormal / Student-t / Tweedie / mixtures are deferred to a later iteration.

y_p2_int <- matrix(
  rnbinom(n * 2L, size = 5, mu = exp(0.5 + 0.2 * df$x)),
  n, 2L
)
df$y_p2_int <- y_p2_int
fit_eb_C <- gdpar_eb(
  formula        = y_p2_int ~ x,
  family         = gdpar_family("neg_binomial_2"),
  amm            = list(
    mu  = amm_spec(p = 2L, dims = dimwise(a = ~ x)),
    phi = amm_spec(p = 2L, dims = dimwise(a = ~ x))
  ),
  data           = df,
  iter_warmup    = 500L,
  iter_sampling  = 500L,
  chains         = 2L,
  refresh        = 0L,
  seed           = 4L,
  skip_id_check  = TRUE
)
print(fit_eb_C)

The fit object carries new Path C-specific slots:

• theta_ref_kp_hat: 3D numeric array of shape [J_groups, K, p].
• theta_ref_kp_se: same shape; per-coordinate marginal standard errors derived from each slot’s Laplace covariance block.
• theta_ref_kp_cov_per_slot: named list of K matrices each p × p; the per-slot blocks of the joint Laplace covariance after the block-diagonal extraction of decision D43 = (a).
• correction_tensor_constant: 3D array of shape [K, p, p] with the Proposition 7B* tensor correction of decision D37 = (i).
• correction_tensor_dispositions: named character vector reporting whether each slot’s correction was applied ("ok"), "non_finite", "non_psd", "missing", or "disabled".

7. Numerical diagnostics: how to read them

Every gdpar_eb_fit carries a diagnostics_numerical slot with four entries derived from the multi-start Laplace strategy of Charter Section 2.8:

• kappa (base regime / Path A / Path B) or kappa_per_slot (Path C): the condition number of the marginal Hessian at the chosen MAP. The guard laplace_control$kappa_threshold (default 1e10) aborts the fit with gdpar_eb_numerical_error when exceeded. A kappa between 1e6 and 1e9 is a warning sign; consider tightening the prior on theta_ref or moving the anchor closer to the data via anchor = "empirical_y".

• lm_perturbation: the Levenberg-Marquardt ridge actually added to the marginal covariance when the bare Hessian was singular or non-PSD. A non-zero value means the Laplace approximation needed numerical stabilization; the result is still valid but the effective sample size for the EB anchor is smaller than the nominal Laplace draws would suggest.

• multi_start_dispersion: standard deviation of the log marginal across the multi_start_M independent inits, normalized by the absolute mean. A value above 0.05 triggers a diagnostic warning because it suggests multi-modality of the marginal likelihood (open question O5*-EBFB of v07b Section 9.5). Consider raising multi_start_M from the default 5 to 10–20 to stress-test the optimum.

• marginal_log_lik_history: the per-init log marginal achieved by each optimize() call. Useful for forensics when the dispersion is large.

8. EB vs FB comparison via gdpar_compare_eb_fb()

The companion function gdpar_compare_eb_fb() (canonized in Sub-phase 8.6.E) takes a gdpar_eb_fit and a gdpar_fit fitted on the same dataset and reports three operational diagnostics of the EB-vs-FB theory of v07:

fit_fb <- gdpar(
  formula        = y_scalar ~ x,
  family         = gdpar_family("gaussian"),
  amm            = amm_spec(a = ~ x),
  data           = df,
  iter_warmup    = 500L,
  iter_sampling  = 500L,
  chains         = 2L,
  refresh        = 0L,
  seed           = 1L
)
cmp <- gdpar_compare_eb_fb(fit_eb, fit_fb, level = 0.95,
                            tv_bins = 30L)
print(cmp)
summary(cmp)

The output carries three tables:

• theta_diff_table: per-anchor cell comparison of \(\widehat\theta_{\text{ref}}^{\text{EB}}\) and the FB posterior mean \(E_{\text{FB}}[\theta_{\text{ref}}]\), with the difference and the difference normalized by the FB standard error. Under the standing hypotheses of v07 Section 4 (EB-MARG-ID + PRIOR-FB-WEAK + HIER-COMPLEX), Theorem 7A predicts diff_rel close to zero up to \(O(n^{-1/2})\).

• tv_table: marginal empirical total variation distance per common \(\xi\) parameter, computed via histogram plug-in over the shared support. Under Theorem 7A, marginal TV \(\to 0\) in probability as \(n \to \infty\). A persistent large TV across many parameters suggests one of the discrepancy conditions of Proposition 7D (multi-modality of the marginal likelihood, near-singular Fisher information, informative prior, deep hierarchy).

• coverage_table: per anchor cell, the EB credible interval width (with inflation applied when eb_correction = TRUE), the FB credible interval width, and the ratio width_eb / width_fb. This operationally verifies the \(O(n^{-1})\) under-cover prediction of Proposition 7B (and its matricial / tensorial extensions). A ratio systematically below 1 confirms the EB under-cover; a ratio above 1 after correction suggests the inflation is over-correcting, which is acceptable in finite samples and consistent with the asymptotic guarantee.

9. Troubleshooting

10. Where to go next

• For the theoretical canonization of EB-vs-FB at the scalar level, read vignette("v07_eb_vs_fb").
• For the multivariate / multi-slot extension of the theory (Theorem 7A*, Proposition 7B*, Theorem 7C* compound multi-slot, open question O5*-EBFB), read vignette("v07b_eb_multivariate").
• For the AMM canonization on which both EB and FB operate, read the AMM family of canonical vignettes (v01–v06).
• For CATE / ITE applications and the EB-vs-FB shrinkage discussion in the context of heterogeneous treatment effects, read vignette("v08_amm_for_cate_ite_positioning") and the implementation vignettes v08b_cate_ite_bridge_implementation and v08c_meta_learner_comparison.
• If your data are dependent — temporally (a time series) or spatially — the EB posterior understates uncertainty because the likelihood assumes conditional independence. Measure the violation with gdpar_dependence_diagnostic() (lag-1 autocorrelation, Durbin-Watson, Ljung-Box) or its spatial sibling gdpar_spatial_dependence_diagnostic() (Moran’s I), and obtain dependence-robust standard errors and percentile intervals with gdpar_dependence_robust() (temporal block bootstrap) or gdpar_spatial_dependence_robust() (spatial block bootstrap), Block 9, Axis 2. This makes the uncertainty robust to dependence; gdpar does not model the dependence. The full recipe is in vignette("vop09_dependence_robust").

References

The canonical theoretical references for EB-vs-FB in the AMM are Petrone, Rousseau, and Scricciolo (2014) and Rousseau and Szabo (2017) for the asymptotic equivalence and merging results; Carlin and Gelfand (1990) for the higher-order coverage discrepancy and the inflation correction; and Robbins (1956) / Efron (2010) for the compound decision framework. Full bibliographic details are in vignette("v07_eb_vs_fb") Section “References Cited in This Block”.