diff --git a/NEWS.md b/NEWS.md
index 476f2d60..0f668e1c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -27,7 +27,7 @@
 * Corrected inclusion-probability updates for `beta_bernoulli_prior()` on continuous (GGM) models: under the determinant-tilted precision prior, the conjugate Beta update omits a normalizing-constant factor and biases the sampled inclusion probability toward sparsity (at 5 variables with a uniform hyperprior its prior mean lands near 0.37 instead of 0.5). The update now draws from the corrected conditional using a table built from the prior distribution at the first fit of a model configuration and cached on disk (`tools::R_user_dir("bgms", "cache")`; a one-time cost of the order of minutes, announced when `verbose = TRUE`). The sampled inclusion probability is returned per chain in `fit$inclusion_parameter_samples`.
 * Corrected stochastic block model updates for `sbm_prior()` on continuous (GGM) models, using the same cached table: the block-probability draws, the block-allocation weights, and the new-cluster weight all carry the normalizing-constant correction. Without it the sampled partition collapses toward one block (at 20 variables the prior mean number of blocks lands near 1.0 instead of 1.87); with it a prior-only chain reproduces the model's partition prior (total variation 0.01-0.02 at 5-20 variables). `sample_ggm_prior()` accepts `edge_prior` objects (`bernoulli_prior()`, `beta_bernoulli_prior()`, `sbm_prior()`) and returns sampled allocations under the block model.
 * The normalizing-constant correction extends to mixed models with `beta_bernoulli_prior()` or `sbm_prior()`: the determinant tilt acts on the continuous precision block, so the correction table is built for the continuous variables and the block-structure corrections read continuous-continuous edges only. With fewer than two continuous variables no edge is tilted and the plain conjugate updates apply unchanged; with exactly two, the single tilted pair supports the inclusion-probability correction but not the block-model slope curve, so `sbm_prior()` warns and keeps the plain conjugate updates there. Prior-only mixed chains reproduce the Beta hyperprior on the inclusion probability and the partition prior on the number of blocks; without the correction the inclusion probability biases toward sparsity and the partition toward fewer blocks.
-* Hierarchical graph-prior specification for continuous (GGM) models: `bgm(graph_prior_spec = "hierarchical")` composes the edge prior and the precision prior as `p(Gamma) p(K | Gamma)` with `p(K | Gamma)` normalized per graph, so the graph marginal is exactly the edge prior (under the joint specification it is reweighted by the per-graph normalizer). Each between-edge move evaluates the normalizer ratio by a deterministic local Z-ratio approximation, calibrated online against a block-Gibbs oracle in an appended warm-up window (`calibration_window`) and frozen with a hull clamp before sampling. Requires `edge_selection = TRUE`, a `normal_prior()` interaction prior, and a shape-1 `gamma_prior()` (or `exponential_prior()`) precision scale prior. The same specification is available in `sample_ggm_prior(spec = "hierarchical")`. With-data simulation-based calibration at 50 variables, n = 25, Beta-Bernoulli(2, 4): posterior inclusion-probability bias +0.002 with uniform ranks (the additive approximation alone reads -0.03). Prior-only chains with a sampled inclusion probability at large dimension remain out of envelope (the alarm suite flags them). On mixed models the specification normalizes the continuous block `p(K_yy | Gamma_yy)`: the Z-ratio enters the continuous-continuous edge moves only (counts on the continuous subgraph), and at least two continuous variables are required.
+* Hierarchical graph-prior specification for continuous (GGM) models: `bgm(precision_graph_prior = "hierarchical")` composes the edge prior and the precision prior as `p(Gamma) p(K | Gamma)` with `p(K | Gamma)` normalized per graph, so the graph marginal is exactly the edge prior (under the joint specification it is reweighted by the per-graph normalizer). Each between-edge move evaluates the normalizer ratio by a deterministic local Z-ratio approximation, calibrated online against a block-Gibbs oracle in an appended warm-up window (`calibration_window`) and frozen with a hull clamp before sampling. Requires `edge_selection = TRUE`, a `normal_prior()` interaction prior, and a shape-1 `gamma_prior()` (or `exponential_prior()`) precision scale prior. The same specification is available in `sample_ggm_prior(spec = "hierarchical")`. With-data simulation-based calibration at 50 variables, n = 25, Beta-Bernoulli(2, 4): posterior inclusion-probability bias +0.002 with uniform ranks (the additive approximation alone reads -0.03). Prior-only chains with a sampled inclusion probability at large dimension remain out of envelope (the alarm suite flags them). On mixed models the specification normalizes the continuous block `p(K_yy | Gamma_yy)`: the Z-ratio enters the continuous-continuous edge moves only (counts on the continuous subgraph), and at least two continuous variables are required.
 * Z-ratio alarm suite: `summarize_zratio_diagnostics()` audits the frozen Z-ratio kernel on graphs each chain visited (targeted, random, and additive-zone channels against a measurement-only block-Gibbs oracle), checks the calibration stream for end-of-warmup drift, and reports regime context. The verdict compares the maximum audit error to a per-regime threshold and prints like other sampler warnings; the summary is attached as `fit$zratio_diag` (and `$zratio_diagnostics` on `sample_ggm_prior()` output).
 * The NUTS diagnostics summary now prints the `warmup_incomplete` flag (energy not stationary) it already computed.
 * `extract_prior_inclusion_probabilities()`: prior edge-inclusion probabilities in the same matrix layout as `extract_posterior_inclusion_probabilities()`, for prior/posterior inclusion-odds computations. Under the joint spike-and-slab prior on a continuous block the graph marginal is reweighted by the per-graph normalizer (positive-definite-cone mass shaped by the determinant tilt), so continuous-continuous edges do not keep the edge-prior marginal at any `delta` — with a uniform Beta-Bernoulli hyperprior at 3 variables the prior edge probability is about 0.37 rather than 0.5, and a fixed `bernoulli_prior(0.5)` at `delta = 0` lands near 0.27. The values are read from the cached correction table (`bernoulli_prior()`, `beta_bernoulli_prior()`) or estimated by a prior-only chain with the fit's own correction settings (`sbm_prior()`, cached on the fit). Mixed models report per-edge-class values (discrete-discrete, continuous-continuous, cross); ordinal models use the analytic edge-prior marginals, including the exchangeable partition mixture for `sbm_prior()`.
