R Scaffolding

The R layer sits between user-facing functions (bgm(), bgmCompare()) and the C++ sampler. It validates inputs, constructs a specification object, dispatches to the appropriate C++ sampler, and assembles the returned output into an S3 fit object.

Data flow

flowchart TD
    A["<b>bgm(x, ...)</b>"]
    B["<b>bgm_spec()</b><br/>validates, assembles sub-lists"]
    C["<b>new_bgm_spec()</b><br/>asserts types and field presence"]
    D["<b>validate_bgm_spec()</b><br/>cross-field invariant checks"]
    E["<b>run_sampler(spec)</b><br/>dispatches to C++ via model_type"]
    F["<b>build_output(spec, raw)</b><br/>normalizes raw C++ output to S3"]
    G["<b>bgm fit object</b>"]
    A --> B --> C --> D --> E --> F --> G

Both bgm() and bgmCompare() follow this pipeline. The bgmCompare() path adds group-specific preprocessing (projection matrices, group indices, precomputed sufficient statistics) before the spec is constructed.

Validation pipeline

Validation is split across four files, each with a focused scope.

Data checks (validate_data.R)

data_check() is the entry point. It enforces:

  • Input is a data frame or matrix (converts to integer matrix for ordinal data)
  • No constant columns
  • At least two variables
  • Column names exist and are unique
  • For ordinal data: all values are non-negative integers

Missing data handling branches on na_action:

  • "listwise" — rows with any NA are dropped
  • "impute" — missing indices are recorded for Bayesian imputation in C++

Variable type detection (validate_model.R)

validate_variable_types() classifies each column as "ordinal" or "blume-capel" (for discrete variables) based on the variable_type argument. When variable_type = "ordinal" (the default), all discrete variables are treated as ordinal. When a named vector is provided, each variable is classified individually.

validate_baseline_category() sets the reference category for Blume-Capel variables. It is required whenever any variable is Blume-Capel: the function errors if baseline_category was not supplied. A single integer is recycled across variables. When every variable is ordinal the function returns zeros, which nothing reads.

Sampler validation (validate_sampler.R)

validate_sampler() resolves the update_method argument into a concrete sampler type ("nuts", "adaptive-metropolis", or "gibbs") and sets defaults for tuning parameters (target_accept, nuts_max_depth). The Gibbs sampler is gated here on data type only: it requires all-continuous data. The prior-family requirement (a Normal or Cauchy interaction prior with a Gamma-family scale prior) is checked C++-side instead, in sample_ggm.cpp via GGMModel::row_block_gibbs_eligible(). Its target_accept is recorded as NA_real_ — an exact draw has no acceptance target, and the C++ side ignores the value for Gibbs. The function also issues the progressive short-warmup warnings (see Warmup Schedule), detects the number of available cores, and resolves the progress display type.

Prior validation (R/priors.R, validate_model.R)

Priors are user-facing S3 objects (bgms_parameter_prior, bgms_scale_prior, bgms_indicator_prior). Each prior constructor (cauchy_prior(), gamma_prior(), bernoulli_prior(), …) validates its own hyperparameters at the call site, so by the time the spec is built every prior object is already self-consistent.

unpack_parameter_prior(), unpack_scale_prior(), unpack_interaction_prior(), unpack_threshold_prior(), and unpack_indicator_prior() flatten each prior object into the (family, hyperparameters) representation that the C++ bridge expects. validate_edge_prior() and validate_difference_prior() accept either a prior object or a legacy character string ("Bernoulli", "Beta-Bernoulli", "Stochastic-Block"); strings are forwarded to the matching constructor with a lifecycle warning.

Spec construction (bgm_spec.R)

bgm_spec() is the internal constructor. It calls all validators, then delegates to one of four model-specific builders:

