From 23695fec2397b7f1d38004fb7fa7cd170acb381a Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Sun, 5 Jul 2026 20:22:27 +0200 Subject: [PATCH 1/5] Accept a Cauchy slab under the hierarchical graph-prior specification The eligibility check rejected every slab except the normal, but the Cauchy slab is a scale mixture of normals: conditional on the per-edge mixture weight it is a normal slab, so the Z-ratio normalizer applies unchanged. The Gibbs edge move already conditions on the weights, and the Metropolis and NUTS between-edge moves evaluate the marginal Cauchy density, which integrates the weights out; all three samplers target the same posterior. On a synthetic check their inclusion probabilities agree to within Monte Carlo noise (max difference 0.009 across methods and edges). The eligibility messages now name the received prior and the supported alternatives. --- R/bgm.R | 9 ++++---- R/bgm_spec.R | 36 ++++++++++++++++++----------- R/sample_ggm_prior.R | 30 ++++++++++++++++-------- man/bgm.Rd | 9 ++++---- tests/testthat/test-bgm-hier-spec.R | 27 ++++++++++++++++++++-- 5 files changed, 77 insertions(+), 34 deletions(-) diff --git a/R/bgm.R b/R/bgm.R index 862c12a5..7095bf43 100644 --- a/R/bgm.R +++ b/R/bgm.R @@ -189,10 +189,11 @@ #' 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 +#' (the Cauchy slab enters through its scale mixture of normals), 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 diff --git a/R/bgm_spec.R b/R/bgm_spec.R index 1b878eb4..92d23e40 100644 --- a/R/bgm_spec.R +++ b/R/bgm_spec.R @@ -403,22 +403,30 @@ bgm_spec = function(x, "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( + "graph_prior_spec = \"hierarchical\" supports a normal or Cauchy ", + "interaction (slab) prior: the Z-ratio normalizer is derived for ", + "the normal slab, and the Cauchy slab enters through its scale ", + "mixture of normals. Got %s_prior(). Use interaction_prior = ", + "normal_prior() or cauchy_prior(), or keep graph_prior_spec = ", + "\"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( + "graph_prior_spec = \"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 ", + "graph_prior_spec = \"joint\"." + ), + format(scale_shape) + )) } } diff --git a/R/sample_ggm_prior.R b/R/sample_ggm_prior.R index 57e0198e..e550ea80 100644 --- a/R/sample_ggm_prior.R +++ b/R/sample_ggm_prior.R @@ -315,18 +315,28 @@ 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: the Z-ratio normalizer is derived for the normal ", + "slab, and the Cauchy slab enters through its scale mixture of ", + "normals. 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..557d28ba 100644 --- a/man/bgm.Rd +++ b/man/bgm.Rd @@ -202,10 +202,11 @@ 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 +(the Cauchy slab enters through its scale mixture of normals), 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 diff --git a/tests/testthat/test-bgm-hier-spec.R b/tests/testthat/test-bgm-hier-spec.R index 2ae962b5..40bbe00d 100644 --- a/tests/testthat/test-bgm-hier-spec.R +++ b/tests/testthat/test-bgm-hier-spec.R @@ -12,11 +12,11 @@ test_that("hierarchical spec eligibility is validated", { expect_error( bgm( x = Y, variable_type = "continuous", - interaction_prior = cauchy_prior(scale = 0.5), + interaction_prior = beta_prime_prior(), graph_prior_spec = "hierarchical", update_method = "gibbs", display_progress = "none", verbose = FALSE ), - "normal" + "normal or Cauchy" ) expect_error( bgm( @@ -50,6 +50,29 @@ test_that("hierarchical spec eligibility is validated", { ) }) +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), + graph_prior_spec = "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$graph_prior_spec, "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) From 6cae5c5e4f9923418efc49e01e91e16e26883e2d Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Sun, 5 Jul 2026 20:34:50 +0200 Subject: [PATCH 2/5] Rename graph_prior_spec to precision_graph_prior The argument selects how the prior on the continuous precision block composes with the graph: a joint spike-and-slab prior on (K, Gamma) or the hierarchical factorization p(Gamma) p(K | Gamma). The old name read as configuring the prior on the graph itself, which is edge_prior's role. The argument has not appeared in a CRAN release, so the rename carries no deprecation shim. Also updates the prior-chain eligibility test for the Cauchy slab now being accepted. --- NEWS.md | 2 +- R/bgm.R | 92 +++---- R/bgm_spec.R | 264 ++++++++++---------- R/build_arguments.R | 136 +++++------ R/build_output_bgm.R | 272 ++++++++++----------- R/build_output_mixed_mrf.R | 256 +++++++++---------- R/build_spec.R | 152 ++++++------ R/run_sampler.R | 60 ++--- man/bgm.Rd | 6 +- tests/testthat/test-bgm-hier-spec.R | 58 ++--- tests/testthat/test-build-arguments.R | 8 +- tests/testthat/test-hier-zratio-identity.R | 9 +- 12 files changed, 654 insertions(+), 661 deletions(-) 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 7095bf43..3baf6ec6 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 @@ -202,7 +202,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 @@ -395,7 +395,7 @@ #' } #' #' @export -bgm = function( +bgm <- function( x, variable_type = "ordinal", baseline_category, @@ -408,7 +408,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"), @@ -441,87 +441,87 @@ bgm = function( ) { # Set verbose option for internal functions, restore on exit - old_verbose = getOption("bgms.verbose") + old_verbose <- getOption("bgms.verbose") options(bgms.verbose = verbose) on.exit(options(bgms.verbose = old_verbose), add = TRUE) # --- Legacy deprecation: v0.1.6.0 renames ----------------------------------- - if(hasArg(interaction_scale)) { + if (hasArg(interaction_scale)) { lifecycle::deprecate_warn( "0.1.6.0", "bgm(interaction_scale =)", "bgm(interaction_prior =)" ) - if(!hasArg(pairwise_scale) && + if (!hasArg(pairwise_scale) && identical(interaction_prior, cauchy_prior(scale = 1))) { - interaction_prior = cauchy_prior(scale = interaction_scale) + interaction_prior <- cauchy_prior(scale = interaction_scale) } } - if(hasArg(burnin)) { + if (hasArg(burnin)) { lifecycle::deprecate_warn("0.1.6.0", "bgm(burnin =)", "bgm(warmup =)") - if(!hasArg(warmup)) warmup = burnin + if (!hasArg(warmup)) warmup <- burnin } - if(hasArg(save)) { + if (hasArg(save)) { lifecycle::deprecate_warn("0.1.6.0", "bgm(save =)") } - if(hasArg(threshold_alpha) || hasArg(threshold_beta)) { + if (hasArg(threshold_alpha) || hasArg(threshold_beta)) { lifecycle::deprecate_warn( "0.1.6.0", "bgm(threshold_alpha =, threshold_beta =)", "bgm(threshold_prior =)" ) - if(identical(threshold_prior, beta_prime_prior(0.5, 0.5))) { - ta = if(hasArg(threshold_alpha)) threshold_alpha else 0.5 - tb = if(hasArg(threshold_beta)) threshold_beta else 0.5 - threshold_prior = beta_prime_prior(alpha = ta, beta = tb) + if (identical(threshold_prior, beta_prime_prior(0.5, 0.5))) { + ta <- if (hasArg(threshold_alpha)) threshold_alpha else 0.5 + tb <- if (hasArg(threshold_beta)) threshold_beta else 0.5 + threshold_prior <- beta_prime_prior(alpha = ta, beta = tb) } } # --- Legacy deprecation: scalar prior parameters ---------------------------- - if(hasArg(pairwise_scale)) { + if (hasArg(pairwise_scale)) { lifecycle::deprecate_warn( "0.2.0", "bgm(pairwise_scale =)", "bgm(interaction_prior =)" ) - if(identical(interaction_prior, cauchy_prior(scale = 1))) { - interaction_prior = cauchy_prior(scale = pairwise_scale) + if (identical(interaction_prior, cauchy_prior(scale = 1))) { + interaction_prior <- cauchy_prior(scale = pairwise_scale) } } - if(hasArg(main_alpha) || hasArg(main_beta)) { + if (hasArg(main_alpha) || hasArg(main_beta)) { lifecycle::deprecate_warn( "0.2.0", "bgm(main_alpha =)", "bgm(threshold_prior =)" ) - if(identical(threshold_prior, beta_prime_prior(0.5, 0.5))) { - ma = if(hasArg(main_alpha)) main_alpha else 0.5 - mb = if(hasArg(main_beta)) main_beta else 0.5 - threshold_prior = beta_prime_prior(alpha = ma, beta = mb) + if (identical(threshold_prior, beta_prime_prior(0.5, 0.5))) { + ma <- if (hasArg(main_alpha)) main_alpha else 0.5 + mb <- if (hasArg(main_beta)) main_beta else 0.5 + threshold_prior <- beta_prime_prior(alpha = ma, beta = mb) } } # Handle edge_prior: accept both string (deprecated) and object (new) - if(is.character(edge_prior)) { + if (is.character(edge_prior)) { lifecycle::deprecate_warn( "0.2.0", "bgm(edge_prior = 'must be a prior object')", "bgm(edge_prior = 'bernoulli_prior()')" ) - edge_prior_str = match.arg(edge_prior, + edge_prior_str <- match.arg(edge_prior, choices = c("Bernoulli", "Beta-Bernoulli", "Stochastic-Block") ) - ip = if(hasArg(inclusion_probability)) inclusion_probability else 0.5 - bba = if(hasArg(beta_bernoulli_alpha)) beta_bernoulli_alpha else 1 - bbb = if(hasArg(beta_bernoulli_beta)) beta_bernoulli_beta else 1 - bbab = if(hasArg(beta_bernoulli_alpha_between)) beta_bernoulli_alpha_between else 1 - bbbb = if(hasArg(beta_bernoulli_beta_between)) beta_bernoulli_beta_between else 1 - da = if(hasArg(dirichlet_alpha)) dirichlet_alpha else 1 - lam = if(hasArg(lambda)) lambda else 1 + ip <- if (hasArg(inclusion_probability)) inclusion_probability else 0.5 + bba <- if (hasArg(beta_bernoulli_alpha)) beta_bernoulli_alpha else 1 + bbb <- if (hasArg(beta_bernoulli_beta)) beta_bernoulli_beta else 1 + bbab <- if (hasArg(beta_bernoulli_alpha_between)) beta_bernoulli_alpha_between else 1 + bbbb <- if (hasArg(beta_bernoulli_beta_between)) beta_bernoulli_beta_between else 1 + da <- if (hasArg(dirichlet_alpha)) dirichlet_alpha else 1 + lam <- if (hasArg(lambda)) lambda else 1 - edge_prior = switch(edge_prior_str, + edge_prior <- switch(edge_prior_str, "Bernoulli" = bernoulli_prior(inclusion_probability = ip), "Beta-Bernoulli" = beta_bernoulli_prior(alpha = bba, beta = bbb), "Stochastic-Block" = sbm_prior( @@ -532,13 +532,13 @@ bgm = function( ) } else { # Warn if loose edge params are also supplied alongside an object - if(hasArg(inclusion_probability)) { + if (hasArg(inclusion_probability)) { lifecycle::deprecate_warn( "0.2.0", "bgm(inclusion_probability =)", "bgm(edge_prior = 'bernoulli_prior()')" ) } - if(hasArg(beta_bernoulli_alpha)) { + if (hasArg(beta_bernoulli_alpha)) { lifecycle::deprecate_warn( "0.2.0", "bgm(beta_bernoulli_alpha =)", "bgm(edge_prior = 'beta_bernoulli_prior()')" @@ -547,17 +547,17 @@ bgm = function( } # --- Unpack prior objects to flat parameters for bgm_spec ------------------- - ip = unpack_interaction_prior(interaction_prior) - tp = unpack_threshold_prior(threshold_prior) - mp = unpack_parameter_prior(means_prior) - sp = unpack_scale_prior(precision_scale_prior) + ip <- unpack_interaction_prior(interaction_prior) + tp <- unpack_threshold_prior(threshold_prior) + mp <- unpack_parameter_prior(means_prior) + sp <- unpack_scale_prior(precision_scale_prior) # --- Build spec, sample, build output ---------------------------------------- - spec = bgm_spec( + spec <- bgm_spec( x = x, model_type = "omrf", variable_type = variable_type, - baseline_category = if(hasArg(baseline_category)) baseline_category else NULL, + baseline_category = if (hasArg(baseline_category)) baseline_category else NULL, na_action = na_action, interaction_prior_type = ip$interaction_prior_type, pairwise_scale = ip$pairwise_scale, @@ -578,10 +578,10 @@ 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, + target_accept = if (hasArg(target_accept)) target_accept else NULL, iter = iter, warmup = warmup, nuts_max_depth = nuts_max_depth, @@ -594,7 +594,7 @@ bgm = function( progress_callback = progress_callback ) - raw = run_sampler(spec) - output = build_output(spec, raw) + raw <- run_sampler(spec) + output <- build_output(spec, raw) return(output) } diff --git a/R/bgm_spec.R b/R/bgm_spec.R index 92d23e40..f4b4f6f1 100644 --- a/R/bgm_spec.R +++ b/R/bgm_spec.R @@ -238,79 +238,79 @@ validate_bgm_spec = function(spec) { # # Parameters mirror the union of bgm() and bgmCompare() arguments. # ============================================================================== -bgm_spec = function(x, - model_type = c("omrf", "ggm", "compare", "mixed_mrf"), - # Variable specification - variable_type = "ordinal", - baseline_category = NULL, - # Data (compare-specific) - y = NULL, - group_indicator = NULL, - # Missing data - na_action = c("listwise", "impute"), - # Priors (new: prior objects unpacked by bgm()) - interaction_prior_type = "cauchy", - pairwise_scale = 1, - interaction_alpha = NA_real_, - interaction_beta = NA_real_, - threshold_prior_type = "beta-prime", - main_alpha = 0.5, - main_beta = 0.5, - threshold_scale = NA_real_, - means_prior_type = "normal", - means_scale = 1, - means_alpha = NA_real_, - means_beta = NA_real_, - scale_prior_type = "gamma", - scale_shape = 1, - scale_rate = 1, - scale_eta = NA_real_, - delta = NULL, - edge_selection = TRUE, - edge_prior = bernoulli_prior(0.5), - graph_prior_spec = c("joint", "hierarchical"), - calibration_window = NULL, - # Legacy edge prior params (accepted for backward compat) - inclusion_probability = 0.5, - beta_bernoulli_alpha_between = 1, - beta_bernoulli_beta_between = 1, - dirichlet_alpha = 1, - lambda = 1, - # Priors (compare-specific) - difference_selection = TRUE, - main_difference_selection = FALSE, - difference_prior = c( - "Bernoulli", "Beta-Bernoulli", "Stochastic-Block" - ), - difference_scale = 1, - difference_prior_type = "cauchy", - difference_probability = 0.5, - # Compare difference prior hyperparameters - beta_bernoulli_alpha = 1, - beta_bernoulli_beta = 1, - difference_beta_bernoulli_alpha_between = 1, - difference_beta_bernoulli_beta_between = 1, - difference_dirichlet_alpha = 1, - difference_lambda = 1, - # Sampler - update_method = c( - "nuts", - "adaptive-metropolis", - "gibbs" - ), - target_accept = NULL, - iter = 10000L, - warmup = 1000L, - nuts_max_depth = 10L, - learn_mass_matrix = TRUE, - chains = 4L, - cores = parallel::detectCores(), - seed = NULL, - display_progress = c("per-chain", "total", "none"), - verbose = TRUE, - progress_callback = NULL) { - model_type = match.arg(model_type) - na_action = tryCatch(match.arg(na_action), error = function(e) { +bgm_spec <- function(x, + model_type = c("omrf", "ggm", "compare", "mixed_mrf"), + # Variable specification + variable_type = "ordinal", + baseline_category = NULL, + # Data (compare-specific) + y = NULL, + group_indicator = NULL, + # Missing data + na_action = c("listwise", "impute"), + # Priors (new: prior objects unpacked by bgm()) + interaction_prior_type = "cauchy", + pairwise_scale = 1, + interaction_alpha = NA_real_, + interaction_beta = NA_real_, + threshold_prior_type = "beta-prime", + main_alpha = 0.5, + main_beta = 0.5, + threshold_scale = NA_real_, + means_prior_type = "normal", + means_scale = 1, + means_alpha = NA_real_, + means_beta = NA_real_, + scale_prior_type = "gamma", + scale_shape = 1, + scale_rate = 1, + scale_eta = NA_real_, + delta = NULL, + edge_selection = TRUE, + edge_prior = bernoulli_prior(0.5), + precision_graph_prior = c("joint", "hierarchical"), + calibration_window = NULL, + # Legacy edge prior params (accepted for backward compat) + inclusion_probability = 0.5, + beta_bernoulli_alpha_between = 1, + beta_bernoulli_beta_between = 1, + dirichlet_alpha = 1, + lambda = 1, + # Priors (compare-specific) + difference_selection = TRUE, + main_difference_selection = FALSE, + difference_prior = c( + "Bernoulli", "Beta-Bernoulli", "Stochastic-Block" + ), + difference_scale = 1, + difference_prior_type = "cauchy", + difference_probability = 0.5, + # Compare difference prior hyperparameters + beta_bernoulli_alpha = 1, + beta_bernoulli_beta = 1, + difference_beta_bernoulli_alpha_between = 1, + difference_beta_bernoulli_beta_between = 1, + difference_dirichlet_alpha = 1, + difference_lambda = 1, + # Sampler + update_method = c( + "nuts", + "adaptive-metropolis", + "gibbs" + ), + target_accept = NULL, + iter = 10000L, + warmup = 1000L, + nuts_max_depth = 10L, + learn_mass_matrix = TRUE, + chains = 4L, + cores = parallel::detectCores(), + seed = NULL, + display_progress = c("per-chain", "total", "none"), + verbose = TRUE, + progress_callback = NULL) { + model_type <- match.arg(model_type) + na_action <- tryCatch(match.arg(na_action), error = function(e) { stop(paste0( "The na_action argument should be one of \"listwise\" or \"impute\", not \"", na_action, "\"." @@ -318,44 +318,44 @@ bgm_spec = function(x, }) # --- Data validation -------------------------------------------------------- - x = data_check(x, "x") - data_columnnames = if(is.null(colnames(x))) { + x <- data_check(x, "x") + data_columnnames <- if (is.null(colnames(x))) { paste0("Variable ", seq_len(ncol(x))) } else { colnames(x) } - num_variables = ncol(x) + num_variables <- ncol(x) # --- Variable types --------------------------------------------------------- - allow_continuous = (model_type != "compare") - vt = validate_variable_types( + allow_continuous <- (model_type != "compare") + vt <- validate_variable_types( variable_type = variable_type, num_variables = num_variables, allow_continuous = allow_continuous, allow_mixed = (model_type != "compare"), - caller = if(model_type == "compare") "bgmCompare" else "bgm" + caller = if (model_type == "compare") "bgmCompare" else "bgm" ) - variable_type = vt$variable_type - is_ordinal = vt$variable_bool - is_continuous = vt$is_continuous - is_mixed = vt$is_mixed + variable_type <- vt$variable_type + is_ordinal <- vt$variable_bool + is_continuous <- vt$is_continuous + is_mixed <- vt$is_mixed # Resolve model_type if "omrf" default was kept but data is continuous - if(model_type == "omrf" && is_continuous) { - model_type = "ggm" + if (model_type == "omrf" && is_continuous) { + model_type <- "ggm" } - if(model_type == "omrf" && is_mixed) { - model_type = "mixed_mrf" + if (model_type == "omrf" && is_mixed) { + model_type <- "mixed_mrf" } # Auto-resolve delta = NULL to the dimension-adaptive default 0.5 * log(p), # where p is the dimension of the continuous precision matrix. For models # without a continuous block (omrf, compare) the tilt has no target, so # NULL resolves to 0. - if(is.null(delta)) { - delta = if(model_type == "ggm") { + if (is.null(delta)) { + delta <- if (model_type == "ggm") { 0.5 * log(max(num_variables, 1)) - } else if(model_type == "mixed_mrf") { + } else if (model_type == "mixed_mrf") { 0.5 * log(max(sum(variable_type == "continuous"), 1)) } else { 0 @@ -363,11 +363,11 @@ bgm_spec = function(x, } # Validate determinant-tilt exponent and reject for pure-ordinal models - if(!is.numeric(delta) || length(delta) != 1L || is.na(delta) || + if (!is.numeric(delta) || length(delta) != 1L || is.na(delta) || !is.finite(delta) || delta < 0) { stop("'delta' must be a single finite non-negative numeric, or NULL.") } - if(delta > 0 && model_type %in% c("omrf", "compare")) { + if (delta > 0 && model_type %in% c("omrf", "compare")) { stop( "'delta' (determinant tilt) requires continuous variables; the ", "current model_type is '", model_type, "', which has no precision ", @@ -379,51 +379,51 @@ 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") { - if(!model_type %in% c("ggm", "mixed_mrf")) { + 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) { + 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." ) } - if(!edge_selection) { + 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(!interaction_prior_type %in% c("normal", "cauchy")) { + if (!interaction_prior_type %in% c("normal", "cauchy")) { stop(sprintf( paste0( - "graph_prior_spec = \"hierarchical\" supports a normal or Cauchy ", + "precision_graph_prior = \"hierarchical\" supports a normal or Cauchy ", "interaction (slab) prior: the Z-ratio normalizer is derived for ", "the normal slab, and the Cauchy slab enters through its scale ", "mixture of normals. Got %s_prior(). Use interaction_prior = ", - "normal_prior() or cauchy_prior(), or keep graph_prior_spec = ", + "normal_prior() or cauchy_prior(), or keep precision_graph_prior = ", "\"joint\"." ), interaction_prior_type )) } - if(abs(scale_shape - 1) > 1e-12) { + if (abs(scale_shape - 1) > 1e-12) { stop(sprintf( paste0( - "graph_prior_spec = \"hierarchical\" requires shape = 1 on the ", + "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 ", - "graph_prior_spec = \"joint\"." + "precision_graph_prior = \"joint\"." ), format(scale_shape) )) @@ -431,7 +431,7 @@ bgm_spec = function(x, } # --- Sampler (needs is_continuous and edge_selection early) ------------------ - sampler = validate_sampler( + sampler <- validate_sampler( update_method = update_method, target_accept = target_accept, iter = iter, @@ -443,20 +443,20 @@ bgm_spec = function(x, seed = seed, display_progress = display_progress, is_continuous = is_continuous, - edge_selection = if(model_type == "compare") FALSE else edge_selection, + edge_selection = if (model_type == "compare") FALSE else edge_selection, verbose = verbose, progress_callback = progress_callback ) # --- Resolve edge prior object ----------------------------------------------- - if(inherits(edge_prior, "bgms_indicator_prior")) { - ep_flat = unpack_indicator_prior(edge_prior, num_variables) - } else if(is.character(edge_prior)) { + if (inherits(edge_prior, "bgms_indicator_prior")) { + ep_flat <- unpack_indicator_prior(edge_prior, num_variables) + } else if (is.character(edge_prior)) { # Legacy string path (tests and bgmCompare may call bgm_spec directly) - edge_prior_str = match.arg(edge_prior, + edge_prior_str <- match.arg(edge_prior, choices = c("Bernoulli", "Beta-Bernoulli", "Stochastic-Block") ) - ep_flat = validate_edge_prior( + ep_flat <- validate_edge_prior( edge_selection = edge_selection, edge_prior = edge_prior_str, inclusion_probability = inclusion_probability, num_variables = num_variables, @@ -466,12 +466,12 @@ bgm_spec = function(x, beta_bernoulli_beta_between = beta_bernoulli_beta_between, dirichlet_alpha = dirichlet_alpha, lambda = lambda ) - ep_flat$beta_bernoulli_alpha = beta_bernoulli_alpha - ep_flat$beta_bernoulli_beta = beta_bernoulli_beta - ep_flat$beta_bernoulli_alpha_between = beta_bernoulli_alpha_between - ep_flat$beta_bernoulli_beta_between = beta_bernoulli_beta_between - ep_flat$dirichlet_alpha = dirichlet_alpha - ep_flat$lambda = lambda + ep_flat$beta_bernoulli_alpha <- beta_bernoulli_alpha + ep_flat$beta_bernoulli_beta <- beta_bernoulli_beta + ep_flat$beta_bernoulli_alpha_between <- beta_bernoulli_alpha_between + ep_flat$beta_bernoulli_beta_between <- beta_bernoulli_beta_between + ep_flat$dirichlet_alpha <- dirichlet_alpha + ep_flat$lambda <- lambda } else { stop( "'edge_prior' must be a bgms_indicator_prior object.", @@ -479,15 +479,15 @@ bgm_spec = function(x, ) } # Override edge_selection if explicitly FALSE - if(!edge_selection) { - ep_flat$edge_selection = FALSE - ep_flat$edge_prior = "Not Applicable" - ep_flat$inclusion_probability = matrix(0.5, nrow = 1, ncol = 1) + if (!edge_selection) { + ep_flat$edge_selection <- FALSE + ep_flat$edge_prior <- "Not Applicable" + ep_flat$inclusion_probability <- matrix(0.5, nrow = 1, ncol = 1) } # --- Build by model type ---------------------------------------------------- - if(model_type == "ggm") { - spec = build_spec_ggm( + if (model_type == "ggm") { + spec <- build_spec_ggm( x = x, data_columnnames = data_columnnames, num_variables = num_variables, variable_type = variable_type, is_ordinal = is_ordinal, @@ -503,12 +503,12 @@ 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 ) - } else if(model_type == "mixed_mrf") { - spec = build_spec_mixed_mrf( + } else if (model_type == "mixed_mrf") { + spec <- build_spec_mixed_mrf( x = x, data_columnnames = data_columnnames, num_variables = num_variables, variable_type = variable_type, is_ordinal = is_ordinal, @@ -530,12 +530,12 @@ 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 ) - } else if(model_type == "omrf") { - spec = build_spec_omrf( + } else if (model_type == "omrf") { + spec <- build_spec_omrf( x = x, data_columnnames = data_columnnames, num_variables = num_variables, variable_type = variable_type, is_ordinal = is_ordinal, @@ -552,7 +552,7 @@ bgm_spec = function(x, edge_prior_flat = ep_flat ) } else { - spec = build_spec_compare( + spec <- build_spec_compare( x = x, y = y, group_indicator = group_indicator, data_columnnames = data_columnnames, num_variables = num_variables, diff --git a/R/build_arguments.R b/R/build_arguments.R index ba473000..94b7c7f7 100644 --- a/R/build_arguments.R +++ b/R/build_arguments.R @@ -28,35 +28,35 @@ build_arguments = function(spec) { -build_arguments_ggm = 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" ) } @@ -103,49 +103,49 @@ build_arguments_omrf = function(spec) { } -build_arguments_mixed_mrf = 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..94674d28 100644 --- a/R/build_output_bgm.R +++ b/R/build_output_bgm.R @@ -11,22 +11,22 @@ # 1. Parameter naming: GGM uses "Var (precision)", OMRF uses "Var (k)" # 2. Main posterior mean shape: GGM = px1, OMRF = pxmax_categories # ============================================================================== -build_output_bgm = function(spec, raw) { - d = spec$data - v = spec$variables - p = spec$prior - s = spec$sampler +build_output_bgm <- function(spec, raw) { + d <- spec$data + v <- spec$variables + p <- spec$prior + s <- spec$sampler - is_continuous = v$is_continuous - num_variables = d$num_variables - data_columnnames = d$data_columnnames - edge_selection = p$edge_selection - edge_prior = p$edge_prior + is_continuous <- v$is_continuous + num_variables <- d$num_variables + data_columnnames <- d$data_columnnames + edge_selection <- p$edge_selection + edge_prior <- p$edge_prior # 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 @@ -37,42 +37,42 @@ build_output_bgm = function(spec, raw) { # via convert_results_to_list(). Split into main and pairwise components and # transpose to (iters x params) --- the layout that MCMC summary functions # expect. - if(is_continuous) { + if (is_continuous) { # GGM: samples contain the upper triangle of the precision matrix # (row-major). Diagonal entries are "main"; off-diagonal are "pairwise". - diag_idx = integer(num_variables) - offdiag_idx = integer(num_variables * (num_variables - 1L) / 2L) - pos = 0L - di = 0L - oi = 0L - for(i in seq_len(num_variables)) { - for(j in i:num_variables) { - pos = pos + 1L - if(i == j) { - di = di + 1L - diag_idx[di] = pos + diag_idx <- integer(num_variables) + offdiag_idx <- integer(num_variables * (num_variables - 1L) / 2L) + pos <- 0L + di <- 0L + oi <- 0L + for (i in seq_len(num_variables)) { + for (j in i:num_variables) { + pos <- pos + 1L + if (i == j) { + di <- di + 1L + diag_idx[di] <- pos } else { - oi = oi + 1L - offdiag_idx[oi] = pos + oi <- oi + 1L + offdiag_idx[oi] <- pos } } } - raw = lapply(raw, function(chain) { - samples_t = t(chain$samples) - res = list( + raw <- lapply(raw, function(chain) { + samples_t <- t(chain$samples) + res <- list( main_samples = samples_t[, diag_idx, drop = FALSE], pairwise_samples = samples_t[, offdiag_idx, drop = FALSE], userInterrupt = isTRUE(chain$userInterrupt), chain_id = chain$chain_id ) - if(!is.null(chain$indicator_samples)) { - res$indicator_samples = t(chain$indicator_samples)[, offdiag_idx, drop = FALSE] + if (!is.null(chain$indicator_samples)) { + res$indicator_samples <- t(chain$indicator_samples)[, offdiag_idx, drop = FALSE] } - if(!is.null(chain$allocation_samples)) { - res$allocations = t(chain$allocation_samples) + if (!is.null(chain$allocation_samples)) { + res$allocations <- t(chain$allocation_samples) } - if(!is.null(chain$inclusion_parameter_samples)) { - res$inclusion_parameter = as.numeric(chain$inclusion_parameter_samples) + if (!is.null(chain$inclusion_parameter_samples)) { + res$inclusion_parameter <- as.numeric(chain$inclusion_parameter_samples) } attach_diagnostic_traces(res, chain) }) @@ -80,52 +80,52 @@ build_output_bgm = function(spec, raw) { # OMRF: the first num_thresholds params are main effects, the rest are # pairwise. NUTS diagnostics use bare names from C++; rename to the # trailing-__ convention expected by summarize_nuts_diagnostics(). - num_thresholds = spec$precomputed$num_thresholds - raw = lapply(raw, function(chain) { - samples_t = t(chain$samples) - n_params = ncol(samples_t) - res = list( + num_thresholds <- spec$precomputed$num_thresholds + raw <- lapply(raw, function(chain) { + samples_t <- t(chain$samples) + n_params <- ncol(samples_t) + res <- list( main_samples = samples_t[, seq_len(num_thresholds), drop = FALSE], pairwise_samples = samples_t[, seq(num_thresholds + 1L, n_params), drop = FALSE], userInterrupt = isTRUE(chain$userInterrupt), chain_id = chain$chain_id ) - if(!is.null(chain$indicator_samples)) { - res$indicator_samples = t(chain$indicator_samples) + if (!is.null(chain$indicator_samples)) { + res$indicator_samples <- t(chain$indicator_samples) } - if(!is.null(chain$allocation_samples)) { - res$allocations = t(chain$allocation_samples) + if (!is.null(chain$allocation_samples)) { + res$allocations <- t(chain$allocation_samples) } - if(!is.null(chain$inclusion_parameter_samples)) { - res$inclusion_parameter = as.numeric(chain$inclusion_parameter_samples) + if (!is.null(chain$inclusion_parameter_samples)) { + res$inclusion_parameter <- as.numeric(chain$inclusion_parameter_samples) } attach_diagnostic_traces(res, chain) }) } # --- Parameter names -------------------------------------------------------- - if(is_continuous) { + if (is_continuous) { # GGM: raw samples are on precision scale; label accordingly - names_main = paste0(data_columnnames, " (precision)") - is_ordinal_variable = NULL - num_categories = NULL + names_main <- paste0(data_columnnames, " (precision)") + is_ordinal_variable <- NULL + num_categories <- NULL } else { # OMRF: per-category thresholds or BC linear/quadratic - is_ordinal_variable = v$is_ordinal - num_categories = d$num_categories - names_main = character() - for(vi in seq_len(num_variables)) { - if(is_ordinal_variable[vi]) { - cats = ordinal_threshold_labels( + is_ordinal_variable <- v$is_ordinal + num_categories <- d$num_categories + names_main <- character() + for (vi in seq_len(num_variables)) { + if (is_ordinal_variable[vi]) { + cats <- ordinal_threshold_labels( num_categories[vi], - if(!is.null(d$category_levels)) d$category_levels[[vi]] else NULL + if (!is.null(d$category_levels)) d$category_levels[[vi]] else NULL ) - names_main = c( + names_main <- c( names_main, paste0(data_columnnames[vi], " (", cats, ")") ) } else { - names_main = c( + names_main <- c( names_main, paste0(data_columnnames[vi], " (linear)"), paste0(data_columnnames[vi], " (quadratic)") @@ -134,10 +134,10 @@ build_output_bgm = function(spec, raw) { } } - edge_names = character() - for(i in seq_len(num_variables - 1)) { - for(j in seq(i + 1, num_variables)) { - edge_names = c( + edge_names <- character() + for (i in seq_len(num_variables - 1)) { + for (j in seq(i + 1, num_variables)) { + edge_names <- c( edge_names, paste0(data_columnnames[i], "-", data_columnnames[j]) ) @@ -148,129 +148,129 @@ build_output_bgm = function(spec, raw) { # Store normalized raw chains and metadata so that ensure_summaries() can # compute ESS/Rhat/MCSE on first demand. Posterior_summary_* fields start # as NULL and are populated lazily. - cache = new.env(parent = emptyenv()) - cache$raw = raw - cache$edge_selection = edge_selection - cache$names_main = names_main - cache$edge_names = edge_names - cache$is_continuous = is_continuous - cache$model_type = if(is_continuous) "ggm" else "omrf" - cache$summaries_computed = FALSE + cache <- new.env(parent = emptyenv()) + cache$raw <- raw + cache$edge_selection <- edge_selection + cache$names_main <- names_main + cache$edge_names <- edge_names + cache$is_continuous <- is_continuous + cache$model_type <- if (is_continuous) "ggm" else "omrf" + cache$summaries_computed <- FALSE # --- Compute posterior means from raw samples (cheap) ----------------------- - pooled_main = do.call(rbind, lapply(raw, function(ch) ch$main_samples)) - pooled_pair = do.call(rbind, lapply(raw, function(ch) ch$pairwise_samples)) - main_means = colMeans(pooled_main) - pair_means = colMeans(pooled_pair) + pooled_main <- do.call(rbind, lapply(raw, function(ch) ch$main_samples)) + pooled_pair <- do.call(rbind, lapply(raw, function(ch) ch$pairwise_samples)) + main_means <- colMeans(pooled_main) + pair_means <- colMeans(pooled_pair) - results = list() + results <- list() # --- Edge selection summaries (deferred) ------------------------------------ - has_sbm = FALSE - if(edge_selection) { - has_sbm = identical(edge_prior, "Stochastic-Block") && + has_sbm <- FALSE + if (edge_selection) { + has_sbm <- identical(edge_prior, "Stochastic-Block") && "allocations" %in% names(raw[[1]]) - if(has_sbm) { - sbm_convergence = summarize_alloc_pairs( + if (has_sbm) { + sbm_convergence <- summarize_alloc_pairs( allocations = lapply(raw, `[[`, "allocations"), node_names = data_columnnames ) - results$posterior_summary_pairwise_allocations = sbm_convergence$sbm_summary - co_occur_matrix = sbm_convergence$co_occur_matrix + results$posterior_summary_pairwise_allocations <- sbm_convergence$sbm_summary + co_occur_matrix <- sbm_convergence$co_occur_matrix } - if("inclusion_parameter" %in% names(raw[[1]])) { - results$inclusion_parameter_samples = + if ("inclusion_parameter" %in% names(raw[[1]])) { + results$inclusion_parameter_samples <- lapply(raw, `[[`, "inclusion_parameter") } } # --- Posterior mean: main --------------------------------------------------- - if(is_continuous) { + if (is_continuous) { # GGM has no main effects - results$posterior_mean_main = NULL + results$posterior_mean_main <- NULL } else { # OMRF: p x max_categories matrix - num_params = ifelse(is_ordinal_variable, num_categories, 2L) - max_num_categories = max(num_params) + num_params <- ifelse(is_ordinal_variable, num_categories, 2L) + max_num_categories <- max(num_params) - pmm = matrix(NA, nrow = num_variables, ncol = max_num_categories) - start = 0L - stop = 0L - for(vi in seq_len(num_variables)) { - if(is_ordinal_variable[vi]) { - start = stop + 1L - stop = start + num_categories[vi] - 1L - pmm[vi, seq_len(num_categories[vi])] = main_means[start:stop] + pmm <- matrix(NA, nrow = num_variables, ncol = max_num_categories) + start <- 0L + stop <- 0L + for (vi in seq_len(num_variables)) { + if (is_ordinal_variable[vi]) { + start <- stop + 1L + stop <- start + num_categories[vi] - 1L + pmm[vi, seq_len(num_categories[vi])] <- main_means[start:stop] } else { - start = stop + 1L - stop = start + 1L - pmm[vi, 1:2] = main_means[start:stop] + start <- stop + 1L + stop <- start + 1L + pmm[vi, 1:2] <- main_means[start:stop] } } - results$posterior_mean_main = pmm - rownames(results$posterior_mean_main) = data_columnnames - colnames(results$posterior_mean_main) = paste0("cat (", seq_len(ncol(pmm)), ")") + results$posterior_mean_main <- pmm + rownames(results$posterior_mean_main) <- data_columnnames + colnames(results$posterior_mean_main) <- paste0("cat (", seq_len(ncol(pmm)), ")") } # --- Posterior mean: associations ------------------------------------------- # For GGM: C++ stores precision; convert to association scale (* -0.5). # For OMRF: C++ already stores association-scale values. - associations = matrix(0, + associations <- matrix(0, nrow = num_variables, ncol = num_variables, dimnames = list(data_columnnames, data_columnnames) ) - associations[lower.tri(associations)] = pair_means - associations = associations + t(associations) - if(is_continuous) { - associations = -0.5 * associations + associations[lower.tri(associations)] <- pair_means + associations <- associations + t(associations) + if (is_continuous) { + associations <- -0.5 * associations } - results$posterior_mean_pairwise = associations + results$posterior_mean_pairwise <- associations # --- Residual variance (GGM only) ------------------------------------------- # C++ stores precision diagonal; residual variance = 1 / precision. # Average per-sample inversions to avoid Jensen's inequality bias. - if(is_continuous) { - results$posterior_mean_residual_variance = colMeans(1 / pooled_main) - names(results$posterior_mean_residual_variance) = data_columnnames + if (is_continuous) { + results$posterior_mean_residual_variance <- colMeans(1 / pooled_main) + names(results$posterior_mean_residual_variance) <- data_columnnames } # --- Posterior mean: indicator + SBM ---------------------------------------- - if(edge_selection) { - pooled_ind = do.call(rbind, lapply(raw, function(ch) ch$indicator_samples)) - indicator_means = colMeans(pooled_ind) - results$posterior_mean_indicator = matrix(0, + if (edge_selection) { + pooled_ind <- do.call(rbind, lapply(raw, function(ch) ch$indicator_samples)) + indicator_means <- colMeans(pooled_ind) + results$posterior_mean_indicator <- matrix(0, nrow = num_variables, ncol = num_variables, dimnames = list(data_columnnames, data_columnnames) ) - results$posterior_mean_indicator[lower.tri(results$posterior_mean_indicator)] = + results$posterior_mean_indicator[lower.tri(results$posterior_mean_indicator)] <- indicator_means - results$posterior_mean_indicator = results$posterior_mean_indicator + + results$posterior_mean_indicator <- results$posterior_mean_indicator + t(results$posterior_mean_indicator) - if(has_sbm) { - results$posterior_mean_coclustering_matrix = co_occur_matrix - results = attach_sbm_posterior_summary(results, raw, build_arguments(spec)) + if (has_sbm) { + results$posterior_mean_coclustering_matrix <- co_occur_matrix + results <- attach_sbm_posterior_summary(results, raw, build_arguments(spec)) } } # --- arguments + class ------------------------------------------------------ - results$arguments = build_arguments(spec) + results$arguments <- build_arguments(spec) # Report the number of chains actually kept; failed chains are dropped # upstream, so raw holds only the survivors. - results$arguments$num_chains = length(raw) - class(results) = "bgms" + results$arguments$num_chains <- length(raw) + class(results) <- "bgms" # --- raw_samples ------------------------------------------------------------ # Allocation samples are per node (one column per variable), so they carry # node names, not edge names. - alloc_names = if(identical(edge_prior, "Stochastic-Block")) { + alloc_names <- if (identical(edge_prior, "Stochastic-Block")) { data_columnnames } else { NULL } - results$raw_samples = build_raw_samples_list( + results$raw_samples <- build_raw_samples_list( raw, edge_selection, edge_prior, names_main, edge_names, alloc_names ) @@ -278,28 +278,28 @@ build_output_bgm = function(spec, raw) { # NULL placeholders ensure names(fit) lists these fields for easybgm compat. # The $.bgms method routes actual access through the cache. # Use list(NULL) because results$x = NULL removes the element in R. - results["posterior_summary_main"] = list(NULL) - results["posterior_summary_pairwise"] = list(NULL) - if(edge_selection) { - results["posterior_summary_indicator"] = list(NULL) + results["posterior_summary_main"] <- list(NULL) + results["posterior_summary_pairwise"] <- list(NULL) + if (edge_selection) { + results["posterior_summary_indicator"] <- list(NULL) } - results$cache = cache + results$cache <- cache # --- Sampler diagnostics ---------------------------------------------------- - results = attach_sampler_diagnostics( + results <- attach_sampler_diagnostics( results, raw, s$update_method, s$nuts_max_depth, names_main = names_main, names_pairwise = edge_names, target_accept = s$target_accept ) # --- Z-ratio alarm suite (hierarchical graph-prior spec) ---------------------- - if(!is.null(zratio_chains)) { - zc = zratio_constants( + if (!is.null(zratio_chains)) { + zc <- zratio_constants( delta = p$delta, sigma = 2 * p$pairwise_scale, beta = p$scale_rate / 2 ) - results$zratio_diag = summarize_zratio_diagnostics( + results$zratio_diag <- summarize_zratio_diagnostics( zratio_chains, zc, num_nodes = num_variables, seed = s$seed, @@ -307,8 +307,8 @@ build_output_bgm = function(spec, raw) { ) } - results$.bgm_spec = spec - if(needs_easybgm_s3_compat()) { + results$.bgm_spec <- spec + if (needs_easybgm_s3_compat()) { results } else { s3_list_to_bgms(results) diff --git a/R/build_output_mixed_mrf.R b/R/build_output_mixed_mrf.R index 35439ad9..29a24330 100644 --- a/R/build_output_mixed_mrf.R +++ b/R/build_output_mixed_mrf.R @@ -17,42 +17,42 @@ # stored separately in posterior_mean_residual_variance; # the diagonal of the pairwise interaction matrix is zero. # ============================================================================== -build_output_mixed_mrf = function(spec, raw) { - d = spec$data - v = spec$variables - pr = spec$prior - s = spec$sampler +build_output_mixed_mrf <- function(spec, raw) { + d <- spec$data + v <- spec$variables + pr <- spec$prior + s <- spec$sampler - p = d$num_discrete - q = d$num_continuous - num_variables = d$num_variables - data_columnnames = d$data_columnnames - disc_names = d$data_columnnames_discrete - cont_names = d$data_columnnames_continuous - disc_idx = d$discrete_indices - cont_idx = d$continuous_indices - is_ordinal = v$is_ordinal - num_categories = d$num_categories - edge_selection = pr$edge_selection + p <- d$num_discrete + q <- d$num_continuous + num_variables <- d$num_variables + data_columnnames <- d$data_columnnames + disc_names <- d$data_columnnames_discrete + cont_names <- d$data_columnnames_continuous + disc_idx <- d$discrete_indices + cont_idx <- d$continuous_indices + is_ordinal <- v$is_ordinal + num_categories <- d$num_categories + edge_selection <- pr$edge_selection # 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 } # --- Compute index layout in flat parameter vector -------------------------- - layout = compute_mixed_parameter_indices( + layout <- compute_mixed_parameter_indices( num_thresholds = spec$precomputed$num_thresholds, p = p, q = q ) - nt = layout$num_thresholds - main_idx = layout$main_idx - pairwise_idx = layout$pairwise_idx + nt <- layout$num_thresholds + main_idx <- layout$main_idx + pairwise_idx <- layout$pairwise_idx # --- Indicator index layout ------------------------------------------------- # C++ indicator vector: [Gxx_ut | Gyy_ut | Gxy] @@ -60,63 +60,63 @@ build_output_mixed_mrf = function(spec, raw) { # Discrete, continuous, cross edges --- same order as pairwise_idx above. # --- Normalize raw output per chain ----------------------------------------- - raw = lapply(raw, function(chain) { - samples_t = t(chain$samples) - res = list( + raw <- lapply(raw, function(chain) { + samples_t <- t(chain$samples) + res <- list( main_samples = samples_t[, main_idx, drop = FALSE], pairwise_samples = samples_t[, pairwise_idx, drop = FALSE], userInterrupt = isTRUE(chain$userInterrupt), chain_id = chain$chain_id ) - if(!is.null(chain$indicator_samples)) { - res$indicator_samples = t(chain$indicator_samples) + if (!is.null(chain$indicator_samples)) { + res$indicator_samples <- t(chain$indicator_samples) } - if(!is.null(chain$allocation_samples)) { - res$allocations = t(chain$allocation_samples) + if (!is.null(chain$allocation_samples)) { + res$allocations <- t(chain$allocation_samples) } - if(!is.null(chain$inclusion_parameter_samples)) { - res$inclusion_parameter = as.numeric(chain$inclusion_parameter_samples) + if (!is.null(chain$inclusion_parameter_samples)) { + res$inclusion_parameter <- as.numeric(chain$inclusion_parameter_samples) } attach_diagnostic_traces(res, chain) }) # --- Parameter names -------------------------------------------------------- # Main effect names (in internal order: discrete first, continuous second) - names_main = character() - for(si in seq_len(p)) { - if(is_ordinal[si]) { - cats = ordinal_threshold_labels( + names_main <- character() + for (si in seq_len(p)) { + if (is_ordinal[si]) { + cats <- ordinal_threshold_labels( num_categories[si], - if(!is.null(d$category_levels)) d$category_levels[[si]] else NULL + if (!is.null(d$category_levels)) d$category_levels[[si]] else NULL ) - names_main = c(names_main, paste0(disc_names[si], " (", cats, ")")) + names_main <- c(names_main, paste0(disc_names[si], " (", cats, ")")) } else { - names_main = c( + names_main <- c( names_main, paste0(disc_names[si], " (linear)"), paste0(disc_names[si], " (quadratic)") ) } } - for(ji in seq_len(q)) { - names_main = c(names_main, paste0(cont_names[ji], " (mean)")) + for (ji in seq_len(q)) { + names_main <- c(names_main, paste0(cont_names[ji], " (mean)")) } - for(ji in seq_len(q)) { - names_main = c(names_main, paste0(cont_names[ji], " (precision diag)")) + for (ji in seq_len(q)) { + names_main <- c(names_main, paste0(cont_names[ji], " (precision diag)")) } # Pairwise edge names --- internal order, mapped to original column names # We need a mapping from internal index to original variable name # Internal variables: [disc_1, ..., disc_p, cont_1, ..., cont_q] # Their original names: c(disc_names, cont_names) - all_internal_names = c(disc_names, cont_names) + all_internal_names <- c(disc_names, cont_names) - edge_names = character() + edge_names <- character() # Discrete-discrete edges - if(p > 1) { - for(i in seq_len(p - 1)) { - for(j in seq(i + 1, p)) { - edge_names = c( + if (p > 1) { + for (i in seq_len(p - 1)) { + for (j in seq(i + 1, p)) { + edge_names <- c( edge_names, paste0(disc_names[i], "-", disc_names[j]) ) @@ -124,10 +124,10 @@ build_output_mixed_mrf = function(spec, raw) { } } # Continuous-continuous edges (off-diagonal) - if(q > 1) { - for(i in seq_len(q - 1)) { - for(j in seq(i + 1, q)) { - edge_names = c( + if (q > 1) { + for (i in seq_len(q - 1)) { + for (j in seq(i + 1, q)) { + edge_names <- c( edge_names, paste0(cont_names[i], "-", cont_names[j]) ) @@ -135,10 +135,10 @@ build_output_mixed_mrf = function(spec, raw) { } } # Cross edges (discrete-continuous) - if(p > 0 && q > 0) { - for(i in seq_len(p)) { - for(j in seq_len(q)) { - edge_names = c( + if (p > 0 && q > 0) { + for (i in seq_len(p)) { + for (j in seq_len(q)) { + edge_names <- c( edge_names, paste0(disc_names[i], "-", cont_names[j]) ) @@ -147,139 +147,139 @@ build_output_mixed_mrf = function(spec, raw) { } # --- Lazy MCMC diagnostics cache -------------------------------------------- - cache = new.env(parent = emptyenv()) - cache$raw = raw - cache$edge_selection = edge_selection - cache$names_main = names_main - cache$edge_names = edge_names - cache$is_continuous = FALSE - cache$model_type = "mixed_mrf" - cache$summaries_computed = FALSE + cache <- new.env(parent = emptyenv()) + cache$raw <- raw + cache$edge_selection <- edge_selection + cache$names_main <- names_main + cache$edge_names <- edge_names + cache$is_continuous <- FALSE + cache$model_type <- "mixed_mrf" + cache$summaries_computed <- FALSE # --- Compute posterior means from raw samples (cheap) ----------------------- - pooled_main = do.call(rbind, lapply(raw, function(ch) ch$main_samples)) - pooled_pair = do.call(rbind, lapply(raw, function(ch) ch$pairwise_samples)) - main_means = colMeans(pooled_main) - pair_means = colMeans(pooled_pair) + pooled_main <- do.call(rbind, lapply(raw, function(ch) ch$main_samples)) + pooled_pair <- do.call(rbind, lapply(raw, function(ch) ch$pairwise_samples)) + main_means <- colMeans(pooled_main) + pair_means <- colMeans(pooled_pair) # Split main_means into true main effects and quadratic (precision diagonal) - n_main = nt + q # thresholds + continuous means - n_quad = layout$num_quadratic # precision diagonal entries + n_main <- nt + q # thresholds + continuous means + n_quad <- layout$num_quadratic # precision diagonal entries - cache$n_main = n_main - cache$n_quad = n_quad + cache$n_main <- n_main + cache$n_quad <- n_quad - results = list() + results <- list() # --- Edge selection summaries (deferred) ------------------------------------ - edge_prior = pr$edge_prior - has_sbm = FALSE - if(edge_selection) { - has_sbm = identical(edge_prior, "Stochastic-Block") && + edge_prior <- pr$edge_prior + has_sbm <- FALSE + if (edge_selection) { + has_sbm <- identical(edge_prior, "Stochastic-Block") && "allocations" %in% names(raw[[1]]) - if(has_sbm) { - sbm_convergence = summarize_alloc_pairs( + if (has_sbm) { + sbm_convergence <- summarize_alloc_pairs( allocations = lapply(raw, `[[`, "allocations"), node_names = all_internal_names ) - results$posterior_summary_pairwise_allocations = sbm_convergence$sbm_summary - co_occur_matrix = sbm_convergence$co_occur_matrix + results$posterior_summary_pairwise_allocations <- sbm_convergence$sbm_summary + co_occur_matrix <- sbm_convergence$co_occur_matrix } - if("inclusion_parameter" %in% names(raw[[1]])) { - results$inclusion_parameter_samples = + if ("inclusion_parameter" %in% names(raw[[1]])) { + results$inclusion_parameter_samples <- lapply(raw, `[[`, "inclusion_parameter") } } # --- Posterior mean: main --------------------------------------------------- # Discrete main effects: p x max_cats matrix (like OMRF) - num_params_disc = ifelse(is_ordinal, num_categories, 2L) - max_num_cats = max(num_params_disc) - pmm_disc = matrix(NA, nrow = p, ncol = max_num_cats) - start = 0L - stop = 0L - for(si in seq_len(p)) { - if(is_ordinal[si]) { - start = stop + 1L - stop = start + num_categories[si] - 1L - pmm_disc[si, seq_len(num_categories[si])] = main_means[start:stop] + num_params_disc <- ifelse(is_ordinal, num_categories, 2L) + max_num_cats <- max(num_params_disc) + pmm_disc <- matrix(NA, nrow = p, ncol = max_num_cats) + start <- 0L + stop <- 0L + for (si in seq_len(p)) { + if (is_ordinal[si]) { + start <- stop + 1L + stop <- start + num_categories[si] - 1L + pmm_disc[si, seq_len(num_categories[si])] <- main_means[start:stop] } else { - start = stop + 1L - stop = start + 1L - pmm_disc[si, 1:2] = main_means[start:stop] + start <- stop + 1L + stop <- start + 1L + pmm_disc[si, 1:2] <- main_means[start:stop] } } - rownames(pmm_disc) = disc_names - colnames(pmm_disc) = paste0("cat (", seq_len(max_num_cats), ")") + rownames(pmm_disc) <- disc_names + colnames(pmm_disc) <- paste0("cat (", seq_len(max_num_cats), ")") # Continuous main effects: q x 1 matrix (means only) - pmm_cont = matrix(main_means[nt + seq_len(q)], + pmm_cont <- matrix(main_means[nt + seq_len(q)], nrow = q, ncol = 1, dimnames = list(cont_names, "mean") ) - results$posterior_mean_main = list( + results$posterior_mean_main <- list( discrete = pmm_disc, continuous = pmm_cont ) # --- Posterior mean: associations (all blocks) ----------------------------- - dn = list(data_columnnames, data_columnnames) - results$posterior_mean_pairwise = fill_mixed_symmetric( + dn <- list(data_columnnames, data_columnnames) + results$posterior_mean_pairwise <- fill_mixed_symmetric( pair_means, p, q, disc_idx, cont_idx, dn ) # --- Residual variance (continuous diagonal) -------------------------------- # C++ stores negative association diagonal; residual variance = -1/(2*diag). # Average per-sample inversions to avoid Jensen's inequality bias. - pooled_quad = pooled_main[, nt + q + seq_len(q), drop = FALSE] - rv = colMeans(-1 / (2 * pooled_quad)) - names(rv) = cont_names - results$posterior_mean_residual_variance = rv + pooled_quad <- pooled_main[, nt + q + seq_len(q), drop = FALSE] + rv <- colMeans(-1 / (2 * pooled_quad)) + names(rv) <- cont_names + results$posterior_mean_residual_variance <- rv # --- Posterior mean: indicator ----------------------------------------------- - if(edge_selection) { - pooled_ind = do.call(rbind, lapply(raw, function(ch) ch$indicator_samples)) - indicator_means = colMeans(pooled_ind) - results$posterior_mean_indicator = fill_mixed_symmetric( + if (edge_selection) { + pooled_ind <- do.call(rbind, lapply(raw, function(ch) ch$indicator_samples)) + indicator_means <- colMeans(pooled_ind) + results$posterior_mean_indicator <- fill_mixed_symmetric( indicator_means, p, q, disc_idx, cont_idx, dn ) - if(has_sbm) { - results$posterior_mean_coclustering_matrix = co_occur_matrix - results = attach_sbm_posterior_summary(results, raw, build_arguments(spec)) + if (has_sbm) { + results$posterior_mean_coclustering_matrix <- co_occur_matrix + results <- attach_sbm_posterior_summary(results, raw, build_arguments(spec)) } } # --- arguments + class ------------------------------------------------------ - results$arguments = build_arguments(spec) + results$arguments <- build_arguments(spec) # Report the number of chains actually kept; failed chains are dropped # upstream, so raw holds only the survivors. - results$arguments$num_chains = length(raw) + results$arguments$num_chains <- length(raw) # NULL placeholders ensure names(fit) lists these fields for easybgm compat. # Use list(NULL) because results$x = NULL removes the element in R. - results["posterior_summary_main"] = list(NULL) - results["posterior_summary_pairwise"] = list(NULL) - if(edge_selection) { - results["posterior_summary_indicator"] = list(NULL) + results["posterior_summary_main"] <- list(NULL) + results["posterior_summary_pairwise"] <- list(NULL) + if (edge_selection) { + results["posterior_summary_indicator"] <- list(NULL) } - results$cache = cache - class(results) = "bgms" + results$cache <- cache + class(results) <- "bgms" # --- raw_samples ------------------------------------------------------------ - alloc_names = if(identical(edge_prior, "Stochastic-Block")) { + alloc_names <- if (identical(edge_prior, "Stochastic-Block")) { all_internal_names } else { NULL } - results$raw_samples = build_raw_samples_list( + results$raw_samples <- build_raw_samples_list( raw, edge_selection, edge_prior, names_main, edge_names, alloc_names ) # --- Sampler diagnostics ---------------------------------------------------- - results = attach_sampler_diagnostics( + results <- attach_sampler_diagnostics( results, raw, s$update_method, s$nuts_max_depth, names_main = names_main, names_pairwise = edge_names, target_accept = s$target_accept @@ -288,13 +288,13 @@ build_output_mixed_mrf = function(spec, raw) { # --- Z-ratio alarm suite (hierarchical spec on the continuous block) --------- # The indicator vector is [Gxx upper, Gyy upper, Gxy row-major], both # triangles without diagonals; the audit reads the Gyy segment. - if(!is.null(zratio_chains)) { - zc = zratio_constants( + if (!is.null(zratio_chains)) { + zc <- zratio_constants( delta = pr$delta, sigma = 2 * pr$pairwise_scale, beta = pr$scale_rate / 2 ) - results$zratio_diag = summarize_zratio_diagnostics( + results$zratio_diag <- summarize_zratio_diagnostics( zratio_chains, zc, num_nodes = q, seed = s$seed, @@ -304,8 +304,8 @@ build_output_mixed_mrf = function(spec, raw) { ) } - results$.bgm_spec = spec - if(needs_easybgm_s3_compat()) { + results$.bgm_spec <- spec + if (needs_easybgm_s3_compat()) { results } else { s3_list_to_bgms(results) diff --git a/R/build_spec.R b/R/build_spec.R index b89e4a19..e986459c 100644 --- a/R/build_spec.R +++ b/R/build_spec.R @@ -52,34 +52,34 @@ sampler_sublist = function(s) { -build_spec_ggm = function(x, data_columnnames, num_variables, - variable_type, is_ordinal, is_continuous, - baseline_category, - na_action, sampler, - interaction_prior_type, pairwise_scale, - interaction_alpha, interaction_beta, - scale_prior_type, scale_shape, scale_rate, - scale_eta = NA_real_, - delta = 0, - graph_prior_spec = "joint", - calibration_window = NULL, - edge_prior_flat) { +build_spec_ggm <- function(x, data_columnnames, num_variables, + variable_type, is_ordinal, is_continuous, + baseline_category, + na_action, sampler, + interaction_prior_type, pairwise_scale, + interaction_alpha, interaction_beta, + scale_prior_type, scale_shape, scale_rate, + scale_eta = NA_real_, + delta = 0, + precision_graph_prior = "joint", + calibration_window = NULL, + edge_prior_flat) { # Missing data - md = validate_missing_data( + md <- validate_missing_data( x = x, na_action = na_action, is_continuous = TRUE ) - x = md$x + x <- md$x # Center continuous data (GGM likelihood assumes zero mean). The column # means are kept so prediction can center newdata on the training scale. - column_means = colMeans(x) - x = center_continuous_data(x) + column_means <- colMeans(x) + x <- center_continuous_data(x) # Standardized-frame scale prior: derive the raw diagonal rate eta / s - scale_rate = resolve_scale_rate(scale_rate, scale_eta, pairwise_scale) + scale_rate <- resolve_scale_rate(scale_rate, scale_eta, pairwise_scale) - ep = edge_prior_flat + ep <- edge_prior_flat new_bgm_spec( model_type = "ggm", @@ -112,7 +112,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) @@ -215,53 +215,53 @@ build_spec_omrf = function(x, data_columnnames, num_variables, # the spec with metadata needed by sample_mixed_mrf() and # build_output_mixed_mrf(). # ------------------------------------------------------------------ -build_spec_mixed_mrf = function(x, data_columnnames, num_variables, - variable_type, is_ordinal, - baseline_category, - na_action, sampler, - interaction_prior_type, pairwise_scale, - interaction_alpha, interaction_beta, - threshold_prior_type, main_alpha, main_beta, - threshold_scale, - means_prior_type, means_scale, - means_alpha, means_beta, - scale_prior_type, scale_shape, scale_rate, - scale_eta = NA_real_, - delta = 0, - graph_prior_spec = "joint", - calibration_window = NULL, - edge_prior_flat) { +build_spec_mixed_mrf <- function(x, data_columnnames, num_variables, + variable_type, is_ordinal, + baseline_category, + na_action, sampler, + interaction_prior_type, pairwise_scale, + interaction_alpha, interaction_beta, + threshold_prior_type, main_alpha, main_beta, + threshold_scale, + means_prior_type, means_scale, + means_alpha, means_beta, + scale_prior_type, scale_shape, scale_rate, + scale_eta = NA_real_, + delta = 0, + precision_graph_prior = "joint", + calibration_window = NULL, + edge_prior_flat) { # Standardized-frame scale prior: derive the raw diagonal rate eta / s - scale_rate = resolve_scale_rate(scale_rate, scale_eta, pairwise_scale) + scale_rate <- resolve_scale_rate(scale_rate, scale_eta, pairwise_scale) # Identify discrete vs continuous columns - cont_idx = which(variable_type == "continuous") - disc_idx = which(variable_type != "continuous") - p = length(disc_idx) - q = length(cont_idx) + cont_idx <- which(variable_type == "continuous") + disc_idx <- which(variable_type != "continuous") + p <- length(disc_idx) + q <- length(cont_idx) # Split data - x_disc = x[, disc_idx, drop = FALSE] - x_cont = x[, cont_idx, drop = FALSE] + x_disc <- x[, disc_idx, drop = FALSE] + x_cont <- x[, cont_idx, drop = FALSE] # Ensure integer matrix for discrete data - storage.mode(x_disc) = "integer" + storage.mode(x_disc) <- "integer" # Ensure numeric matrix for continuous data - storage.mode(x_cont) = "double" + storage.mode(x_cont) <- "double" # Discrete variable properties (subset to discrete columns) - is_ordinal_disc = is_ordinal[disc_idx] - vtype_disc = variable_type[disc_idx] + is_ordinal_disc <- is_ordinal[disc_idx] + vtype_disc <- variable_type[disc_idx] # Subset baseline_category to discrete columns when the user supplies a # full-length vector (one entry per variable, including continuous ones). - if(length(baseline_category) == num_variables && num_variables != p) { - baseline_category = baseline_category[disc_idx] + if (length(baseline_category) == num_variables && num_variables != p) { + baseline_category <- baseline_category[disc_idx] } # Baseline category for discrete variables - bc = validate_baseline_category( + bc <- validate_baseline_category( baseline_category = baseline_category, baseline_category_provided = !is.null(baseline_category), x = x_disc, @@ -269,30 +269,30 @@ build_spec_mixed_mrf = function(x, data_columnnames, num_variables, ) # Missing data handling - na_impute = FALSE - missing_index_discrete = NULL - missing_index_continuous = NULL + na_impute <- FALSE + missing_index_discrete <- NULL + missing_index_continuous <- NULL - if(na_action == "listwise") { - missing_rows = apply(x_disc, 1, anyNA) | apply(x_cont, 1, anyNA) - if(all(missing_rows)) { + if (na_action == "listwise") { + missing_rows <- apply(x_disc, 1, anyNA) | apply(x_cont, 1, anyNA) + if (all(missing_rows)) { stop(paste0( "All rows in x contain at least one missing response.\n", "You could try option na_action = \"impute\"." )) } - n_removed = sum(missing_rows) - if(n_removed > 0 && isTRUE(getOption("bgms.verbose", TRUE))) { - n_remaining = nrow(x_disc) - n_removed + n_removed <- sum(missing_rows) + if (n_removed > 0 && isTRUE(getOption("bgms.verbose", TRUE))) { + n_remaining <- nrow(x_disc) - n_removed message( - n_removed, " row", if(n_removed > 1) "s" else "", + n_removed, " row", if (n_removed > 1) "s" else "", " with missing values excluded (n = ", n_remaining, " remaining).\n", "To impute missing values instead, use na_action = \"impute\"." ) } - x_disc = x_disc[!missing_rows, , drop = FALSE] - x_cont = x_cont[!missing_rows, , drop = FALSE] - if(nrow(x_disc) < 2) { + x_disc <- x_disc[!missing_rows, , drop = FALSE] + x_cont <- x_cont[!missing_rows, , drop = FALSE] + if (nrow(x_disc) < 2) { stop(paste0( "After removing missing observations from the input matrix x,\n", "there were less than two rows left in x." @@ -300,27 +300,27 @@ build_spec_mixed_mrf = function(x, data_columnnames, num_variables, } } else { # Impute path: handle discrete and continuous sub-matrices separately - md_disc = handle_impute(x_disc) - md_cont = handle_impute(x_cont) - x_disc = md_disc$x - x_cont = md_cont$x - na_impute = md_disc$na_impute || md_cont$na_impute - if(md_disc$na_impute) missing_index_discrete = md_disc$missing_index - if(md_cont$na_impute) missing_index_continuous = md_cont$missing_index + md_disc <- handle_impute(x_disc) + md_cont <- handle_impute(x_cont) + x_disc <- md_disc$x + x_cont <- md_cont$x + na_impute <- md_disc$na_impute || md_cont$na_impute + if (md_disc$na_impute) missing_index_discrete <- md_disc$missing_index + if (md_cont$na_impute) missing_index_continuous <- md_cont$missing_index } # Ordinal recoding (reformat discrete data) - ord = reformat_ordinal_data( + ord <- reformat_ordinal_data( x = x_disc, is_ordinal = is_ordinal_disc, baseline_category = bc ) - x_disc_recoded = ord$x - num_categories = ord$num_categories - bc_final = ord$baseline_category + x_disc_recoded <- ord$x + num_categories <- ord$num_categories + bc_final <- ord$baseline_category - ep = edge_prior_flat + ep <- edge_prior_flat - num_thresholds = sum(ifelse(is_ordinal_disc, num_categories, 2L)) + num_thresholds <- sum(ifelse(is_ordinal_disc, num_categories, 2L)) new_bgm_spec( model_type = "mixed_mrf", @@ -375,7 +375,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/run_sampler.R b/R/run_sampler.R index bba332ef..673b68b6 100644 --- a/R/run_sampler.R +++ b/R/run_sampler.R @@ -61,29 +61,29 @@ run_sampler = function(spec) { # ============================================================================== # run_sampler_ggm() # ============================================================================== -run_sampler_ggm = function(spec) { - d = spec$data - p = spec$prior - s = spec$sampler - m = spec$missing +run_sampler_ggm <- function(spec) { + d <- spec$data + p <- spec$prior + s <- spec$sampler + m <- spec$missing - bb_alpha_between = bb_between_or_sentinel(p$beta_bernoulli_alpha_between) - bb_beta_between = bb_between_or_sentinel(p$beta_bernoulli_beta_between) + bb_alpha_between <- bb_between_or_sentinel(p$beta_bernoulli_alpha_between) + bb_beta_between <- bb_between_or_sentinel(p$beta_bernoulli_beta_between) # Graph-prior specification: the joint path corrects the hyperparameter # updates for the untracked normalizer Z(Gamma); the hierarchical path # tracks Z(Gamma) itself in the between-edge moves (per-edge Z-ratio # engine) and keeps the hyperparameter updates clean conjugate. The two # are mutually exclusive. - correction = NULL - zratio = NULL - if(identical(p$graph_prior_spec, "hierarchical")) { - zc = zratio_constants( + correction <- NULL + zratio <- NULL + if (identical(p$precision_graph_prior, "hierarchical")) { + zc <- zratio_constants( delta = p$delta, sigma = 2 * p$pairwise_scale, beta = p$scale_rate / 2 ) - zratio = list( + zratio <- list( addc = zc$addc, tg = zc$tg, ihat = zc$ihat, ghat = zc$ghat, wt = zc$wt, psi0 = zc$psi0, delta = zc$delta, sigma = zc$sigma, beta = zc$beta, @@ -92,10 +92,10 @@ run_sampler_ggm = function(spec) { ) ) } else { - correction = ggm_edge_prior_correction(p, s, d$num_variables) + correction <- ggm_edge_prior_correction(p, s, d$num_variables) } - out_raw = sample_ggm( + out_raw <- sample_ggm( inputFromR = list( X = d$x, pairwise_scale = p$pairwise_scale, @@ -206,28 +206,28 @@ run_sampler_omrf = function(spec) { # ============================================================================== # run_sampler_mixed_mrf() # ============================================================================== -run_sampler_mixed_mrf = function(spec) { - d = spec$data - v = spec$variables - m = spec$missing - p = spec$prior - s = spec$sampler +run_sampler_mixed_mrf <- function(spec) { + d <- spec$data + v <- spec$variables + m <- spec$missing + p <- spec$prior + s <- spec$sampler - bb_alpha_between = bb_between_or_sentinel(p$beta_bernoulli_alpha_between) - bb_beta_between = bb_between_or_sentinel(p$beta_bernoulli_beta_between) + bb_alpha_between <- bb_between_or_sentinel(p$beta_bernoulli_alpha_between) + bb_beta_between <- bb_between_or_sentinel(p$beta_bernoulli_beta_between) # Graph-prior specification on the continuous block: same dichotomy as the # GGM path (see run_sampler_ggm). The Z-ratio constants and the Stage-3d # window are sized on the continuous subgraph. - correction = NULL - zratio = NULL - if(identical(p$graph_prior_spec, "hierarchical")) { - zc = zratio_constants( + correction <- NULL + zratio <- NULL + if (identical(p$precision_graph_prior, "hierarchical")) { + zc <- zratio_constants( delta = p$delta, sigma = 2 * p$pairwise_scale, beta = p$scale_rate / 2 ) - zratio = list( + zratio <- list( addc = zc$addc, tg = zc$tg, ihat = zc$ihat, ghat = zc$ghat, wt = zc$wt, psi0 = zc$psi0, delta = zc$delta, sigma = zc$sigma, beta = zc$beta, @@ -236,12 +236,12 @@ run_sampler_mixed_mrf = function(spec) { ) ) } else { - correction = ggm_edge_prior_correction( + correction <- ggm_edge_prior_correction( p, s, d$num_variables, d$num_continuous ) } - input_list = list( + input_list <- list( discrete_observations = d$x_discrete, continuous_observations = d$x_continuous, num_categories = d$num_categories, @@ -264,7 +264,7 @@ run_sampler_mixed_mrf = function(spec) { scale_rate = p$scale_rate ) - out_raw = sample_mixed_mrf( + out_raw <- sample_mixed_mrf( inputFromR = input_list, prior_inclusion_prob = p$inclusion_probability, initial_edge_indicators = matrix(1L, diff --git a/man/bgm.Rd b/man/bgm.Rd index 557d28ba..c120096a 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 @@ -215,7 +215,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 40bbe00d..b4363fb6 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) { @@ -7,13 +7,13 @@ hier_test_data = function(q = 10, n = 40, seed = 4) { } test_that("hierarchical spec eligibility is validated", { - Y = hier_test_data() + Y <- hier_test_data() expect_error( bgm( x = Y, variable_type = "continuous", interaction_prior = beta_prime_prior(), - graph_prior_spec = "hierarchical", + precision_graph_prior = "hierarchical", update_method = "gibbs", display_progress = "none", verbose = FALSE ), "normal or Cauchy" @@ -23,7 +23,7 @@ test_that("hierarchical spec eligibility is validated", { 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,7 +43,7 @@ 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" @@ -52,48 +52,48 @@ test_that("hierarchical spec eligibility is validated", { 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( + 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), - graph_prior_spec = "hierarchical", calibration_window = 50, + precision_graph_prior = "hierarchical", calibration_window = 50, update_method = method, chains = 1, cores = 1, seed = 7, display_progress = "none", verbose = FALSE ) - s = summary(fit) + 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$graph_prior_spec, "hierarchical", 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) - fit = bgm( + Y <- hier_test_data(q = 12) + fit <- bgm( x = Y, variable_type = "continuous", iter = 150, warmup = 250, 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 ) - zd = fit@zratio_diag + zd <- fit@zratio_diag expect_false(is.null(zd)) expect_equal(nrow(zd$per_chain), 2L) expect_true(all(zd$per_chain$frozen)) 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) @@ -101,14 +101,14 @@ test_that("bgm fits the hierarchical spec and attaches the alarm suite", { test_that("the joint default is unchanged", { skip_on_cran() - Y = hier_test_data(q = 6) - fit = bgm( + Y <- hier_test_data(q = 6) + fit <- bgm( x = Y, variable_type = "continuous", iter = 100, warmup = 150, 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) }) @@ -135,38 +135,38 @@ test_that("mixed indicator layout rebuilds the continuous subgraph", { test_that("mixed data supports the hierarchical spec on the continuous block", { skip_on_cran() set.seed(9) - n = 60 - X = cbind( + n <- 60 + X <- cbind( matrix(sample(0:2, n * 3, replace = TRUE), n, 3), matrix(rnorm(n * 8), n, 8) ) - colnames(X) = paste0("V", seq_len(11)) - vt = c(rep("ordinal", 3), rep("continuous", 8)) + colnames(X) <- paste0("V", seq_len(11)) + vt <- c(rep("ordinal", 3), rep("continuous", 8)) expect_error( bgm( 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" ) - fit = bgm( + fit <- bgm( x = X, variable_type = vt, iter = 120, warmup = 200, 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 ) - zd = fit@zratio_diag + zd <- fit@zratio_diag 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..cddb82fa 100644 --- a/tests/testthat/test-build-arguments.R +++ b/tests/testthat/test-build-arguments.R @@ -89,12 +89,12 @@ spec_compare = function(...) { # ============================================================================== test_that("GGM build_arguments: all expected field names present", { - s = spec_ggm() - a = build_arguments(s) - expected = c( + s <- spec_ggm() + a <- build_arguments(s) + 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", From af4404baf0d26bb8b5e0838c9682de2cebc34a27 Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Sun, 5 Jul 2026 21:10:29 +0200 Subject: [PATCH 3/5] Format the renamed files with the package code style The rename commit reformatted these files with the default styler settings (<- assignment, space after if). Reapply the package styler configuration so they use = assignment and no space after if, matching the rest of the codebase. --- R/bgm.R | 84 ++++---- R/bgm_spec.R | 246 +++++++++++------------ R/build_arguments.R | 5 +- R/build_output_bgm.R | 272 +++++++++++++------------- R/build_output_mixed_mrf.R | 256 ++++++++++++------------ R/build_spec.R | 150 +++++++------- R/run_sampler.R | 60 +++--- tests/testthat/test-bgm-hier-spec.R | 32 +-- tests/testthat/test-build-arguments.R | 6 +- 9 files changed, 554 insertions(+), 557 deletions(-) diff --git a/R/bgm.R b/R/bgm.R index 3baf6ec6..9b04e73e 100644 --- a/R/bgm.R +++ b/R/bgm.R @@ -395,7 +395,7 @@ #' } #' #' @export -bgm <- function( +bgm = function( x, variable_type = "ordinal", baseline_category, @@ -441,87 +441,87 @@ bgm <- function( ) { # Set verbose option for internal functions, restore on exit - old_verbose <- getOption("bgms.verbose") + old_verbose = getOption("bgms.verbose") options(bgms.verbose = verbose) on.exit(options(bgms.verbose = old_verbose), add = TRUE) # --- Legacy deprecation: v0.1.6.0 renames ----------------------------------- - if (hasArg(interaction_scale)) { + if(hasArg(interaction_scale)) { lifecycle::deprecate_warn( "0.1.6.0", "bgm(interaction_scale =)", "bgm(interaction_prior =)" ) - if (!hasArg(pairwise_scale) && + if(!hasArg(pairwise_scale) && identical(interaction_prior, cauchy_prior(scale = 1))) { - interaction_prior <- cauchy_prior(scale = interaction_scale) + interaction_prior = cauchy_prior(scale = interaction_scale) } } - if (hasArg(burnin)) { + if(hasArg(burnin)) { lifecycle::deprecate_warn("0.1.6.0", "bgm(burnin =)", "bgm(warmup =)") - if (!hasArg(warmup)) warmup <- burnin + if(!hasArg(warmup)) warmup = burnin } - if (hasArg(save)) { + if(hasArg(save)) { lifecycle::deprecate_warn("0.1.6.0", "bgm(save =)") } - if (hasArg(threshold_alpha) || hasArg(threshold_beta)) { + if(hasArg(threshold_alpha) || hasArg(threshold_beta)) { lifecycle::deprecate_warn( "0.1.6.0", "bgm(threshold_alpha =, threshold_beta =)", "bgm(threshold_prior =)" ) - if (identical(threshold_prior, beta_prime_prior(0.5, 0.5))) { - ta <- if (hasArg(threshold_alpha)) threshold_alpha else 0.5 - tb <- if (hasArg(threshold_beta)) threshold_beta else 0.5 - threshold_prior <- beta_prime_prior(alpha = ta, beta = tb) + if(identical(threshold_prior, beta_prime_prior(0.5, 0.5))) { + ta = if(hasArg(threshold_alpha)) threshold_alpha else 0.5 + tb = if(hasArg(threshold_beta)) threshold_beta else 0.5 + threshold_prior = beta_prime_prior(alpha = ta, beta = tb) } } # --- Legacy deprecation: scalar prior parameters ---------------------------- - if (hasArg(pairwise_scale)) { + if(hasArg(pairwise_scale)) { lifecycle::deprecate_warn( "0.2.0", "bgm(pairwise_scale =)", "bgm(interaction_prior =)" ) - if (identical(interaction_prior, cauchy_prior(scale = 1))) { - interaction_prior <- cauchy_prior(scale = pairwise_scale) + if(identical(interaction_prior, cauchy_prior(scale = 1))) { + interaction_prior = cauchy_prior(scale = pairwise_scale) } } - if (hasArg(main_alpha) || hasArg(main_beta)) { + if(hasArg(main_alpha) || hasArg(main_beta)) { lifecycle::deprecate_warn( "0.2.0", "bgm(main_alpha =)", "bgm(threshold_prior =)" ) - if (identical(threshold_prior, beta_prime_prior(0.5, 0.5))) { - ma <- if (hasArg(main_alpha)) main_alpha else 0.5 - mb <- if (hasArg(main_beta)) main_beta else 0.5 - threshold_prior <- beta_prime_prior(alpha = ma, beta = mb) + if(identical(threshold_prior, beta_prime_prior(0.5, 0.5))) { + ma = if(hasArg(main_alpha)) main_alpha else 0.5 + mb = if(hasArg(main_beta)) main_beta else 0.5 + threshold_prior = beta_prime_prior(alpha = ma, beta = mb) } } # Handle edge_prior: accept both string (deprecated) and object (new) - if (is.character(edge_prior)) { + if(is.character(edge_prior)) { lifecycle::deprecate_warn( "0.2.0", "bgm(edge_prior = 'must be a prior object')", "bgm(edge_prior = 'bernoulli_prior()')" ) - edge_prior_str <- match.arg(edge_prior, + edge_prior_str = match.arg(edge_prior, choices = c("Bernoulli", "Beta-Bernoulli", "Stochastic-Block") ) - ip <- if (hasArg(inclusion_probability)) inclusion_probability else 0.5 - bba <- if (hasArg(beta_bernoulli_alpha)) beta_bernoulli_alpha else 1 - bbb <- if (hasArg(beta_bernoulli_beta)) beta_bernoulli_beta else 1 - bbab <- if (hasArg(beta_bernoulli_alpha_between)) beta_bernoulli_alpha_between else 1 - bbbb <- if (hasArg(beta_bernoulli_beta_between)) beta_bernoulli_beta_between else 1 - da <- if (hasArg(dirichlet_alpha)) dirichlet_alpha else 1 - lam <- if (hasArg(lambda)) lambda else 1 + ip = if(hasArg(inclusion_probability)) inclusion_probability else 0.5 + bba = if(hasArg(beta_bernoulli_alpha)) beta_bernoulli_alpha else 1 + bbb = if(hasArg(beta_bernoulli_beta)) beta_bernoulli_beta else 1 + bbab = if(hasArg(beta_bernoulli_alpha_between)) beta_bernoulli_alpha_between else 1 + bbbb = if(hasArg(beta_bernoulli_beta_between)) beta_bernoulli_beta_between else 1 + da = if(hasArg(dirichlet_alpha)) dirichlet_alpha else 1 + lam = if(hasArg(lambda)) lambda else 1 - edge_prior <- switch(edge_prior_str, + edge_prior = switch(edge_prior_str, "Bernoulli" = bernoulli_prior(inclusion_probability = ip), "Beta-Bernoulli" = beta_bernoulli_prior(alpha = bba, beta = bbb), "Stochastic-Block" = sbm_prior( @@ -532,13 +532,13 @@ bgm <- function( ) } else { # Warn if loose edge params are also supplied alongside an object - if (hasArg(inclusion_probability)) { + if(hasArg(inclusion_probability)) { lifecycle::deprecate_warn( "0.2.0", "bgm(inclusion_probability =)", "bgm(edge_prior = 'bernoulli_prior()')" ) } - if (hasArg(beta_bernoulli_alpha)) { + if(hasArg(beta_bernoulli_alpha)) { lifecycle::deprecate_warn( "0.2.0", "bgm(beta_bernoulli_alpha =)", "bgm(edge_prior = 'beta_bernoulli_prior()')" @@ -547,17 +547,17 @@ bgm <- function( } # --- Unpack prior objects to flat parameters for bgm_spec ------------------- - ip <- unpack_interaction_prior(interaction_prior) - tp <- unpack_threshold_prior(threshold_prior) - mp <- unpack_parameter_prior(means_prior) - sp <- unpack_scale_prior(precision_scale_prior) + ip = unpack_interaction_prior(interaction_prior) + tp = unpack_threshold_prior(threshold_prior) + mp = unpack_parameter_prior(means_prior) + sp = unpack_scale_prior(precision_scale_prior) # --- Build spec, sample, build output ---------------------------------------- - spec <- bgm_spec( + spec = bgm_spec( x = x, model_type = "omrf", variable_type = variable_type, - baseline_category = if (hasArg(baseline_category)) baseline_category else NULL, + baseline_category = if(hasArg(baseline_category)) baseline_category else NULL, na_action = na_action, interaction_prior_type = ip$interaction_prior_type, pairwise_scale = ip$pairwise_scale, @@ -581,7 +581,7 @@ bgm <- function( precision_graph_prior = precision_graph_prior, calibration_window = calibration_window, update_method = update_method, - target_accept = if (hasArg(target_accept)) target_accept else NULL, + target_accept = if(hasArg(target_accept)) target_accept else NULL, iter = iter, warmup = warmup, nuts_max_depth = nuts_max_depth, @@ -594,7 +594,7 @@ bgm <- function( progress_callback = progress_callback ) - raw <- run_sampler(spec) - output <- build_output(spec, raw) + raw = run_sampler(spec) + output = build_output(spec, raw) return(output) } diff --git a/R/bgm_spec.R b/R/bgm_spec.R index f4b4f6f1..0c73f475 100644 --- a/R/bgm_spec.R +++ b/R/bgm_spec.R @@ -238,79 +238,79 @@ validate_bgm_spec = function(spec) { # # Parameters mirror the union of bgm() and bgmCompare() arguments. # ============================================================================== -bgm_spec <- function(x, - model_type = c("omrf", "ggm", "compare", "mixed_mrf"), - # Variable specification - variable_type = "ordinal", - baseline_category = NULL, - # Data (compare-specific) - y = NULL, - group_indicator = NULL, - # Missing data - na_action = c("listwise", "impute"), - # Priors (new: prior objects unpacked by bgm()) - interaction_prior_type = "cauchy", - pairwise_scale = 1, - interaction_alpha = NA_real_, - interaction_beta = NA_real_, - threshold_prior_type = "beta-prime", - main_alpha = 0.5, - main_beta = 0.5, - threshold_scale = NA_real_, - means_prior_type = "normal", - means_scale = 1, - means_alpha = NA_real_, - means_beta = NA_real_, - scale_prior_type = "gamma", - scale_shape = 1, - scale_rate = 1, - scale_eta = NA_real_, - delta = NULL, - edge_selection = TRUE, - edge_prior = bernoulli_prior(0.5), - precision_graph_prior = c("joint", "hierarchical"), - calibration_window = NULL, - # Legacy edge prior params (accepted for backward compat) - inclusion_probability = 0.5, - beta_bernoulli_alpha_between = 1, - beta_bernoulli_beta_between = 1, - dirichlet_alpha = 1, - lambda = 1, - # Priors (compare-specific) - difference_selection = TRUE, - main_difference_selection = FALSE, - difference_prior = c( - "Bernoulli", "Beta-Bernoulli", "Stochastic-Block" - ), - difference_scale = 1, - difference_prior_type = "cauchy", - difference_probability = 0.5, - # Compare difference prior hyperparameters - beta_bernoulli_alpha = 1, - beta_bernoulli_beta = 1, - difference_beta_bernoulli_alpha_between = 1, - difference_beta_bernoulli_beta_between = 1, - difference_dirichlet_alpha = 1, - difference_lambda = 1, - # Sampler - update_method = c( - "nuts", - "adaptive-metropolis", - "gibbs" - ), - target_accept = NULL, - iter = 10000L, - warmup = 1000L, - nuts_max_depth = 10L, - learn_mass_matrix = TRUE, - chains = 4L, - cores = parallel::detectCores(), - seed = NULL, - display_progress = c("per-chain", "total", "none"), - verbose = TRUE, - progress_callback = NULL) { - model_type <- match.arg(model_type) - na_action <- tryCatch(match.arg(na_action), error = function(e) { +bgm_spec = function(x, + model_type = c("omrf", "ggm", "compare", "mixed_mrf"), + # Variable specification + variable_type = "ordinal", + baseline_category = NULL, + # Data (compare-specific) + y = NULL, + group_indicator = NULL, + # Missing data + na_action = c("listwise", "impute"), + # Priors (new: prior objects unpacked by bgm()) + interaction_prior_type = "cauchy", + pairwise_scale = 1, + interaction_alpha = NA_real_, + interaction_beta = NA_real_, + threshold_prior_type = "beta-prime", + main_alpha = 0.5, + main_beta = 0.5, + threshold_scale = NA_real_, + means_prior_type = "normal", + means_scale = 1, + means_alpha = NA_real_, + means_beta = NA_real_, + scale_prior_type = "gamma", + scale_shape = 1, + scale_rate = 1, + scale_eta = NA_real_, + delta = NULL, + edge_selection = TRUE, + edge_prior = bernoulli_prior(0.5), + precision_graph_prior = c("joint", "hierarchical"), + calibration_window = NULL, + # Legacy edge prior params (accepted for backward compat) + inclusion_probability = 0.5, + beta_bernoulli_alpha_between = 1, + beta_bernoulli_beta_between = 1, + dirichlet_alpha = 1, + lambda = 1, + # Priors (compare-specific) + difference_selection = TRUE, + main_difference_selection = FALSE, + difference_prior = c( + "Bernoulli", "Beta-Bernoulli", "Stochastic-Block" + ), + difference_scale = 1, + difference_prior_type = "cauchy", + difference_probability = 0.5, + # Compare difference prior hyperparameters + beta_bernoulli_alpha = 1, + beta_bernoulli_beta = 1, + difference_beta_bernoulli_alpha_between = 1, + difference_beta_bernoulli_beta_between = 1, + difference_dirichlet_alpha = 1, + difference_lambda = 1, + # Sampler + update_method = c( + "nuts", + "adaptive-metropolis", + "gibbs" + ), + target_accept = NULL, + iter = 10000L, + warmup = 1000L, + nuts_max_depth = 10L, + learn_mass_matrix = TRUE, + chains = 4L, + cores = parallel::detectCores(), + seed = NULL, + display_progress = c("per-chain", "total", "none"), + verbose = TRUE, + progress_callback = NULL) { + model_type = match.arg(model_type) + na_action = tryCatch(match.arg(na_action), error = function(e) { stop(paste0( "The na_action argument should be one of \"listwise\" or \"impute\", not \"", na_action, "\"." @@ -318,44 +318,44 @@ bgm_spec <- function(x, }) # --- Data validation -------------------------------------------------------- - x <- data_check(x, "x") - data_columnnames <- if (is.null(colnames(x))) { + x = data_check(x, "x") + data_columnnames = if(is.null(colnames(x))) { paste0("Variable ", seq_len(ncol(x))) } else { colnames(x) } - num_variables <- ncol(x) + num_variables = ncol(x) # --- Variable types --------------------------------------------------------- - allow_continuous <- (model_type != "compare") - vt <- validate_variable_types( + allow_continuous = (model_type != "compare") + vt = validate_variable_types( variable_type = variable_type, num_variables = num_variables, allow_continuous = allow_continuous, allow_mixed = (model_type != "compare"), - caller = if (model_type == "compare") "bgmCompare" else "bgm" + caller = if(model_type == "compare") "bgmCompare" else "bgm" ) - variable_type <- vt$variable_type - is_ordinal <- vt$variable_bool - is_continuous <- vt$is_continuous - is_mixed <- vt$is_mixed + variable_type = vt$variable_type + is_ordinal = vt$variable_bool + is_continuous = vt$is_continuous + is_mixed = vt$is_mixed # Resolve model_type if "omrf" default was kept but data is continuous - if (model_type == "omrf" && is_continuous) { - model_type <- "ggm" + if(model_type == "omrf" && is_continuous) { + model_type = "ggm" } - if (model_type == "omrf" && is_mixed) { - model_type <- "mixed_mrf" + if(model_type == "omrf" && is_mixed) { + model_type = "mixed_mrf" } # Auto-resolve delta = NULL to the dimension-adaptive default 0.5 * log(p), # where p is the dimension of the continuous precision matrix. For models # without a continuous block (omrf, compare) the tilt has no target, so # NULL resolves to 0. - if (is.null(delta)) { - delta <- if (model_type == "ggm") { + if(is.null(delta)) { + delta = if(model_type == "ggm") { 0.5 * log(max(num_variables, 1)) - } else if (model_type == "mixed_mrf") { + } else if(model_type == "mixed_mrf") { 0.5 * log(max(sum(variable_type == "continuous"), 1)) } else { 0 @@ -363,11 +363,11 @@ bgm_spec <- function(x, } # Validate determinant-tilt exponent and reject for pure-ordinal models - if (!is.numeric(delta) || length(delta) != 1L || is.na(delta) || + if(!is.numeric(delta) || length(delta) != 1L || is.na(delta) || !is.finite(delta) || delta < 0) { stop("'delta' must be a single finite non-negative numeric, or NULL.") } - if (delta > 0 && model_type %in% c("omrf", "compare")) { + if(delta > 0 && model_type %in% c("omrf", "compare")) { stop( "'delta' (determinant tilt) requires continuous variables; the ", "current model_type is '", model_type, "', which has no precision ", @@ -379,16 +379,16 @@ 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. - precision_graph_prior <- match.arg(precision_graph_prior) - if (precision_graph_prior == "hierarchical") { - if (!model_type %in% c("ggm", "mixed_mrf")) { + precision_graph_prior = match.arg(precision_graph_prior) + if(precision_graph_prior == "hierarchical") { + if(!model_type %in% c("ggm", "mixed_mrf")) { stop( "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) { + if(model_type == "mixed_mrf" && sum(variable_type == "continuous") < 2) { stop( "precision_graph_prior = \"hierarchical\" on mixed data needs at least ", "two continuous variables (the specification normalizes the ", @@ -396,14 +396,14 @@ bgm_spec <- function(x, "specification." ) } - if (!edge_selection) { + if(!edge_selection) { stop( "precision_graph_prior = \"hierarchical\" normalizes p(K | Gamma) ", "across graphs and needs edge_selection = TRUE; with a fixed ", "graph the specifications coincide." ) } - if (!interaction_prior_type %in% c("normal", "cauchy")) { + if(!interaction_prior_type %in% c("normal", "cauchy")) { stop(sprintf( paste0( "precision_graph_prior = \"hierarchical\" supports a normal or Cauchy ", @@ -416,7 +416,7 @@ bgm_spec <- function(x, interaction_prior_type )) } - if (abs(scale_shape - 1) > 1e-12) { + if(abs(scale_shape - 1) > 1e-12) { stop(sprintf( paste0( "precision_graph_prior = \"hierarchical\" requires shape = 1 on the ", @@ -431,7 +431,7 @@ bgm_spec <- function(x, } # --- Sampler (needs is_continuous and edge_selection early) ------------------ - sampler <- validate_sampler( + sampler = validate_sampler( update_method = update_method, target_accept = target_accept, iter = iter, @@ -443,20 +443,20 @@ bgm_spec <- function(x, seed = seed, display_progress = display_progress, is_continuous = is_continuous, - edge_selection = if (model_type == "compare") FALSE else edge_selection, + edge_selection = if(model_type == "compare") FALSE else edge_selection, verbose = verbose, progress_callback = progress_callback ) # --- Resolve edge prior object ----------------------------------------------- - if (inherits(edge_prior, "bgms_indicator_prior")) { - ep_flat <- unpack_indicator_prior(edge_prior, num_variables) - } else if (is.character(edge_prior)) { + if(inherits(edge_prior, "bgms_indicator_prior")) { + ep_flat = unpack_indicator_prior(edge_prior, num_variables) + } else if(is.character(edge_prior)) { # Legacy string path (tests and bgmCompare may call bgm_spec directly) - edge_prior_str <- match.arg(edge_prior, + edge_prior_str = match.arg(edge_prior, choices = c("Bernoulli", "Beta-Bernoulli", "Stochastic-Block") ) - ep_flat <- validate_edge_prior( + ep_flat = validate_edge_prior( edge_selection = edge_selection, edge_prior = edge_prior_str, inclusion_probability = inclusion_probability, num_variables = num_variables, @@ -466,12 +466,12 @@ bgm_spec <- function(x, beta_bernoulli_beta_between = beta_bernoulli_beta_between, dirichlet_alpha = dirichlet_alpha, lambda = lambda ) - ep_flat$beta_bernoulli_alpha <- beta_bernoulli_alpha - ep_flat$beta_bernoulli_beta <- beta_bernoulli_beta - ep_flat$beta_bernoulli_alpha_between <- beta_bernoulli_alpha_between - ep_flat$beta_bernoulli_beta_between <- beta_bernoulli_beta_between - ep_flat$dirichlet_alpha <- dirichlet_alpha - ep_flat$lambda <- lambda + ep_flat$beta_bernoulli_alpha = beta_bernoulli_alpha + ep_flat$beta_bernoulli_beta = beta_bernoulli_beta + ep_flat$beta_bernoulli_alpha_between = beta_bernoulli_alpha_between + ep_flat$beta_bernoulli_beta_between = beta_bernoulli_beta_between + ep_flat$dirichlet_alpha = dirichlet_alpha + ep_flat$lambda = lambda } else { stop( "'edge_prior' must be a bgms_indicator_prior object.", @@ -479,15 +479,15 @@ bgm_spec <- function(x, ) } # Override edge_selection if explicitly FALSE - if (!edge_selection) { - ep_flat$edge_selection <- FALSE - ep_flat$edge_prior <- "Not Applicable" - ep_flat$inclusion_probability <- matrix(0.5, nrow = 1, ncol = 1) + if(!edge_selection) { + ep_flat$edge_selection = FALSE + ep_flat$edge_prior = "Not Applicable" + ep_flat$inclusion_probability = matrix(0.5, nrow = 1, ncol = 1) } # --- Build by model type ---------------------------------------------------- - if (model_type == "ggm") { - spec <- build_spec_ggm( + if(model_type == "ggm") { + spec = build_spec_ggm( x = x, data_columnnames = data_columnnames, num_variables = num_variables, variable_type = variable_type, is_ordinal = is_ordinal, @@ -507,8 +507,8 @@ bgm_spec <- function(x, calibration_window = calibration_window, edge_prior_flat = ep_flat ) - } else if (model_type == "mixed_mrf") { - spec <- build_spec_mixed_mrf( + } else if(model_type == "mixed_mrf") { + spec = build_spec_mixed_mrf( x = x, data_columnnames = data_columnnames, num_variables = num_variables, variable_type = variable_type, is_ordinal = is_ordinal, @@ -534,8 +534,8 @@ bgm_spec <- function(x, calibration_window = calibration_window, edge_prior_flat = ep_flat ) - } else if (model_type == "omrf") { - spec <- build_spec_omrf( + } else if(model_type == "omrf") { + spec = build_spec_omrf( x = x, data_columnnames = data_columnnames, num_variables = num_variables, variable_type = variable_type, is_ordinal = is_ordinal, @@ -552,7 +552,7 @@ bgm_spec <- function(x, edge_prior_flat = ep_flat ) } else { - spec <- build_spec_compare( + spec = build_spec_compare( x = x, y = y, group_indicator = group_indicator, data_columnnames = data_columnnames, num_variables = num_variables, diff --git a/R/build_arguments.R b/R/build_arguments.R index 94b7c7f7..2568f908 100644 --- a/R/build_arguments.R +++ b/R/build_arguments.R @@ -27,8 +27,7 @@ build_arguments = function(spec) { } - -build_arguments_ggm <- function(spec) { +build_arguments_ggm = function(spec) { list( num_variables = spec$data$num_variables, num_cases = spec$data$num_cases, @@ -103,7 +102,7 @@ build_arguments_omrf = function(spec) { } -build_arguments_mixed_mrf <- function(spec) { +build_arguments_mixed_mrf = function(spec) { list( num_variables = spec$data$num_variables, num_discrete = spec$data$num_discrete, diff --git a/R/build_output_bgm.R b/R/build_output_bgm.R index 94674d28..b632f4d7 100644 --- a/R/build_output_bgm.R +++ b/R/build_output_bgm.R @@ -11,22 +11,22 @@ # 1. Parameter naming: GGM uses "Var (precision)", OMRF uses "Var (k)" # 2. Main posterior mean shape: GGM = px1, OMRF = pxmax_categories # ============================================================================== -build_output_bgm <- function(spec, raw) { - d <- spec$data - v <- spec$variables - p <- spec$prior - s <- spec$sampler +build_output_bgm = function(spec, raw) { + d = spec$data + v = spec$variables + p = spec$prior + s = spec$sampler - is_continuous <- v$is_continuous - num_variables <- d$num_variables - data_columnnames <- d$data_columnnames - edge_selection <- p$edge_selection - edge_prior <- p$edge_prior + is_continuous = v$is_continuous + num_variables = d$num_variables + data_columnnames = d$data_columnnames + edge_selection = p$edge_selection + edge_prior = p$edge_prior # 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$precision_graph_prior, "hierarchical")) { + zratio_chains = if(identical(p$precision_graph_prior, "hierarchical")) { raw } else { NULL @@ -37,42 +37,42 @@ build_output_bgm <- function(spec, raw) { # via convert_results_to_list(). Split into main and pairwise components and # transpose to (iters x params) --- the layout that MCMC summary functions # expect. - if (is_continuous) { + if(is_continuous) { # GGM: samples contain the upper triangle of the precision matrix # (row-major). Diagonal entries are "main"; off-diagonal are "pairwise". - diag_idx <- integer(num_variables) - offdiag_idx <- integer(num_variables * (num_variables - 1L) / 2L) - pos <- 0L - di <- 0L - oi <- 0L - for (i in seq_len(num_variables)) { - for (j in i:num_variables) { - pos <- pos + 1L - if (i == j) { - di <- di + 1L - diag_idx[di] <- pos + diag_idx = integer(num_variables) + offdiag_idx = integer(num_variables * (num_variables - 1L) / 2L) + pos = 0L + di = 0L + oi = 0L + for(i in seq_len(num_variables)) { + for(j in i:num_variables) { + pos = pos + 1L + if(i == j) { + di = di + 1L + diag_idx[di] = pos } else { - oi <- oi + 1L - offdiag_idx[oi] <- pos + oi = oi + 1L + offdiag_idx[oi] = pos } } } - raw <- lapply(raw, function(chain) { - samples_t <- t(chain$samples) - res <- list( + raw = lapply(raw, function(chain) { + samples_t = t(chain$samples) + res = list( main_samples = samples_t[, diag_idx, drop = FALSE], pairwise_samples = samples_t[, offdiag_idx, drop = FALSE], userInterrupt = isTRUE(chain$userInterrupt), chain_id = chain$chain_id ) - if (!is.null(chain$indicator_samples)) { - res$indicator_samples <- t(chain$indicator_samples)[, offdiag_idx, drop = FALSE] + if(!is.null(chain$indicator_samples)) { + res$indicator_samples = t(chain$indicator_samples)[, offdiag_idx, drop = FALSE] } - if (!is.null(chain$allocation_samples)) { - res$allocations <- t(chain$allocation_samples) + if(!is.null(chain$allocation_samples)) { + res$allocations = t(chain$allocation_samples) } - if (!is.null(chain$inclusion_parameter_samples)) { - res$inclusion_parameter <- as.numeric(chain$inclusion_parameter_samples) + if(!is.null(chain$inclusion_parameter_samples)) { + res$inclusion_parameter = as.numeric(chain$inclusion_parameter_samples) } attach_diagnostic_traces(res, chain) }) @@ -80,52 +80,52 @@ build_output_bgm <- function(spec, raw) { # OMRF: the first num_thresholds params are main effects, the rest are # pairwise. NUTS diagnostics use bare names from C++; rename to the # trailing-__ convention expected by summarize_nuts_diagnostics(). - num_thresholds <- spec$precomputed$num_thresholds - raw <- lapply(raw, function(chain) { - samples_t <- t(chain$samples) - n_params <- ncol(samples_t) - res <- list( + num_thresholds = spec$precomputed$num_thresholds + raw = lapply(raw, function(chain) { + samples_t = t(chain$samples) + n_params = ncol(samples_t) + res = list( main_samples = samples_t[, seq_len(num_thresholds), drop = FALSE], pairwise_samples = samples_t[, seq(num_thresholds + 1L, n_params), drop = FALSE], userInterrupt = isTRUE(chain$userInterrupt), chain_id = chain$chain_id ) - if (!is.null(chain$indicator_samples)) { - res$indicator_samples <- t(chain$indicator_samples) + if(!is.null(chain$indicator_samples)) { + res$indicator_samples = t(chain$indicator_samples) } - if (!is.null(chain$allocation_samples)) { - res$allocations <- t(chain$allocation_samples) + if(!is.null(chain$allocation_samples)) { + res$allocations = t(chain$allocation_samples) } - if (!is.null(chain$inclusion_parameter_samples)) { - res$inclusion_parameter <- as.numeric(chain$inclusion_parameter_samples) + if(!is.null(chain$inclusion_parameter_samples)) { + res$inclusion_parameter = as.numeric(chain$inclusion_parameter_samples) } attach_diagnostic_traces(res, chain) }) } # --- Parameter names -------------------------------------------------------- - if (is_continuous) { + if(is_continuous) { # GGM: raw samples are on precision scale; label accordingly - names_main <- paste0(data_columnnames, " (precision)") - is_ordinal_variable <- NULL - num_categories <- NULL + names_main = paste0(data_columnnames, " (precision)") + is_ordinal_variable = NULL + num_categories = NULL } else { # OMRF: per-category thresholds or BC linear/quadratic - is_ordinal_variable <- v$is_ordinal - num_categories <- d$num_categories - names_main <- character() - for (vi in seq_len(num_variables)) { - if (is_ordinal_variable[vi]) { - cats <- ordinal_threshold_labels( + is_ordinal_variable = v$is_ordinal + num_categories = d$num_categories + names_main = character() + for(vi in seq_len(num_variables)) { + if(is_ordinal_variable[vi]) { + cats = ordinal_threshold_labels( num_categories[vi], - if (!is.null(d$category_levels)) d$category_levels[[vi]] else NULL + if(!is.null(d$category_levels)) d$category_levels[[vi]] else NULL ) - names_main <- c( + names_main = c( names_main, paste0(data_columnnames[vi], " (", cats, ")") ) } else { - names_main <- c( + names_main = c( names_main, paste0(data_columnnames[vi], " (linear)"), paste0(data_columnnames[vi], " (quadratic)") @@ -134,10 +134,10 @@ build_output_bgm <- function(spec, raw) { } } - edge_names <- character() - for (i in seq_len(num_variables - 1)) { - for (j in seq(i + 1, num_variables)) { - edge_names <- c( + edge_names = character() + for(i in seq_len(num_variables - 1)) { + for(j in seq(i + 1, num_variables)) { + edge_names = c( edge_names, paste0(data_columnnames[i], "-", data_columnnames[j]) ) @@ -148,129 +148,129 @@ build_output_bgm <- function(spec, raw) { # Store normalized raw chains and metadata so that ensure_summaries() can # compute ESS/Rhat/MCSE on first demand. Posterior_summary_* fields start # as NULL and are populated lazily. - cache <- new.env(parent = emptyenv()) - cache$raw <- raw - cache$edge_selection <- edge_selection - cache$names_main <- names_main - cache$edge_names <- edge_names - cache$is_continuous <- is_continuous - cache$model_type <- if (is_continuous) "ggm" else "omrf" - cache$summaries_computed <- FALSE + cache = new.env(parent = emptyenv()) + cache$raw = raw + cache$edge_selection = edge_selection + cache$names_main = names_main + cache$edge_names = edge_names + cache$is_continuous = is_continuous + cache$model_type = if(is_continuous) "ggm" else "omrf" + cache$summaries_computed = FALSE # --- Compute posterior means from raw samples (cheap) ----------------------- - pooled_main <- do.call(rbind, lapply(raw, function(ch) ch$main_samples)) - pooled_pair <- do.call(rbind, lapply(raw, function(ch) ch$pairwise_samples)) - main_means <- colMeans(pooled_main) - pair_means <- colMeans(pooled_pair) + pooled_main = do.call(rbind, lapply(raw, function(ch) ch$main_samples)) + pooled_pair = do.call(rbind, lapply(raw, function(ch) ch$pairwise_samples)) + main_means = colMeans(pooled_main) + pair_means = colMeans(pooled_pair) - results <- list() + results = list() # --- Edge selection summaries (deferred) ------------------------------------ - has_sbm <- FALSE - if (edge_selection) { - has_sbm <- identical(edge_prior, "Stochastic-Block") && + has_sbm = FALSE + if(edge_selection) { + has_sbm = identical(edge_prior, "Stochastic-Block") && "allocations" %in% names(raw[[1]]) - if (has_sbm) { - sbm_convergence <- summarize_alloc_pairs( + if(has_sbm) { + sbm_convergence = summarize_alloc_pairs( allocations = lapply(raw, `[[`, "allocations"), node_names = data_columnnames ) - results$posterior_summary_pairwise_allocations <- sbm_convergence$sbm_summary - co_occur_matrix <- sbm_convergence$co_occur_matrix + results$posterior_summary_pairwise_allocations = sbm_convergence$sbm_summary + co_occur_matrix = sbm_convergence$co_occur_matrix } - if ("inclusion_parameter" %in% names(raw[[1]])) { - results$inclusion_parameter_samples <- + if("inclusion_parameter" %in% names(raw[[1]])) { + results$inclusion_parameter_samples = lapply(raw, `[[`, "inclusion_parameter") } } # --- Posterior mean: main --------------------------------------------------- - if (is_continuous) { + if(is_continuous) { # GGM has no main effects - results$posterior_mean_main <- NULL + results$posterior_mean_main = NULL } else { # OMRF: p x max_categories matrix - num_params <- ifelse(is_ordinal_variable, num_categories, 2L) - max_num_categories <- max(num_params) + num_params = ifelse(is_ordinal_variable, num_categories, 2L) + max_num_categories = max(num_params) - pmm <- matrix(NA, nrow = num_variables, ncol = max_num_categories) - start <- 0L - stop <- 0L - for (vi in seq_len(num_variables)) { - if (is_ordinal_variable[vi]) { - start <- stop + 1L - stop <- start + num_categories[vi] - 1L - pmm[vi, seq_len(num_categories[vi])] <- main_means[start:stop] + pmm = matrix(NA, nrow = num_variables, ncol = max_num_categories) + start = 0L + stop = 0L + for(vi in seq_len(num_variables)) { + if(is_ordinal_variable[vi]) { + start = stop + 1L + stop = start + num_categories[vi] - 1L + pmm[vi, seq_len(num_categories[vi])] = main_means[start:stop] } else { - start <- stop + 1L - stop <- start + 1L - pmm[vi, 1:2] <- main_means[start:stop] + start = stop + 1L + stop = start + 1L + pmm[vi, 1:2] = main_means[start:stop] } } - results$posterior_mean_main <- pmm - rownames(results$posterior_mean_main) <- data_columnnames - colnames(results$posterior_mean_main) <- paste0("cat (", seq_len(ncol(pmm)), ")") + results$posterior_mean_main = pmm + rownames(results$posterior_mean_main) = data_columnnames + colnames(results$posterior_mean_main) = paste0("cat (", seq_len(ncol(pmm)), ")") } # --- Posterior mean: associations ------------------------------------------- # For GGM: C++ stores precision; convert to association scale (* -0.5). # For OMRF: C++ already stores association-scale values. - associations <- matrix(0, + associations = matrix(0, nrow = num_variables, ncol = num_variables, dimnames = list(data_columnnames, data_columnnames) ) - associations[lower.tri(associations)] <- pair_means - associations <- associations + t(associations) - if (is_continuous) { - associations <- -0.5 * associations + associations[lower.tri(associations)] = pair_means + associations = associations + t(associations) + if(is_continuous) { + associations = -0.5 * associations } - results$posterior_mean_pairwise <- associations + results$posterior_mean_pairwise = associations # --- Residual variance (GGM only) ------------------------------------------- # C++ stores precision diagonal; residual variance = 1 / precision. # Average per-sample inversions to avoid Jensen's inequality bias. - if (is_continuous) { - results$posterior_mean_residual_variance <- colMeans(1 / pooled_main) - names(results$posterior_mean_residual_variance) <- data_columnnames + if(is_continuous) { + results$posterior_mean_residual_variance = colMeans(1 / pooled_main) + names(results$posterior_mean_residual_variance) = data_columnnames } # --- Posterior mean: indicator + SBM ---------------------------------------- - if (edge_selection) { - pooled_ind <- do.call(rbind, lapply(raw, function(ch) ch$indicator_samples)) - indicator_means <- colMeans(pooled_ind) - results$posterior_mean_indicator <- matrix(0, + if(edge_selection) { + pooled_ind = do.call(rbind, lapply(raw, function(ch) ch$indicator_samples)) + indicator_means = colMeans(pooled_ind) + results$posterior_mean_indicator = matrix(0, nrow = num_variables, ncol = num_variables, dimnames = list(data_columnnames, data_columnnames) ) - results$posterior_mean_indicator[lower.tri(results$posterior_mean_indicator)] <- + results$posterior_mean_indicator[lower.tri(results$posterior_mean_indicator)] = indicator_means - results$posterior_mean_indicator <- results$posterior_mean_indicator + + results$posterior_mean_indicator = results$posterior_mean_indicator + t(results$posterior_mean_indicator) - if (has_sbm) { - results$posterior_mean_coclustering_matrix <- co_occur_matrix - results <- attach_sbm_posterior_summary(results, raw, build_arguments(spec)) + if(has_sbm) { + results$posterior_mean_coclustering_matrix = co_occur_matrix + results = attach_sbm_posterior_summary(results, raw, build_arguments(spec)) } } # --- arguments + class ------------------------------------------------------ - results$arguments <- build_arguments(spec) + results$arguments = build_arguments(spec) # Report the number of chains actually kept; failed chains are dropped # upstream, so raw holds only the survivors. - results$arguments$num_chains <- length(raw) - class(results) <- "bgms" + results$arguments$num_chains = length(raw) + class(results) = "bgms" # --- raw_samples ------------------------------------------------------------ # Allocation samples are per node (one column per variable), so they carry # node names, not edge names. - alloc_names <- if (identical(edge_prior, "Stochastic-Block")) { + alloc_names = if(identical(edge_prior, "Stochastic-Block")) { data_columnnames } else { NULL } - results$raw_samples <- build_raw_samples_list( + results$raw_samples = build_raw_samples_list( raw, edge_selection, edge_prior, names_main, edge_names, alloc_names ) @@ -278,28 +278,28 @@ build_output_bgm <- function(spec, raw) { # NULL placeholders ensure names(fit) lists these fields for easybgm compat. # The $.bgms method routes actual access through the cache. # Use list(NULL) because results$x = NULL removes the element in R. - results["posterior_summary_main"] <- list(NULL) - results["posterior_summary_pairwise"] <- list(NULL) - if (edge_selection) { - results["posterior_summary_indicator"] <- list(NULL) + results["posterior_summary_main"] = list(NULL) + results["posterior_summary_pairwise"] = list(NULL) + if(edge_selection) { + results["posterior_summary_indicator"] = list(NULL) } - results$cache <- cache + results$cache = cache # --- Sampler diagnostics ---------------------------------------------------- - results <- attach_sampler_diagnostics( + results = attach_sampler_diagnostics( results, raw, s$update_method, s$nuts_max_depth, names_main = names_main, names_pairwise = edge_names, target_accept = s$target_accept ) # --- Z-ratio alarm suite (hierarchical graph-prior spec) ---------------------- - if (!is.null(zratio_chains)) { - zc <- zratio_constants( + if(!is.null(zratio_chains)) { + zc = zratio_constants( delta = p$delta, sigma = 2 * p$pairwise_scale, beta = p$scale_rate / 2 ) - results$zratio_diag <- summarize_zratio_diagnostics( + results$zratio_diag = summarize_zratio_diagnostics( zratio_chains, zc, num_nodes = num_variables, seed = s$seed, @@ -307,8 +307,8 @@ build_output_bgm <- function(spec, raw) { ) } - results$.bgm_spec <- spec - if (needs_easybgm_s3_compat()) { + results$.bgm_spec = spec + if(needs_easybgm_s3_compat()) { results } else { s3_list_to_bgms(results) diff --git a/R/build_output_mixed_mrf.R b/R/build_output_mixed_mrf.R index 29a24330..0b42c0b3 100644 --- a/R/build_output_mixed_mrf.R +++ b/R/build_output_mixed_mrf.R @@ -17,42 +17,42 @@ # stored separately in posterior_mean_residual_variance; # the diagonal of the pairwise interaction matrix is zero. # ============================================================================== -build_output_mixed_mrf <- function(spec, raw) { - d <- spec$data - v <- spec$variables - pr <- spec$prior - s <- spec$sampler +build_output_mixed_mrf = function(spec, raw) { + d = spec$data + v = spec$variables + pr = spec$prior + s = spec$sampler - p <- d$num_discrete - q <- d$num_continuous - num_variables <- d$num_variables - data_columnnames <- d$data_columnnames - disc_names <- d$data_columnnames_discrete - cont_names <- d$data_columnnames_continuous - disc_idx <- d$discrete_indices - cont_idx <- d$continuous_indices - is_ordinal <- v$is_ordinal - num_categories <- d$num_categories - edge_selection <- pr$edge_selection + p = d$num_discrete + q = d$num_continuous + num_variables = d$num_variables + data_columnnames = d$data_columnnames + disc_names = d$data_columnnames_discrete + cont_names = d$data_columnnames_continuous + disc_idx = d$discrete_indices + cont_idx = d$continuous_indices + is_ordinal = v$is_ordinal + num_categories = d$num_categories + edge_selection = pr$edge_selection # 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$precision_graph_prior, "hierarchical")) { + zratio_chains = if(identical(pr$precision_graph_prior, "hierarchical")) { raw } else { NULL } # --- Compute index layout in flat parameter vector -------------------------- - layout <- compute_mixed_parameter_indices( + layout = compute_mixed_parameter_indices( num_thresholds = spec$precomputed$num_thresholds, p = p, q = q ) - nt <- layout$num_thresholds - main_idx <- layout$main_idx - pairwise_idx <- layout$pairwise_idx + nt = layout$num_thresholds + main_idx = layout$main_idx + pairwise_idx = layout$pairwise_idx # --- Indicator index layout ------------------------------------------------- # C++ indicator vector: [Gxx_ut | Gyy_ut | Gxy] @@ -60,63 +60,63 @@ build_output_mixed_mrf <- function(spec, raw) { # Discrete, continuous, cross edges --- same order as pairwise_idx above. # --- Normalize raw output per chain ----------------------------------------- - raw <- lapply(raw, function(chain) { - samples_t <- t(chain$samples) - res <- list( + raw = lapply(raw, function(chain) { + samples_t = t(chain$samples) + res = list( main_samples = samples_t[, main_idx, drop = FALSE], pairwise_samples = samples_t[, pairwise_idx, drop = FALSE], userInterrupt = isTRUE(chain$userInterrupt), chain_id = chain$chain_id ) - if (!is.null(chain$indicator_samples)) { - res$indicator_samples <- t(chain$indicator_samples) + if(!is.null(chain$indicator_samples)) { + res$indicator_samples = t(chain$indicator_samples) } - if (!is.null(chain$allocation_samples)) { - res$allocations <- t(chain$allocation_samples) + if(!is.null(chain$allocation_samples)) { + res$allocations = t(chain$allocation_samples) } - if (!is.null(chain$inclusion_parameter_samples)) { - res$inclusion_parameter <- as.numeric(chain$inclusion_parameter_samples) + if(!is.null(chain$inclusion_parameter_samples)) { + res$inclusion_parameter = as.numeric(chain$inclusion_parameter_samples) } attach_diagnostic_traces(res, chain) }) # --- Parameter names -------------------------------------------------------- # Main effect names (in internal order: discrete first, continuous second) - names_main <- character() - for (si in seq_len(p)) { - if (is_ordinal[si]) { - cats <- ordinal_threshold_labels( + names_main = character() + for(si in seq_len(p)) { + if(is_ordinal[si]) { + cats = ordinal_threshold_labels( num_categories[si], - if (!is.null(d$category_levels)) d$category_levels[[si]] else NULL + if(!is.null(d$category_levels)) d$category_levels[[si]] else NULL ) - names_main <- c(names_main, paste0(disc_names[si], " (", cats, ")")) + names_main = c(names_main, paste0(disc_names[si], " (", cats, ")")) } else { - names_main <- c( + names_main = c( names_main, paste0(disc_names[si], " (linear)"), paste0(disc_names[si], " (quadratic)") ) } } - for (ji in seq_len(q)) { - names_main <- c(names_main, paste0(cont_names[ji], " (mean)")) + for(ji in seq_len(q)) { + names_main = c(names_main, paste0(cont_names[ji], " (mean)")) } - for (ji in seq_len(q)) { - names_main <- c(names_main, paste0(cont_names[ji], " (precision diag)")) + for(ji in seq_len(q)) { + names_main = c(names_main, paste0(cont_names[ji], " (precision diag)")) } # Pairwise edge names --- internal order, mapped to original column names # We need a mapping from internal index to original variable name # Internal variables: [disc_1, ..., disc_p, cont_1, ..., cont_q] # Their original names: c(disc_names, cont_names) - all_internal_names <- c(disc_names, cont_names) + all_internal_names = c(disc_names, cont_names) - edge_names <- character() + edge_names = character() # Discrete-discrete edges - if (p > 1) { - for (i in seq_len(p - 1)) { - for (j in seq(i + 1, p)) { - edge_names <- c( + if(p > 1) { + for(i in seq_len(p - 1)) { + for(j in seq(i + 1, p)) { + edge_names = c( edge_names, paste0(disc_names[i], "-", disc_names[j]) ) @@ -124,10 +124,10 @@ build_output_mixed_mrf <- function(spec, raw) { } } # Continuous-continuous edges (off-diagonal) - if (q > 1) { - for (i in seq_len(q - 1)) { - for (j in seq(i + 1, q)) { - edge_names <- c( + if(q > 1) { + for(i in seq_len(q - 1)) { + for(j in seq(i + 1, q)) { + edge_names = c( edge_names, paste0(cont_names[i], "-", cont_names[j]) ) @@ -135,10 +135,10 @@ build_output_mixed_mrf <- function(spec, raw) { } } # Cross edges (discrete-continuous) - if (p > 0 && q > 0) { - for (i in seq_len(p)) { - for (j in seq_len(q)) { - edge_names <- c( + if(p > 0 && q > 0) { + for(i in seq_len(p)) { + for(j in seq_len(q)) { + edge_names = c( edge_names, paste0(disc_names[i], "-", cont_names[j]) ) @@ -147,139 +147,139 @@ build_output_mixed_mrf <- function(spec, raw) { } # --- Lazy MCMC diagnostics cache -------------------------------------------- - cache <- new.env(parent = emptyenv()) - cache$raw <- raw - cache$edge_selection <- edge_selection - cache$names_main <- names_main - cache$edge_names <- edge_names - cache$is_continuous <- FALSE - cache$model_type <- "mixed_mrf" - cache$summaries_computed <- FALSE + cache = new.env(parent = emptyenv()) + cache$raw = raw + cache$edge_selection = edge_selection + cache$names_main = names_main + cache$edge_names = edge_names + cache$is_continuous = FALSE + cache$model_type = "mixed_mrf" + cache$summaries_computed = FALSE # --- Compute posterior means from raw samples (cheap) ----------------------- - pooled_main <- do.call(rbind, lapply(raw, function(ch) ch$main_samples)) - pooled_pair <- do.call(rbind, lapply(raw, function(ch) ch$pairwise_samples)) - main_means <- colMeans(pooled_main) - pair_means <- colMeans(pooled_pair) + pooled_main = do.call(rbind, lapply(raw, function(ch) ch$main_samples)) + pooled_pair = do.call(rbind, lapply(raw, function(ch) ch$pairwise_samples)) + main_means = colMeans(pooled_main) + pair_means = colMeans(pooled_pair) # Split main_means into true main effects and quadratic (precision diagonal) - n_main <- nt + q # thresholds + continuous means - n_quad <- layout$num_quadratic # precision diagonal entries + n_main = nt + q # thresholds + continuous means + n_quad = layout$num_quadratic # precision diagonal entries - cache$n_main <- n_main - cache$n_quad <- n_quad + cache$n_main = n_main + cache$n_quad = n_quad - results <- list() + results = list() # --- Edge selection summaries (deferred) ------------------------------------ - edge_prior <- pr$edge_prior - has_sbm <- FALSE - if (edge_selection) { - has_sbm <- identical(edge_prior, "Stochastic-Block") && + edge_prior = pr$edge_prior + has_sbm = FALSE + if(edge_selection) { + has_sbm = identical(edge_prior, "Stochastic-Block") && "allocations" %in% names(raw[[1]]) - if (has_sbm) { - sbm_convergence <- summarize_alloc_pairs( + if(has_sbm) { + sbm_convergence = summarize_alloc_pairs( allocations = lapply(raw, `[[`, "allocations"), node_names = all_internal_names ) - results$posterior_summary_pairwise_allocations <- sbm_convergence$sbm_summary - co_occur_matrix <- sbm_convergence$co_occur_matrix + results$posterior_summary_pairwise_allocations = sbm_convergence$sbm_summary + co_occur_matrix = sbm_convergence$co_occur_matrix } - if ("inclusion_parameter" %in% names(raw[[1]])) { - results$inclusion_parameter_samples <- + if("inclusion_parameter" %in% names(raw[[1]])) { + results$inclusion_parameter_samples = lapply(raw, `[[`, "inclusion_parameter") } } # --- Posterior mean: main --------------------------------------------------- # Discrete main effects: p x max_cats matrix (like OMRF) - num_params_disc <- ifelse(is_ordinal, num_categories, 2L) - max_num_cats <- max(num_params_disc) - pmm_disc <- matrix(NA, nrow = p, ncol = max_num_cats) - start <- 0L - stop <- 0L - for (si in seq_len(p)) { - if (is_ordinal[si]) { - start <- stop + 1L - stop <- start + num_categories[si] - 1L - pmm_disc[si, seq_len(num_categories[si])] <- main_means[start:stop] + num_params_disc = ifelse(is_ordinal, num_categories, 2L) + max_num_cats = max(num_params_disc) + pmm_disc = matrix(NA, nrow = p, ncol = max_num_cats) + start = 0L + stop = 0L + for(si in seq_len(p)) { + if(is_ordinal[si]) { + start = stop + 1L + stop = start + num_categories[si] - 1L + pmm_disc[si, seq_len(num_categories[si])] = main_means[start:stop] } else { - start <- stop + 1L - stop <- start + 1L - pmm_disc[si, 1:2] <- main_means[start:stop] + start = stop + 1L + stop = start + 1L + pmm_disc[si, 1:2] = main_means[start:stop] } } - rownames(pmm_disc) <- disc_names - colnames(pmm_disc) <- paste0("cat (", seq_len(max_num_cats), ")") + rownames(pmm_disc) = disc_names + colnames(pmm_disc) = paste0("cat (", seq_len(max_num_cats), ")") # Continuous main effects: q x 1 matrix (means only) - pmm_cont <- matrix(main_means[nt + seq_len(q)], + pmm_cont = matrix(main_means[nt + seq_len(q)], nrow = q, ncol = 1, dimnames = list(cont_names, "mean") ) - results$posterior_mean_main <- list( + results$posterior_mean_main = list( discrete = pmm_disc, continuous = pmm_cont ) # --- Posterior mean: associations (all blocks) ----------------------------- - dn <- list(data_columnnames, data_columnnames) - results$posterior_mean_pairwise <- fill_mixed_symmetric( + dn = list(data_columnnames, data_columnnames) + results$posterior_mean_pairwise = fill_mixed_symmetric( pair_means, p, q, disc_idx, cont_idx, dn ) # --- Residual variance (continuous diagonal) -------------------------------- # C++ stores negative association diagonal; residual variance = -1/(2*diag). # Average per-sample inversions to avoid Jensen's inequality bias. - pooled_quad <- pooled_main[, nt + q + seq_len(q), drop = FALSE] - rv <- colMeans(-1 / (2 * pooled_quad)) - names(rv) <- cont_names - results$posterior_mean_residual_variance <- rv + pooled_quad = pooled_main[, nt + q + seq_len(q), drop = FALSE] + rv = colMeans(-1 / (2 * pooled_quad)) + names(rv) = cont_names + results$posterior_mean_residual_variance = rv # --- Posterior mean: indicator ----------------------------------------------- - if (edge_selection) { - pooled_ind <- do.call(rbind, lapply(raw, function(ch) ch$indicator_samples)) - indicator_means <- colMeans(pooled_ind) - results$posterior_mean_indicator <- fill_mixed_symmetric( + if(edge_selection) { + pooled_ind = do.call(rbind, lapply(raw, function(ch) ch$indicator_samples)) + indicator_means = colMeans(pooled_ind) + results$posterior_mean_indicator = fill_mixed_symmetric( indicator_means, p, q, disc_idx, cont_idx, dn ) - if (has_sbm) { - results$posterior_mean_coclustering_matrix <- co_occur_matrix - results <- attach_sbm_posterior_summary(results, raw, build_arguments(spec)) + if(has_sbm) { + results$posterior_mean_coclustering_matrix = co_occur_matrix + results = attach_sbm_posterior_summary(results, raw, build_arguments(spec)) } } # --- arguments + class ------------------------------------------------------ - results$arguments <- build_arguments(spec) + results$arguments = build_arguments(spec) # Report the number of chains actually kept; failed chains are dropped # upstream, so raw holds only the survivors. - results$arguments$num_chains <- length(raw) + results$arguments$num_chains = length(raw) # NULL placeholders ensure names(fit) lists these fields for easybgm compat. # Use list(NULL) because results$x = NULL removes the element in R. - results["posterior_summary_main"] <- list(NULL) - results["posterior_summary_pairwise"] <- list(NULL) - if (edge_selection) { - results["posterior_summary_indicator"] <- list(NULL) + results["posterior_summary_main"] = list(NULL) + results["posterior_summary_pairwise"] = list(NULL) + if(edge_selection) { + results["posterior_summary_indicator"] = list(NULL) } - results$cache <- cache - class(results) <- "bgms" + results$cache = cache + class(results) = "bgms" # --- raw_samples ------------------------------------------------------------ - alloc_names <- if (identical(edge_prior, "Stochastic-Block")) { + alloc_names = if(identical(edge_prior, "Stochastic-Block")) { all_internal_names } else { NULL } - results$raw_samples <- build_raw_samples_list( + results$raw_samples = build_raw_samples_list( raw, edge_selection, edge_prior, names_main, edge_names, alloc_names ) # --- Sampler diagnostics ---------------------------------------------------- - results <- attach_sampler_diagnostics( + results = attach_sampler_diagnostics( results, raw, s$update_method, s$nuts_max_depth, names_main = names_main, names_pairwise = edge_names, target_accept = s$target_accept @@ -288,13 +288,13 @@ build_output_mixed_mrf <- function(spec, raw) { # --- Z-ratio alarm suite (hierarchical spec on the continuous block) --------- # The indicator vector is [Gxx upper, Gyy upper, Gxy row-major], both # triangles without diagonals; the audit reads the Gyy segment. - if (!is.null(zratio_chains)) { - zc <- zratio_constants( + if(!is.null(zratio_chains)) { + zc = zratio_constants( delta = pr$delta, sigma = 2 * pr$pairwise_scale, beta = pr$scale_rate / 2 ) - results$zratio_diag <- summarize_zratio_diagnostics( + results$zratio_diag = summarize_zratio_diagnostics( zratio_chains, zc, num_nodes = q, seed = s$seed, @@ -304,8 +304,8 @@ build_output_mixed_mrf <- function(spec, raw) { ) } - results$.bgm_spec <- spec - if (needs_easybgm_s3_compat()) { + results$.bgm_spec = spec + if(needs_easybgm_s3_compat()) { results } else { s3_list_to_bgms(results) diff --git a/R/build_spec.R b/R/build_spec.R index e986459c..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,35 +50,34 @@ sampler_sublist = function(s) { } - -build_spec_ggm <- function(x, data_columnnames, num_variables, - variable_type, is_ordinal, is_continuous, - baseline_category, - na_action, sampler, - interaction_prior_type, pairwise_scale, - interaction_alpha, interaction_beta, - scale_prior_type, scale_shape, scale_rate, - scale_eta = NA_real_, - delta = 0, - precision_graph_prior = "joint", - calibration_window = NULL, - edge_prior_flat) { +build_spec_ggm = function(x, data_columnnames, num_variables, + variable_type, is_ordinal, is_continuous, + baseline_category, + na_action, sampler, + interaction_prior_type, pairwise_scale, + interaction_alpha, interaction_beta, + scale_prior_type, scale_shape, scale_rate, + scale_eta = NA_real_, + delta = 0, + precision_graph_prior = "joint", + calibration_window = NULL, + edge_prior_flat) { # Missing data - md <- validate_missing_data( + md = validate_missing_data( x = x, na_action = na_action, is_continuous = TRUE ) - x <- md$x + x = md$x # Center continuous data (GGM likelihood assumes zero mean). The column # means are kept so prediction can center newdata on the training scale. - column_means <- colMeans(x) - x <- center_continuous_data(x) + column_means = colMeans(x) + x = center_continuous_data(x) # Standardized-frame scale prior: derive the raw diagonal rate eta / s - scale_rate <- resolve_scale_rate(scale_rate, scale_eta, pairwise_scale) + scale_rate = resolve_scale_rate(scale_rate, scale_eta, pairwise_scale) - ep <- edge_prior_flat + ep = edge_prior_flat new_bgm_spec( model_type = "ggm", @@ -215,53 +213,53 @@ build_spec_omrf = function(x, data_columnnames, num_variables, # the spec with metadata needed by sample_mixed_mrf() and # build_output_mixed_mrf(). # ------------------------------------------------------------------ -build_spec_mixed_mrf <- function(x, data_columnnames, num_variables, - variable_type, is_ordinal, - baseline_category, - na_action, sampler, - interaction_prior_type, pairwise_scale, - interaction_alpha, interaction_beta, - threshold_prior_type, main_alpha, main_beta, - threshold_scale, - means_prior_type, means_scale, - means_alpha, means_beta, - scale_prior_type, scale_shape, scale_rate, - scale_eta = NA_real_, - delta = 0, - precision_graph_prior = "joint", - calibration_window = NULL, - edge_prior_flat) { +build_spec_mixed_mrf = function(x, data_columnnames, num_variables, + variable_type, is_ordinal, + baseline_category, + na_action, sampler, + interaction_prior_type, pairwise_scale, + interaction_alpha, interaction_beta, + threshold_prior_type, main_alpha, main_beta, + threshold_scale, + means_prior_type, means_scale, + means_alpha, means_beta, + scale_prior_type, scale_shape, scale_rate, + scale_eta = NA_real_, + delta = 0, + precision_graph_prior = "joint", + calibration_window = NULL, + edge_prior_flat) { # Standardized-frame scale prior: derive the raw diagonal rate eta / s - scale_rate <- resolve_scale_rate(scale_rate, scale_eta, pairwise_scale) + scale_rate = resolve_scale_rate(scale_rate, scale_eta, pairwise_scale) # Identify discrete vs continuous columns - cont_idx <- which(variable_type == "continuous") - disc_idx <- which(variable_type != "continuous") - p <- length(disc_idx) - q <- length(cont_idx) + cont_idx = which(variable_type == "continuous") + disc_idx = which(variable_type != "continuous") + p = length(disc_idx) + q = length(cont_idx) # Split data - x_disc <- x[, disc_idx, drop = FALSE] - x_cont <- x[, cont_idx, drop = FALSE] + x_disc = x[, disc_idx, drop = FALSE] + x_cont = x[, cont_idx, drop = FALSE] # Ensure integer matrix for discrete data - storage.mode(x_disc) <- "integer" + storage.mode(x_disc) = "integer" # Ensure numeric matrix for continuous data - storage.mode(x_cont) <- "double" + storage.mode(x_cont) = "double" # Discrete variable properties (subset to discrete columns) - is_ordinal_disc <- is_ordinal[disc_idx] - vtype_disc <- variable_type[disc_idx] + is_ordinal_disc = is_ordinal[disc_idx] + vtype_disc = variable_type[disc_idx] # Subset baseline_category to discrete columns when the user supplies a # full-length vector (one entry per variable, including continuous ones). - if (length(baseline_category) == num_variables && num_variables != p) { - baseline_category <- baseline_category[disc_idx] + if(length(baseline_category) == num_variables && num_variables != p) { + baseline_category = baseline_category[disc_idx] } # Baseline category for discrete variables - bc <- validate_baseline_category( + bc = validate_baseline_category( baseline_category = baseline_category, baseline_category_provided = !is.null(baseline_category), x = x_disc, @@ -269,30 +267,30 @@ build_spec_mixed_mrf <- function(x, data_columnnames, num_variables, ) # Missing data handling - na_impute <- FALSE - missing_index_discrete <- NULL - missing_index_continuous <- NULL + na_impute = FALSE + missing_index_discrete = NULL + missing_index_continuous = NULL - if (na_action == "listwise") { - missing_rows <- apply(x_disc, 1, anyNA) | apply(x_cont, 1, anyNA) - if (all(missing_rows)) { + if(na_action == "listwise") { + missing_rows = apply(x_disc, 1, anyNA) | apply(x_cont, 1, anyNA) + if(all(missing_rows)) { stop(paste0( "All rows in x contain at least one missing response.\n", "You could try option na_action = \"impute\"." )) } - n_removed <- sum(missing_rows) - if (n_removed > 0 && isTRUE(getOption("bgms.verbose", TRUE))) { - n_remaining <- nrow(x_disc) - n_removed + n_removed = sum(missing_rows) + if(n_removed > 0 && isTRUE(getOption("bgms.verbose", TRUE))) { + n_remaining = nrow(x_disc) - n_removed message( - n_removed, " row", if (n_removed > 1) "s" else "", + n_removed, " row", if(n_removed > 1) "s" else "", " with missing values excluded (n = ", n_remaining, " remaining).\n", "To impute missing values instead, use na_action = \"impute\"." ) } - x_disc <- x_disc[!missing_rows, , drop = FALSE] - x_cont <- x_cont[!missing_rows, , drop = FALSE] - if (nrow(x_disc) < 2) { + x_disc = x_disc[!missing_rows, , drop = FALSE] + x_cont = x_cont[!missing_rows, , drop = FALSE] + if(nrow(x_disc) < 2) { stop(paste0( "After removing missing observations from the input matrix x,\n", "there were less than two rows left in x." @@ -300,27 +298,27 @@ build_spec_mixed_mrf <- function(x, data_columnnames, num_variables, } } else { # Impute path: handle discrete and continuous sub-matrices separately - md_disc <- handle_impute(x_disc) - md_cont <- handle_impute(x_cont) - x_disc <- md_disc$x - x_cont <- md_cont$x - na_impute <- md_disc$na_impute || md_cont$na_impute - if (md_disc$na_impute) missing_index_discrete <- md_disc$missing_index - if (md_cont$na_impute) missing_index_continuous <- md_cont$missing_index + md_disc = handle_impute(x_disc) + md_cont = handle_impute(x_cont) + x_disc = md_disc$x + x_cont = md_cont$x + na_impute = md_disc$na_impute || md_cont$na_impute + if(md_disc$na_impute) missing_index_discrete = md_disc$missing_index + if(md_cont$na_impute) missing_index_continuous = md_cont$missing_index } # Ordinal recoding (reformat discrete data) - ord <- reformat_ordinal_data( + ord = reformat_ordinal_data( x = x_disc, is_ordinal = is_ordinal_disc, baseline_category = bc ) - x_disc_recoded <- ord$x - num_categories <- ord$num_categories - bc_final <- ord$baseline_category + x_disc_recoded = ord$x + num_categories = ord$num_categories + bc_final = ord$baseline_category - ep <- edge_prior_flat + ep = edge_prior_flat - num_thresholds <- sum(ifelse(is_ordinal_disc, num_categories, 2L)) + num_thresholds = sum(ifelse(is_ordinal_disc, num_categories, 2L)) new_bgm_spec( model_type = "mixed_mrf", diff --git a/R/run_sampler.R b/R/run_sampler.R index 673b68b6..f48cd716 100644 --- a/R/run_sampler.R +++ b/R/run_sampler.R @@ -61,29 +61,29 @@ run_sampler = function(spec) { # ============================================================================== # run_sampler_ggm() # ============================================================================== -run_sampler_ggm <- function(spec) { - d <- spec$data - p <- spec$prior - s <- spec$sampler - m <- spec$missing +run_sampler_ggm = function(spec) { + d = spec$data + p = spec$prior + s = spec$sampler + m = spec$missing - bb_alpha_between <- bb_between_or_sentinel(p$beta_bernoulli_alpha_between) - bb_beta_between <- bb_between_or_sentinel(p$beta_bernoulli_beta_between) + bb_alpha_between = bb_between_or_sentinel(p$beta_bernoulli_alpha_between) + bb_beta_between = bb_between_or_sentinel(p$beta_bernoulli_beta_between) # Graph-prior specification: the joint path corrects the hyperparameter # updates for the untracked normalizer Z(Gamma); the hierarchical path # tracks Z(Gamma) itself in the between-edge moves (per-edge Z-ratio # engine) and keeps the hyperparameter updates clean conjugate. The two # are mutually exclusive. - correction <- NULL - zratio <- NULL - if (identical(p$precision_graph_prior, "hierarchical")) { - zc <- zratio_constants( + correction = NULL + zratio = NULL + if(identical(p$precision_graph_prior, "hierarchical")) { + zc = zratio_constants( delta = p$delta, sigma = 2 * p$pairwise_scale, beta = p$scale_rate / 2 ) - zratio <- list( + zratio = list( addc = zc$addc, tg = zc$tg, ihat = zc$ihat, ghat = zc$ghat, wt = zc$wt, psi0 = zc$psi0, delta = zc$delta, sigma = zc$sigma, beta = zc$beta, @@ -92,10 +92,10 @@ run_sampler_ggm <- function(spec) { ) ) } else { - correction <- ggm_edge_prior_correction(p, s, d$num_variables) + correction = ggm_edge_prior_correction(p, s, d$num_variables) } - out_raw <- sample_ggm( + out_raw = sample_ggm( inputFromR = list( X = d$x, pairwise_scale = p$pairwise_scale, @@ -206,28 +206,28 @@ run_sampler_omrf = function(spec) { # ============================================================================== # run_sampler_mixed_mrf() # ============================================================================== -run_sampler_mixed_mrf <- function(spec) { - d <- spec$data - v <- spec$variables - m <- spec$missing - p <- spec$prior - s <- spec$sampler +run_sampler_mixed_mrf = function(spec) { + d = spec$data + v = spec$variables + m = spec$missing + p = spec$prior + s = spec$sampler - bb_alpha_between <- bb_between_or_sentinel(p$beta_bernoulli_alpha_between) - bb_beta_between <- bb_between_or_sentinel(p$beta_bernoulli_beta_between) + bb_alpha_between = bb_between_or_sentinel(p$beta_bernoulli_alpha_between) + bb_beta_between = bb_between_or_sentinel(p$beta_bernoulli_beta_between) # Graph-prior specification on the continuous block: same dichotomy as the # GGM path (see run_sampler_ggm). The Z-ratio constants and the Stage-3d # window are sized on the continuous subgraph. - correction <- NULL - zratio <- NULL - if (identical(p$precision_graph_prior, "hierarchical")) { - zc <- zratio_constants( + correction = NULL + zratio = NULL + if(identical(p$precision_graph_prior, "hierarchical")) { + zc = zratio_constants( delta = p$delta, sigma = 2 * p$pairwise_scale, beta = p$scale_rate / 2 ) - zratio <- list( + zratio = list( addc = zc$addc, tg = zc$tg, ihat = zc$ihat, ghat = zc$ghat, wt = zc$wt, psi0 = zc$psi0, delta = zc$delta, sigma = zc$sigma, beta = zc$beta, @@ -236,12 +236,12 @@ run_sampler_mixed_mrf <- function(spec) { ) ) } else { - correction <- ggm_edge_prior_correction( + correction = ggm_edge_prior_correction( p, s, d$num_variables, d$num_continuous ) } - input_list <- list( + input_list = list( discrete_observations = d$x_discrete, continuous_observations = d$x_continuous, num_categories = d$num_categories, @@ -264,7 +264,7 @@ run_sampler_mixed_mrf <- function(spec) { scale_rate = p$scale_rate ) - out_raw <- sample_mixed_mrf( + out_raw = sample_mixed_mrf( inputFromR = input_list, prior_inclusion_prob = p$inclusion_probability, initial_edge_indicators = matrix(1L, diff --git a/tests/testthat/test-bgm-hier-spec.R b/tests/testthat/test-bgm-hier-spec.R index b4363fb6..4c325fca 100644 --- a/tests/testthat/test-bgm-hier-spec.R +++ b/tests/testthat/test-bgm-hier-spec.R @@ -7,7 +7,7 @@ hier_test_data = function(q = 10, n = 40, seed = 4) { } test_that("hierarchical spec eligibility is validated", { - Y <- hier_test_data() + Y = hier_test_data() expect_error( bgm( @@ -52,9 +52,9 @@ test_that("hierarchical spec eligibility is validated", { 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( + 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), @@ -63,7 +63,7 @@ test_that("the hierarchical spec accepts a Cauchy slab on every update method", update_method = method, chains = 1, cores = 1, seed = 7, display_progress = "none", verbose = FALSE ) - s <- summary(fit) + 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 @@ -75,8 +75,8 @@ test_that("the hierarchical spec accepts a Cauchy slab on every update method", test_that("bgm fits the hierarchical spec and attaches the alarm suite", { skip_on_cran() - Y <- hier_test_data(q = 12) - fit <- bgm( + Y = hier_test_data(q = 12) + fit = bgm( x = Y, variable_type = "continuous", iter = 150, warmup = 250, interaction_prior = normal_prior(scale = 0.5), @@ -86,7 +86,7 @@ test_that("bgm fits the hierarchical spec and attaches the alarm suite", { update_method = "gibbs", chains = 2, cores = 2, seed = 11, display_progress = "none", verbose = FALSE ) - zd <- fit@zratio_diag + zd = fit@zratio_diag expect_false(is.null(zd)) expect_equal(nrow(zd$per_chain), 2L) expect_true(all(zd$per_chain$frozen)) @@ -101,8 +101,8 @@ test_that("bgm fits the hierarchical spec and attaches the alarm suite", { test_that("the joint default is unchanged", { skip_on_cran() - Y <- hier_test_data(q = 6) - fit <- bgm( + Y = hier_test_data(q = 6) + fit = bgm( x = Y, variable_type = "continuous", iter = 100, warmup = 150, update_method = "gibbs", chains = 1, cores = 1, seed = 3, @@ -135,13 +135,13 @@ test_that("mixed indicator layout rebuilds the continuous subgraph", { test_that("mixed data supports the hierarchical spec on the continuous block", { skip_on_cran() set.seed(9) - n <- 60 - X <- cbind( + n = 60 + X = cbind( matrix(sample(0:2, n * 3, replace = TRUE), n, 3), matrix(rnorm(n * 8), n, 8) ) - colnames(X) <- paste0("V", seq_len(11)) - vt <- c(rep("ordinal", 3), rep("continuous", 8)) + colnames(X) = paste0("V", seq_len(11)) + vt = c(rep("ordinal", 3), rep("continuous", 8)) expect_error( bgm( @@ -154,7 +154,7 @@ test_that("mixed data supports the hierarchical spec on the continuous block", { "two continuous" ) - fit <- bgm( + fit = bgm( x = X, variable_type = vt, iter = 120, warmup = 200, interaction_prior = normal_prior(scale = 0.5), @@ -164,7 +164,7 @@ test_that("mixed data supports the hierarchical spec on the continuous block", { update_method = "adaptive-metropolis", chains = 1, cores = 1, seed = 5, display_progress = "none", verbose = FALSE ) - zd <- fit@zratio_diag + zd = fit@zratio_diag expect_false(is.null(zd)) expect_true(zd$per_chain$frozen) expect_false(zd$verdict_flagged) diff --git a/tests/testthat/test-build-arguments.R b/tests/testthat/test-build-arguments.R index cddb82fa..56a6a0b0 100644 --- a/tests/testthat/test-build-arguments.R +++ b/tests/testthat/test-build-arguments.R @@ -89,9 +89,9 @@ spec_compare = function(...) { # ============================================================================== test_that("GGM build_arguments: all expected field names present", { - s <- spec_ggm() - a <- build_arguments(s) - expected <- c( + s = spec_ggm() + a = build_arguments(s) + expected = c( "num_variables", "num_cases", "na_impute", "variable_type", "iter", "warmup", "edge_selection", "edge_prior", "precision_graph_prior", "calibration_window", From 59cb1cc203d085f244d3cf625d54e3750d5686c8 Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Sun, 5 Jul 2026 21:12:37 +0200 Subject: [PATCH 4/5] Drop the internal-representation detail from the hierarchical slab message The eligibility error for an unsupported interaction prior explained that the normalizer is derived for the normal slab and that the Cauchy slab enters through its scale mixture of normals. That describes the internal representation rather than what the user should do, so state only that a normal or Cauchy slab is supported. Same change in bgm(), sample_ggm_prior(), and the bgm() documentation. --- R/bgm.R | 3 +-- R/bgm_spec.R | 4 +--- R/sample_ggm_prior.R | 6 ++---- man/bgm.Rd | 3 +-- 4 files changed, 5 insertions(+), 11 deletions(-) diff --git a/R/bgm.R b/R/bgm.R index 9b04e73e..92f5d6de 100644 --- a/R/bgm.R +++ b/R/bgm.R @@ -189,8 +189,7 @@ #' 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()} or \code{cauchy_prior()} interaction prior -#' (the Cauchy slab enters through its scale mixture of normals), a +#' \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 diff --git a/R/bgm_spec.R b/R/bgm_spec.R index 0c73f475..13d65796 100644 --- a/R/bgm_spec.R +++ b/R/bgm_spec.R @@ -407,9 +407,7 @@ bgm_spec = function(x, stop(sprintf( paste0( "precision_graph_prior = \"hierarchical\" supports a normal or Cauchy ", - "interaction (slab) prior: the Z-ratio normalizer is derived for ", - "the normal slab, and the Cauchy slab enters through its scale ", - "mixture of normals. Got %s_prior(). Use interaction_prior = ", + "interaction (slab) prior. Got %s_prior(). Use interaction_prior = ", "normal_prior() or cauchy_prior(), or keep precision_graph_prior = ", "\"joint\"." ), diff --git a/R/sample_ggm_prior.R b/R/sample_ggm_prior.R index e550ea80..18cffb4e 100644 --- a/R/sample_ggm_prior.R +++ b/R/sample_ggm_prior.R @@ -319,10 +319,8 @@ sample_ggm_prior = function( stop(sprintf( paste0( "spec = \"hierarchical\" supports a normal or Cauchy interaction ", - "(slab) prior: the Z-ratio normalizer is derived for the normal ", - "slab, and the Cauchy slab enters through its scale mixture of ", - "normals. Got %s_prior(). Use interaction_prior = normal_prior() ", - "or cauchy_prior()." + "(slab) prior. Got %s_prior(). Use interaction_prior = ", + "normal_prior() or cauchy_prior()." ), ip$interaction_prior_type )) diff --git a/man/bgm.Rd b/man/bgm.Rd index c120096a..9b9ea71f 100644 --- a/man/bgm.Rd +++ b/man/bgm.Rd @@ -202,8 +202,7 @@ 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()} or \code{cauchy_prior()} interaction prior -(the Cauchy slab enters through its scale mixture of normals), a +\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 From e36705edaffbf23ce107daca028b1ce5fcc8f3f8 Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Sun, 5 Jul 2026 21:39:24 +0200 Subject: [PATCH 5/5] Soften the NUTS warmup diagnostic to match the other issue lines The energy-stationarity warmup check is a heuristic, and its message went straight to "increase warmup". The lower-severity divergence and tree-depth lines instead point the user at R-hat and ESS. Align the warmup line with them: report that warmup may be incomplete and direct the user to check R-hat and ESS rather than prescribing a change. --- R/diagnostics_nuts.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 = ", ") ))