diff --git a/R/bgm.R b/R/bgm.R
index 862c12a5..92f5d6de 100644
--- a/R/bgm.R
+++ b/R/bgm.R
@@ -168,7 +168,7 @@
 #'   sampled inclusion probability is returned per chain in
 #'   \code{fit$inclusion_parameter_samples}.
 #'
-#' @param graph_prior_spec Character. How the precision prior composes with
+#' @param precision_graph_prior Character. How the precision prior composes with
 #'   the edge prior under edge selection for continuous (GGM) data:
 #'   \describe{
 #'     \item{"joint"}{(default) The un-normalised joint specification
@@ -189,10 +189,10 @@
 #'       sampling (\code{\link{summarize_zratio_diagnostics}}; the summary
 #'       is returned as \code{fit$zratio_diag} and issues print like other
 #'       sampler warnings). Requires \code{edge_selection = TRUE}, a
-#'       \code{normal_prior()} interaction prior, a shape-1
-#'       \code{gamma_prior()} (or \code{exponential_prior()}) precision
-#'       scale prior, and continuous data — either all-continuous (GGM) or
-#'       mixed with at least two continuous variables. On mixed data the
+#'       \code{normal_prior()} or \code{cauchy_prior()} interaction prior, a
+#'       shape-1 \code{gamma_prior()} (or \code{exponential_prior()})
+#'       precision scale prior, and continuous data — either all-continuous
+#'       (GGM) or mixed with at least two continuous variables. On mixed data the
 #'       normalizer lives on the continuous block \eqn{K_{yy}}, so the
 #'       Z-ratio enters the continuous-continuous edge moves only, with the
 #'       mediating-block counts read off the continuous subgraph; discrete
@@ -201,7 +201,7 @@
 #'   Default: \code{"joint"}.
 #'
 #' @param calibration_window Non-negative integer or \code{NULL} (default).
-#'   Only for \code{graph_prior_spec = "hierarchical"}: length of the
+#'   Only for \code{precision_graph_prior = "hierarchical"}: length of the
 #'   appended warm-up window in which the Z-ratio correction is calibrated
 #'   online against a block-Gibbs oracle and then frozen (with its hull
 #'   clamp) before sampling. \code{NULL} resolves to no window for