Model type Builder Key additions
"ggm" build_spec_ggm() Sufficient statistics (X'X), precision scale
"omrf" build_spec_omrf() Category counts, ordinal/BC flags, scaling factors
"mixed_mrf" build_spec_mixed_mrf() Separate discrete/continuous matrices
"compare" build_spec_compare() Group indices, projection matrix, precomputed pairwise stats

Each builder assembles a spec with seven components:

  • $model_type — model class string ("ggm", "omrf", etc.)
  • $data — observations, dimensions, variable names
  • $variables — variable types, ordinal flags, baseline categories
  • $missing — na_action, imputation flag, missing indices
  • $prior — prior type, hyperparameters, scaling factors
  • $sampler — algorithm, iterations, warmup, chains, seed
  • $precomputed — sufficient statistics and derived quantities (e.g., cross-products for GGM, pairwise stats for compare)

The result passes through new_bgm_spec() (type assertions) and validate_bgm_spec() (cross-field invariants, such as: if edge_selection = TRUE then edge_prior must not be "Not Applicable").

The arguments record (build_arguments.R)

The spec is an internal object: it carries the data matrices, the precomputed sufficient statistics, and everything else the sampler needs, and it is not what downstream code should be reading. So one step sits between the validated spec and the fit object, and it is the step that decides what the rest of the package can see.

build_arguments(spec) produces the $arguments list every fit carries. It is a pure function of the spec, so it belongs to this stage conceptually, but it runs later: each of the three output builders in the next section calls it while assembling the fit object. It dispatches on spec$model_type to one of four builders, mirroring the four spec builders above:

Model type Builder What only it carries
"ggm" build_arguments_ggm() column_means (the training means continuous data was centred on), is_continuous = TRUE
"omrf" build_arguments_omrf() category_levels, blume_capel_shift, baseline_category, nuts_max_depth
"mixed_mrf" build_arguments_mixed_mrf() discrete_indices / continuous_indices, num_discrete / num_continuous, is_mixed = TRUE
"compare" build_arguments_compare() projection, num_groups, group_labels, category_support

Three things about the result are worth knowing, because they explain behaviour elsewhere.

It is the model-class oracle. Nothing downstream re-derives the model type from the data. extract_arguments(fit)$is_continuous and $is_mixed are the flags that simulate() and predict() branch on, and model_type is what print() turns into a label.

It carries the recode maps. category_levels and blume_capel_shift are the only record of how the fitted data were mapped onto zero-based codes, which is what lets predict() accept newdata on the original scale and simulate() return data on it.

It is a compatibility surface, not a clean design. Several fields are duplicated or legacy: no_variables repeats num_variables under the pre-0.1.6.0 name, and version records the package version that produced the fit so that later readers can branch on vintage. Both are kept because reading a saved fit has to keep working.

The ordinal builder also collapses variable_type to a single string when every variable shares one, which is why extract_arguments(fit)$variable_type can be length one or length p and callers pass it through expand_variable_type() before indexing it.

Sampler dispatch (run_sampler.R)

run_sampler() reads spec$model_type and calls the corresponding C++ entry point:

model_type C++ function Entry file
"ggm" sample_ggm() src/sample_ggm.cpp
"omrf" sample_omrf() src/sample_omrf.cpp
"mixed_mrf" sample_mixed_mrf() src/sample_mixed.cpp
"compare" run_bgmCompare_parallel() src/bgmCompare_interface.cpp

Each C++ function receives an R list of arguments, constructs model and prior objects, and calls run_mcmc_sampler() from the chain runner. The return value is a list of per-chain raw output.

For continuous and mixed models with edge selection, run_sampler() also prepares the graph-prior machinery before the C++ call, and the two paths are mutually exclusive:

  • Under the default precision_graph_prior = "hierarchical", it assembles the Z-ratio spec: the analytic cell constants (zratio_cell_constants(), from delta, the slab family and scale, and the diagonal prior’s rate frame), the Option-B absolute-moment surface built once for the analysis (zratio_build_surfaces(); if the build fails or the diagonal’s Gamma shape is non-unit, the engine falls back to the coarser additive path with a message), and the trust-gauge option (options(bgms.zratio_gauge_sweeps = ...)).
  • Under the "joint" composition with a Beta-Bernoulli or SBM edge prior, it builds or loads the normalizing-constant correction table (ggm_edge_prior_correction()), cached on disk via tools::R_user_dir("bgms", "cache") so later fits of the same configuration skip the build.

Output assembly (build_output.R)

build_output() transforms raw C++ output into the S3 objects returned to the user ("bgms" or "bgmCompare" class).

GGM and OMRF (build_output_bgm())

The GGM and OMRF paths share a unified builder. Key operations:

  1. Normalize raw samples — Split the flat parameter vector into main-effect and pairwise-interaction matrices per chain
  2. Compute posterior means — Average across iterations and chains
  3. Compute inclusion probabilities — Average indicator samples (when edge selection is active)
  4. Assemble $raw_samples — Per-chain lists of main, pairwise, indicator, and allocations matrices; under edge selection, also the Rao-Blackwellized inclusion draws and odds accumulators (rb_inclusion, rb_counts — see GGM Internals), mapped with the same transpose and off-diagonal ordering as the indicator draws; with a Beta-Bernoulli edge prior on a continuous or mixed model, also the per-chain sampled inclusion probability (surfaced as fit$inclusion_parameter_samples)
  5. Sampler diagnostics — Under NUTS, summarize_nuts_diagnostics() builds fit$nuts_diag from the tree-depth, divergence, energy, and acceptance traces, including the R-side energy-based warmup completeness check. Under adaptive Metropolis, summarize_am_diagnostics() builds fit$am_diag from the stored samples and the acceptance trace
  6. Z-ratio diagnostics — Under the hierarchical graph prior with the trust gauge enabled, summarize_zratio_gauge() builds fit$zratio_diag from the per-chain gauge blocks; flagged issues print alongside the other sampler warnings

Mixed MRF (build_output_mixed_mrf())

The mixed MRF builder handles the block structure: discrete-discrete, continuous-continuous, and cross-type interactions are stored in separate blocks by C++ and need to be mapped back to the original variable ordering.

Compare (build_output_compare())

The comparison builder splits posterior means into group-level ($group_posterior_means) and contrast ($difference_posterior_means) components, and computes contrast inclusion probabilities.