@@ -407,7 +407,7 @@ bgm = function(
   delta = NULL,
   edge_selection = TRUE,
   edge_prior = bernoulli_prior(0.5),
-  graph_prior_spec = c("joint", "hierarchical"),
+  precision_graph_prior = c("joint", "hierarchical"),
   calibration_window = NULL,
   na_action = c("listwise", "impute"),
   update_method = c("nuts", "adaptive-metropolis", "gibbs"),
@@ -577,7 +577,7 @@ bgm = function(
     delta = delta,
     edge_selection = edge_selection,
     edge_prior = edge_prior,
-    graph_prior_spec = graph_prior_spec,
+    precision_graph_prior = precision_graph_prior,
     calibration_window = calibration_window,
     update_method = update_method,
     target_accept = if(hasArg(target_accept)) target_accept else NULL,
diff --git a/R/bgm_spec.R b/R/bgm_spec.R
index 1b878eb4..13d65796 100644
--- a/R/bgm_spec.R
+++ b/R/bgm_spec.R
@@ -268,7 +268,7 @@ bgm_spec = function(x,
                     delta = NULL,
                     edge_selection = TRUE,
                     edge_prior = bernoulli_prior(0.5),
-                    graph_prior_spec = c("joint", "hierarchical"),
+                    precision_graph_prior = c("joint", "hierarchical"),
                     calibration_window = NULL,
                     # Legacy edge prior params (accepted for backward compat)
                     inclusion_probability = 0.5,
@@ -379,18 +379,18 @@ bgm_spec = function(x,
   # The Z-ratio constants are derived for the Normal slab with an exponential
   # (shape-1 Gamma) diagonal, on the continuous precision matrix, under edge
   # selection. Anything else keeps the joint specification.
-  graph_prior_spec = match.arg(graph_prior_spec)
-  if(graph_prior_spec == "hierarchical") {
+  precision_graph_prior = match.arg(precision_graph_prior)
+  if(precision_graph_prior == "hierarchical") {
     if(!model_type %in% c("ggm", "mixed_mrf")) {
       stop(
-        "graph_prior_spec = \"hierarchical\" needs a continuous precision ",
+        "precision_graph_prior = \"hierarchical\" needs a continuous precision ",
         "block to normalize; the current model_type is '", model_type,
         "'. Use the joint specification, or data with continuous variables."
       )
     }
     if(model_type == "mixed_mrf" && sum(variable_type == "continuous") < 2) {
       stop(
-        "graph_prior_spec = \"hierarchical\" on mixed data needs at least ",
+        "precision_graph_prior = \"hierarchical\" on mixed data needs at least ",
         "two continuous variables (the specification normalizes the ",
         "continuous-block prior across its graphs). Use the joint ",
         "specification."
@@ -398,27 +398,33 @@ bgm_spec = function(x,
     }
     if(!edge_selection) {
       stop(
-        "graph_prior_spec = \"hierarchical\" normalizes p(K | Gamma) ",
+        "precision_graph_prior = \"hierarchical\" normalizes p(K | Gamma) ",
         "across graphs and needs edge_selection = TRUE; with a fixed ",
         "graph the specifications coincide."
       )
     }
-    if(!identical(interaction_prior_type, "normal")) {
-      stop(
-        "graph_prior_spec = \"hierarchical\" requires a normal ",
-        "interaction (slab) prior; the Z-ratio constants are derived for ",
-        "the Normal slab. Use interaction_prior = normal_prior(), or the ",
-        "joint specification."
-      )
+    if(!interaction_prior_type %in% c("normal", "cauchy")) {
+      stop(sprintf(
+        paste0(
+          "precision_graph_prior = \"hierarchical\" supports a normal or Cauchy ",
+          "interaction (slab) prior. Got %s_prior(). Use interaction_prior = ",
+          "normal_prior() or cauchy_prior(), or keep precision_graph_prior = ",
+          "\"joint\"."
+        ),
+        interaction_prior_type
+      ))
     }
     if(abs(scale_shape - 1) > 1e-12) {
-      stop(
-        "graph_prior_spec = \"hierarchical\" requires shape = 1 on the ",
-        "precision scale prior (gamma_prior(shape = 1) or ",
-        "exponential_prior()); the Z-ratio constants are derived for the ",
-        "exponential diagonal. Adjust the prior, or use the joint ",
-        "specification."
-      )
+      stop(sprintf(
+        paste0(
+          "precision_graph_prior = \"hierarchical\" requires shape = 1 on the ",
+          "precision scale prior; the Z-ratio normalizer is derived for ",
+          "the exponential diagonal. Got shape = %s. Use ",
+          "gamma_prior(shape = 1) or exponential_prior(), or keep ",
+          "precision_graph_prior = \"joint\"."
+        ),
+        format(scale_shape)
+      ))
     }
   }
 
@@ -495,7 +501,7 @@ bgm_spec = function(x,
       scale_rate = scale_rate,
       scale_eta = scale_eta,
       delta = delta,
-      graph_prior_spec = graph_prior_spec,
+      precision_graph_prior = precision_graph_prior,
       calibration_window = calibration_window,
       edge_prior_flat = ep_flat
     )
@@ -522,7 +528,7 @@ bgm_spec = function(x,
       scale_rate = scale_rate,
       scale_eta = scale_eta,
       delta = delta,
-      graph_prior_spec = graph_prior_spec,
+      precision_graph_prior = precision_graph_prior,
       calibration_window = calibration_window,
       edge_prior_flat = ep_flat
     )
diff --git a/R/build_arguments.R b/R/build_arguments.R
index ba473000..2568f908 100644
--- a/R/build_arguments.R
+++ b/R/build_arguments.R
@@ -27,36 +27,35 @@ build_arguments = function(spec) {
 }
 
 
-
 build_arguments_ggm = function(spec) {
   list(
-    num_variables                = spec$data$num_variables,
-    num_cases                    = spec$data$num_cases,
-    na_impute                    = spec$missing$na_impute,
-    variable_type                = spec$variables$variable_type,
-    iter                         = spec$sampler$iter,
-    warmup                       = spec$sampler$warmup,
-    edge_selection               = spec$prior$edge_selection,
-    edge_prior                   = spec$prior$edge_prior,
-    graph_prior_spec             = spec$prior$graph_prior_spec,
-    calibration_window           = spec$prior$calibration_window,
-    inclusion_probability        = spec$prior$inclusion_probability,
-    beta_bernoulli_alpha         = spec$prior$beta_bernoulli_alpha,
-    beta_bernoulli_beta          = spec$prior$beta_bernoulli_beta,
+    num_variables = spec$data$num_variables,
+    num_cases = spec$data$num_cases,
+    na_impute = spec$missing$na_impute,
+    variable_type = spec$variables$variable_type,
+    iter = spec$sampler$iter,
+    warmup = spec$sampler$warmup,
+    edge_selection = spec$prior$edge_selection,
+    edge_prior = spec$prior$edge_prior,
+    precision_graph_prior = spec$prior$precision_graph_prior,
+    calibration_window = spec$prior$calibration_window,
+    inclusion_probability = spec$prior$inclusion_probability,
+    beta_bernoulli_alpha = spec$prior$beta_bernoulli_alpha,
+    beta_bernoulli_beta = spec$prior$beta_bernoulli_beta,
     beta_bernoulli_alpha_between = spec$prior$beta_bernoulli_alpha_between,
-    beta_bernoulli_beta_between  = spec$prior$beta_bernoulli_beta_between,
-    dirichlet_alpha              = spec$prior$dirichlet_alpha,
-    lambda                       = spec$prior$lambda,
-    na_action                    = spec$missing$na_action,
-    version                      = packageVersion("bgms"),
-    update_method                = spec$sampler$update_method,
-    target_accept                = spec$sampler$target_accept,
-    num_chains                   = spec$sampler$chains,
-    data_columnnames             = spec$data$data_columnnames,
-    no_variables                 = spec$data$num_variables,
-    column_means                 = spec$data$column_means,
-    is_continuous                = TRUE,
-    model_type                   = "ggm"
+    beta_bernoulli_beta_between = spec$prior$beta_bernoulli_beta_between,
+    dirichlet_alpha = spec$prior$dirichlet_alpha,
+    lambda = spec$prior$lambda,
+    na_action = spec$missing$na_action,
+    version = packageVersion("bgms"),
+    update_method = spec$sampler$update_method,
+    target_accept = spec$sampler$target_accept,
+    num_chains = spec$sampler$chains,
+    data_columnnames = spec$data$data_columnnames,
+    no_variables = spec$data$num_variables,
+    column_means = spec$data$column_means,
+    is_continuous = TRUE,
+    model_type = "ggm"
   )
 }
 
@@ -105,47 +104,47 @@ build_arguments_omrf = function(spec) {
 
 build_arguments_mixed_mrf = function(spec) {
   list(
-    num_variables                = spec$data$num_variables,
-    num_discrete                 = spec$data$num_discrete,
-    num_continuous               = spec$data$num_continuous,
-    num_cases                    = spec$data$num_cases,
-    variable_type                = spec$variables$variable_type,
-    iter                         = spec$sampler$iter,
-    warmup                       = spec$sampler$warmup,
-    pairwise_scale               = spec$prior$pairwise_scale,
-    main_alpha                   = spec$prior$main_alpha,
-    main_beta                    = spec$prior$main_beta,
-    edge_selection               = spec$prior$edge_selection,
-    edge_prior                   = spec$prior$edge_prior,
-    graph_prior_spec             = spec$prior$graph_prior_spec,
-    calibration_window           = spec$prior$calibration_window,
-    inclusion_probability        = spec$prior$inclusion_probability,
-    beta_bernoulli_alpha         = spec$prior$beta_bernoulli_alpha,
-    beta_bernoulli_beta          = spec$prior$beta_bernoulli_beta,
+    num_variables = spec$data$num_variables,
+    num_discrete = spec$data$num_discrete,
+    num_continuous = spec$data$num_continuous,
+    num_cases = spec$data$num_cases,
+    variable_type = spec$variables$variable_type,
+    iter = spec$sampler$iter,
+    warmup = spec$sampler$warmup,
+    pairwise_scale = spec$prior$pairwise_scale,
+    main_alpha = spec$prior$main_alpha,
+    main_beta = spec$prior$main_beta,
+    edge_selection = spec$prior$edge_selection,
+    edge_prior = spec$prior$edge_prior,
+    precision_graph_prior = spec$prior$precision_graph_prior,
+    calibration_window = spec$prior$calibration_window,
+    inclusion_probability = spec$prior$inclusion_probability,
+    beta_bernoulli_alpha = spec$prior$beta_bernoulli_alpha,
+    beta_bernoulli_beta = spec$prior$beta_bernoulli_beta,
     beta_bernoulli_alpha_between = spec$prior$beta_bernoulli_alpha_between,
-    beta_bernoulli_beta_between  = spec$prior$beta_bernoulli_beta_between,
-    dirichlet_alpha              = spec$prior$dirichlet_alpha,
-    lambda                       = spec$prior$lambda,
-    na_action                    = spec$missing$na_action,
-    na_impute                    = spec$missing$na_impute,
-    version                      = packageVersion("bgms"),
-    update_method                = spec$sampler$update_method,
-    target_accept                = spec$sampler$target_accept,
-    nuts_max_depth               = spec$sampler$nuts_max_depth,
-    num_chains                   = spec$sampler$chains,
-    num_categories               = spec$data$num_categories,
-    category_levels              = spec$data$category_levels,
-    blume_capel_shift            = spec$data$blume_capel_shift,
-    data_columnnames             = spec$data$data_columnnames,
-    data_columnnames_discrete    = spec$data$data_columnnames_discrete,
-    data_columnnames_continuous  = spec$data$data_columnnames_continuous,
-    discrete_indices             = spec$data$discrete_indices,
-    continuous_indices           = spec$data$continuous_indices,
-    baseline_category            = spec$variables$baseline_category,
-    is_ordinal                   = spec$variables$is_ordinal,
-    no_variables                 = spec$data$num_variables,
-    is_mixed                     = TRUE,
-    model_type                   = "mixed_mrf"
+    beta_bernoulli_beta_between = spec$prior$beta_bernoulli_beta_between,
+    dirichlet_alpha = spec$prior$dirichlet_alpha,
+    lambda = spec$prior$lambda,
+    na_action = spec$missing$na_action,
+    na_impute = spec$missing$na_impute,
+    version = packageVersion("bgms"),
+    update_method = spec$sampler$update_method,
+    target_accept = spec$sampler$target_accept,
+    nuts_max_depth = spec$sampler$nuts_max_depth,
+    num_chains = spec$sampler$chains,
+    num_categories = spec$data$num_categories,
+    category_levels = spec$data$category_levels,
+    blume_capel_shift = spec$data$blume_capel_shift,
+    data_columnnames = spec$data$data_columnnames,
+    data_columnnames_discrete = spec$data$data_columnnames_discrete,
+    data_columnnames_continuous = spec$data$data_columnnames_continuous,
+    discrete_indices = spec$data$discrete_indices,
+    continuous_indices = spec$data$continuous_indices,
+    baseline_category = spec$variables$baseline_category,
+    is_ordinal = spec$variables$is_ordinal,
+    no_variables = spec$data$num_variables,
+    is_mixed = TRUE,
+    model_type = "mixed_mrf"
   )
 }
 
diff --git a/R/build_output_bgm.R b/R/build_output_bgm.R
index 44d1a856..b632f4d7 100644
--- a/R/build_output_bgm.R
+++ b/R/build_output_bgm.R
@@ -26,7 +26,7 @@ build_output_bgm = function(spec, raw) {
   # Keep the raw chains for the Z-ratio alarm suite: it needs the untouched
   # indicator layout and the per-chain zratio block, both dropped by the
   # normalization below.
-  zratio_chains = if(identical(p$graph_prior_spec, "hierarchical")) {
+  zratio_chains = if(identical(p$precision_graph_prior, "hierarchical")) {
     raw
   } else {
     NULL
diff --git a/R/build_output_mixed_mrf.R b/R/build_output_mixed_mrf.R
index 35439ad9..0b42c0b3 100644
--- a/R/build_output_mixed_mrf.R
+++ b/R/build_output_mixed_mrf.R
@@ -38,7 +38,7 @@ build_output_mixed_mrf = function(spec, raw) {
   # Keep the raw chains for the Z-ratio alarm suite: it needs the untouched
   # indicator layout and the per-chain zratio block, both dropped by the
   # normalization below.
-  zratio_chains = if(identical(pr$graph_prior_spec, "hierarchical")) {
+  zratio_chains = if(identical(pr$precision_graph_prior, "hierarchical")) {
     raw
   } else {
     NULL
diff --git a/R/build_spec.R b/R/build_spec.R
index b89e4a19..b9a99d0c 100644
--- a/R/build_spec.R
+++ b/R/build_spec.R
@@ -29,7 +29,6 @@ edge_prior_spec_fields = function(ep) {
 }
 
 
-
 # ==============================================================================
 # sampler_sublist()  --- extract validated sampler list for new_bgm_spec()
 # ==============================================================================
@@ -51,7 +50,6 @@ sampler_sublist = function(s) {
 }
 
 
-
 build_spec_ggm = function(x, data_columnnames, num_variables,
                           variable_type, is_ordinal, is_continuous,
                           baseline_category,
@@ -61,7 +59,7 @@ build_spec_ggm = function(x, data_columnnames, num_variables,
                           scale_prior_type, scale_shape, scale_rate,
                           scale_eta = NA_real_,
                           delta = 0,
-                          graph_prior_spec = "joint",
+                          precision_graph_prior = "joint",
                           calibration_window = NULL,
                           edge_prior_flat) {
   # Missing data
@@ -112,7 +110,7 @@ build_spec_ggm = function(x, data_columnnames, num_variables,
         scale_rate = scale_rate,
         scale_eta = scale_eta,
         delta = delta,
-        graph_prior_spec = graph_prior_spec,
+        precision_graph_prior = precision_graph_prior,
         calibration_window = calibration_window
       ),
       edge_prior_spec_fields(ep)
@@ -228,7 +226,7 @@ build_spec_mixed_mrf = function(x, data_columnnames, num_variables,
                                 scale_prior_type, scale_shape, scale_rate,
                                 scale_eta = NA_real_,
                                 delta = 0,
-                                graph_prior_spec = "joint",
+                                precision_graph_prior = "joint",
                                 calibration_window = NULL,
                                 edge_prior_flat) {
   # Standardized-frame scale prior: derive the raw diagonal rate eta / s
@@ -375,7 +373,7 @@ build_spec_mixed_mrf = function(x, data_columnnames, num_variables,
         scale_rate = scale_rate,
         scale_eta = scale_eta,
         delta = delta,
-        graph_prior_spec = graph_prior_spec,
+        precision_graph_prior = precision_graph_prior,
         calibration_window = calibration_window
       ),
       edge_prior_spec_fields(ep)
diff --git a/R/diagnostics_nuts.R b/R/diagnostics_nuts.R
index a44bb9e3..c78b28b4 100644
--- a/R/diagnostics_nuts.R
+++ b/R/diagnostics_nuts.R
@@ -206,7 +206,7 @@ summarize_nuts_diagnostics = function(out, nuts_max_depth = 10, verbose = TRUE)
     incomplete_chains = which(warmup_check$warmup_incomplete)
     if(length(incomplete_chains) > 0) {
       issues = c(issues, sprintf(
-        "Warmup incomplete: energy not stationary in chain%s %s - increase warmup",
+        "Warmup may be incomplete: energy not stationary in chain%s %s - check R-hat and ESS",
         if(length(incomplete_chains) > 1) "s" else "",
         paste(incomplete_chains, collapse = ", ")
       ))
diff --git a/R/run_sampler.R b/R/run_sampler.R
index bba332ef..f48cd716 100644
--- a/R/run_sampler.R
+++ b/R/run_sampler.R
@@ -77,7 +77,7 @@ run_sampler_ggm = function(spec) {
   # are mutually exclusive.
   correction = NULL
   zratio = NULL
-  if(identical(p$graph_prior_spec, "hierarchical")) {
+  if(identical(p$precision_graph_prior, "hierarchical")) {
     zc = zratio_constants(
       delta = p$delta,
       sigma = 2 * p$pairwise_scale,
@@ -221,7 +221,7 @@ run_sampler_mixed_mrf = function(spec) {
   # window are sized on the continuous subgraph.
   correction = NULL
   zratio = NULL
-  if(identical(p$graph_prior_spec, "hierarchical")) {
+  if(identical(p$precision_graph_prior, "hierarchical")) {
     zc = zratio_constants(
       delta = p$delta,
       sigma = 2 * p$pairwise_scale,
diff --git a/R/sample_ggm_prior.R b/R/sample_ggm_prior.R
index 57e0198e..18cffb4e 100644
--- a/R/sample_ggm_prior.R
+++ b/R/sample_ggm_prior.R
@@ -315,18 +315,26 @@ sample_ggm_prior = function(
     # Z-ratio engine carries the normalizer into the between-edge moves,
     # and the hyperparameter updates are the clean conjugate draws (no
     # C-correction on this path).
-    if(!identical(ip$interaction_prior_type, "normal")) {
-      stop(
-        "spec = \"hierarchical\" requires a normal interaction (slab) ",
-        "prior; the Z-ratio constants are derived for the Normal slab."
-      )
+    if(!ip$interaction_prior_type %in% c("normal", "cauchy")) {
+      stop(sprintf(
+        paste0(
+          "spec = \"hierarchical\" supports a normal or Cauchy interaction ",
+          "(slab) prior. Got %s_prior(). Use interaction_prior = ",
+          "normal_prior() or cauchy_prior()."
+        ),
+        ip$interaction_prior_type
+      ))
     }
     if(abs(sp$scale_shape - 1) > 1e-12) {
-      stop(
-        "spec = \"hierarchical\" requires shape = 1 on the diagonal scale ",
-        "prior (gamma_prior(shape = 1) or exponential_prior); the Z-ratio ",
-        "constants are derived for the exponential diagonal."
-      )
+      stop(sprintf(
+        paste0(
+          "spec = \"hierarchical\" requires shape = 1 on the diagonal ",
+          "scale prior; the Z-ratio normalizer is derived for the ",
+          "exponential diagonal. Got shape = %s. Use ",
+          "gamma_prior(shape = 1) or exponential_prior()."
+        ),
+        format(sp$scale_shape)
+      ))
     }
     zc = zratio_constants(
       delta = delta,
diff --git a/man/bgm.Rd b/man/bgm.Rd
index 4a042cbe..9b9ea71f 100644
--- a/man/bgm.Rd
+++ b/man/bgm.Rd
@@ -17,7 +17,7 @@ bgm(
   delta = NULL,
   edge_selection = TRUE,
   edge_prior = bernoulli_prior(0.5),
-  graph_prior_spec = c("joint", "hierarchical"),
+  precision_graph_prior = c("joint", "hierarchical"),
   calibration_window = NULL,
   na_action = c("listwise", "impute"),
   update_method = c("nuts", "adaptive-metropolis", "gibbs"),
@@ -181,7 +181,7 @@ configuration skip the build. With \code{beta_bernoulli_prior()} the
 sampled inclusion probability is returned per chain in
 \code{fit$inclusion_parameter_samples}.}
 
-\item{graph_prior_spec}{Character. How the precision prior composes with
+\item{precision_graph_prior}{Character. How the precision prior composes with
 the edge prior under edge selection for continuous (GGM) data:
 \describe{
 \item{"joint"}{(default) The un-normalised joint specification
@@ -202,10 +202,10 @@ warm-up window (see \code{calibration_window}) and audited after
 sampling (\code{\link{summarize_zratio_diagnostics}}; the summary
 is returned as \code{fit$zratio_diag} and issues print like other
 sampler warnings). Requires \code{edge_selection = TRUE}, a
-\code{normal_prior()} interaction prior, a shape-1
-\code{gamma_prior()} (or \code{exponential_prior()}) precision
-scale prior, and continuous data — either all-continuous (GGM) or
-mixed with at least two continuous variables. On mixed data the
+\code{normal_prior()} or \code{cauchy_prior()} interaction prior, a
+shape-1 \code{gamma_prior()} (or \code{exponential_prior()})
+precision scale prior, and continuous data — either all-continuous
+(GGM) or mixed with at least two continuous variables. On mixed data the
 normalizer lives on the continuous block \eqn{K_{yy}}, so the
 Z-ratio enters the continuous-continuous edge moves only, with the
 mediating-block counts read off the continuous subgraph; discrete
@@ -214,7 +214,7 @@ and cross edges are unchanged.}
 Default: \code{"joint"}.}
 
 \item{calibration_window}{Non-negative integer or \code{NULL} (default).
-Only for \code{graph_prior_spec = "hierarchical"}: length of the
+Only for \code{precision_graph_prior = "hierarchical"}: length of the
 appended warm-up window in which the Z-ratio correction is calibrated
 online against a block-Gibbs oracle and then frozen (with its hull
 clamp) before sampling. \code{NULL} resolves to no window for
diff --git a/tests/testthat/test-bgm-hier-spec.R b/tests/testthat/test-bgm-hier-spec.R
index 2ae962b5..4c325fca 100644
--- a/tests/testthat/test-bgm-hier-spec.R
+++ b/tests/testthat/test-bgm-hier-spec.R
@@ -1,4 +1,4 @@
-# Tests for bgm(graph_prior_spec = "hierarchical"): eligibility validation
+# Tests for bgm(precision_graph_prior = "hierarchical"): eligibility validation
 # and the end-to-end fit with the Z-ratio engine and alarm suite attached.
 
 hier_test_data = function(q = 10, n = 40, seed = 4) {
@@ -12,18 +12,18 @@ test_that("hierarchical spec eligibility is validated", {
   expect_error(
     bgm(
       x = Y, variable_type = "continuous",
-      interaction_prior = cauchy_prior(scale = 0.5),
-      graph_prior_spec = "hierarchical",
+      interaction_prior = beta_prime_prior(),
+      precision_graph_prior = "hierarchical",
       update_method = "gibbs", display_progress = "none", verbose = FALSE
     ),
-    "normal"
+    "normal or Cauchy"
   )
   expect_error(
     bgm(
       x = Y, variable_type = "continuous",
       interaction_prior = normal_prior(scale = 0.5),
       precision_scale_prior = gamma_prior(shape = 2, rate = 2),
-      graph_prior_spec = "hierarchical",
+      precision_graph_prior = "hierarchical",
       update_method = "gibbs", display_progress = "none", verbose = FALSE
     ),
     "shape = 1"
@@ -34,7 +34,7 @@ test_that("hierarchical spec eligibility is validated", {
       interaction_prior = normal_prior(scale = 0.5),
       precision_scale_prior = gamma_prior(shape = 1, rate = 2),
       edge_selection = FALSE,
-      graph_prior_spec = "hierarchical",
+      precision_graph_prior = "hierarchical",
       update_method = "gibbs", display_progress = "none", verbose = FALSE
     ),
     "edge_selection"
@@ -43,13 +43,36 @@ test_that("hierarchical spec eligibility is validated", {
     bgm(
       x = matrix(sample(0:3, 200, replace = TRUE), 50, 4),
       interaction_prior = normal_prior(scale = 0.5),
-      graph_prior_spec = "hierarchical",
+      precision_graph_prior = "hierarchical",
       display_progress = "none", verbose = FALSE
     ),
     "continuous"
   )
 })
 
+test_that("the hierarchical spec accepts a Cauchy slab on every update method", {
+  skip_on_cran()
+  Y = hier_test_data(q = 8)
+  for(method in c("nuts", "adaptive-metropolis", "gibbs")) {
+    fit = bgm(
+      x = Y, variable_type = "continuous",
+      iter = 100, warmup = 150,
+      interaction_prior = cauchy_prior(scale = 0.5),
+      precision_scale_prior = gamma_prior(shape = 1, rate = 2),
+      precision_graph_prior = "hierarchical", calibration_window = 50,
+      update_method = method, chains = 1, cores = 1, seed = 7,
+      display_progress = "none", verbose = FALSE
+    )
+    s = summary(fit)
+    expect_true(all(is.finite(s$pairwise$mean)), info = method)
+    expect_true(all(s$indicator$mean >= 0 & s$indicator$mean <= 1),
+      info = method
+    )
+    expect_equal(fit@arguments$precision_graph_prior, "hierarchical", info = method)
+    expect_false(is.null(fit@zratio_diag), info = method)
+  }
+})
+
 test_that("bgm fits the hierarchical spec and attaches the alarm suite", {
   skip_on_cran()
   Y = hier_test_data(q = 12)
@@ -59,7 +82,7 @@ test_that("bgm fits the hierarchical spec and attaches the alarm suite", {
     interaction_prior = normal_prior(scale = 0.5),
     precision_scale_prior = gamma_prior(shape = 1, rate = 2),
     edge_prior = beta_bernoulli_prior(2, 4),
-    graph_prior_spec = "hierarchical", calibration_window = 100,
+    precision_graph_prior = "hierarchical", calibration_window = 100,
     update_method = "gibbs", chains = 2, cores = 2, seed = 11,
     display_progress = "none", verbose = FALSE
   )
@@ -70,7 +93,7 @@ test_that("bgm fits the hierarchical spec and attaches the alarm suite", {
   expect_true(all(zd$per_chain$n_oracle > 0))
   expect_false(zd$verdict_flagged)
   expect_equal(length(zd$audits), 2L)
-  expect_equal(fit@arguments$graph_prior_spec, "hierarchical")
+  expect_equal(fit@arguments$precision_graph_prior, "hierarchical")
   # The joint-path hyperparameter correction must not run on this path;
   # inclusion-parameter samples come from the clean conjugate draw.
   expect_equal(length(fit@inclusion_parameter_samples), 2L)
@@ -85,7 +108,7 @@ test_that("the joint default is unchanged", {
     update_method = "gibbs", chains = 1, cores = 1, seed = 3,
     display_progress = "none", verbose = FALSE
   )
-  expect_equal(fit@arguments$graph_prior_spec, "joint")
+  expect_equal(fit@arguments$precision_graph_prior, "joint")
   expect_null(fit@zratio_diag)
 })
 
@@ -125,7 +148,7 @@ test_that("mixed data supports the hierarchical spec on the continuous block", {
       x = X[, 1:4], variable_type = vt[1:4],
       interaction_prior = normal_prior(scale = 0.5),
       precision_scale_prior = gamma_prior(shape = 1, rate = 2),
-      graph_prior_spec = "hierarchical",
+      precision_graph_prior = "hierarchical",
       display_progress = "none", verbose = FALSE
     ),
     "two continuous"
@@ -137,7 +160,7 @@ test_that("mixed data supports the hierarchical spec on the continuous block", {
     interaction_prior = normal_prior(scale = 0.5),
     precision_scale_prior = gamma_prior(shape = 1, rate = 2),
     edge_prior = beta_bernoulli_prior(2, 4),
-    graph_prior_spec = "hierarchical", calibration_window = 80,
+    precision_graph_prior = "hierarchical", calibration_window = 80,
     update_method = "adaptive-metropolis", chains = 1, cores = 1, seed = 5,
     display_progress = "none", verbose = FALSE
   )
@@ -145,5 +168,5 @@ test_that("mixed data supports the hierarchical spec on the continuous block", {
   expect_false(is.null(zd))
   expect_true(zd$per_chain$frozen)
   expect_false(zd$verdict_flagged)
-  expect_equal(fit@arguments$graph_prior_spec, "hierarchical")
+  expect_equal(fit@arguments$precision_graph_prior, "hierarchical")
 })
diff --git a/tests/testthat/test-build-arguments.R b/tests/testthat/test-build-arguments.R
index 4f17c738..56a6a0b0 100644
--- a/tests/testthat/test-build-arguments.R
+++ b/tests/testthat/test-build-arguments.R
@@ -94,7 +94,7 @@ test_that("GGM build_arguments: all expected field names present", {
   expected = c(
     "num_variables", "num_cases", "na_impute", "variable_type",
     "iter", "warmup", "edge_selection", "edge_prior",
-    "graph_prior_spec", "calibration_window",
+    "precision_graph_prior", "calibration_window",
     "inclusion_probability", "beta_bernoulli_alpha", "beta_bernoulli_beta",
     "beta_bernoulli_alpha_between", "beta_bernoulli_beta_between",
     "dirichlet_alpha", "lambda", "na_action", "version",
diff --git a/tests/testthat/test-hier-zratio-identity.R b/tests/testthat/test-hier-zratio-identity.R
index 01f595c2..abc9c877 100644
--- a/tests/testthat/test-hier-zratio-identity.R
+++ b/tests/testthat/test-hier-zratio-identity.R
@@ -20,14 +20,7 @@ hier_prior_run = function(q, delta, sigma, p_inc, um, edge_prior = NULL,
   )
 }
 
-test_that("hierarchical spec requires the Normal slab and shape-1 diagonal", {
-  expect_error(
-    sample_ggm_prior(
-      p = 4, n_samples = 5, n_warmup = 5, spec = "hierarchical",
-      interaction_prior = cauchy_prior(scale = 1), verbose = FALSE
-    ),
-    "normal interaction"
-  )
+test_that("hierarchical spec requires the shape-1 diagonal", {
   expect_error(
     sample_ggm_prior(
       p = 4, n_samples = 5, n_warmup = 5, spec = "hierarchical",