bgm()

Fits a single-network Bayesian graphical model with bgm() and returns an S7 object of class bgms.

Description

bgm() fits a single network for binary, ordinal, continuous, or mixed data, with optional Bayesian edge selection. This page is the technical reference for function arguments and returned objects. For interpretation and workflow, start with Getting Started, then see Graphical Models and The Bayesian Approach. Methodological background is provided in Marsman et al. (2025).

Usage

bgm(
  x,
  variable_type = "ordinal",
  baseline_category,
  iter = 2000,
  warmup = 2000,
  interaction_prior = normal_prior(scale = 1),
  threshold_prior = beta_prime_prior(alpha = 0.5, beta = 0.5),
  means_prior = normal_prior(scale = 1),
  precision_scale_prior = exponential_prior(eta = 1),
  delta = NULL,
  edge_selection = TRUE,
  edge_prior = bernoulli_prior(0.5),
  precision_graph_prior = c("hierarchical", "joint"),
  na_action = c("listwise", "impute"),
  update_method = c("nuts", "adaptive-metropolis", "gibbs"),
  target_accept,
  nuts_max_depth = 10,
  learn_mass_matrix = TRUE,
  chains = 4,
  cores = parallel::detectCores(),
  display_progress = c("per-chain", "total", "none"),
  seed = NULL,
  verbose = getOption("bgms.verbose", TRUE),
  progress_callback = NULL
)

Arguments

Data

Argument	Description
x	A data frame or matrix with n rows and p columns. Columns may contain binary, ordinal, or continuous variables (see variable_type). Discrete variables are automatically recoded to non-negative integers (0, 1, ..., m); for regular ordinal variables, unobserved categories are collapsed, while Blume–Capel variables retain all categories. Continuous variables are column-centered internally so that the GGM likelihood is formulated with a zero-mean assumption.
variable_type	Character or character vector. Specifies the type of each variable in x. Allowed values: "ordinal", "blume-capel", or "continuous". A single string applies to all variables. A per-variable vector that mixes discrete ("ordinal" / "blume-capel") and "continuous" types fits a mixed MRF. Binary variables are automatically treated as "ordinal". Default: "ordinal".
baseline_category	Integer or vector. Baseline category used in Blume–Capel variables. Can be a single integer (applied to all) or a vector of length p. Required if at least one variable is of type "blume-capel".
na_action	Character. Specifies missing data handling. Either "listwise" (drop rows with missing values) or "impute" (impute within Gibbs sampler, propagating uncertainty). Default: "listwise". See Missing Data for details.

Prior specification

Priors are specified by passing prior-object constructors (see Prior Constructors). The legacy scalar arguments listed below the table (pairwise_scale, main_alpha, main_beta, inclusion_probability, beta_bernoulli_alpha/beta, …) are deprecated as of bgms 0.2.0.0; they still work but emit a lifecycle warning and forward to the equivalent prior object.

Argument	Description
interaction_prior	A bgms_parameter_prior object for pairwise interaction parameters. Allowed: normal_prior() (default), cauchy_prior(), beta_prime_prior(). Default: normal_prior(scale = 1).
threshold_prior	A bgms_parameter_prior object for threshold (main effect) parameters. Allowed: beta_prime_prior() (default), cauchy_prior(), normal_prior(). Default: beta_prime_prior(alpha = 0.5, beta = 0.5).
means_prior	A bgms_parameter_prior object for continuous-variable means (mixed MRF only). Allowed: normal_prior() (default), cauchy_prior(), beta_prime_prior(). Ignored for pure ordinal or pure GGM models. Default: normal_prior(scale = 1).
precision_scale_prior	A bgms_scale_prior object for the diagonal elements of the precision matrix (GGM and mixed MRF). Allowed: exponential_prior() (default), gamma_prior(). Both constructors accept the rate in one of two frames: rate (raw) or eta (standardized). eta is the rate on the diagonal in the coordinate system where the pairwise slab prior has unit scale; the raw rate is derived at fit time as eta / s, with s the scale of the interaction_prior. At fixed eta, graph and partial-correlation inference is invariant to the slab scale. The eta frame requires an interaction_prior with a scale parameter (normal_prior() or cauchy_prior()). Ignored for pure ordinal models. Default: exponential_prior(eta = 1).
delta	Non-negative numeric, or NULL for the dimension-adaptive default. Determinant-tilt exponent on the continuous-block precision matrix \(K\) (GGM) or \(K_{yy}\) (mixed MRF): multiplies the prior by \(|K|^{\delta}\), softly repelling the chain from the boundary of the positive-definite cone. NULL (default) auto-resolves to \(0.5 \log(p)\), with \(p\) the number of continuous variables. Pass delta = 0 for the untilted prior. Ignored for pure ordinal models.
edge_selection	Logical. Whether to perform Bayesian edge selection. If FALSE, the model estimates all edges. Default: TRUE.
edge_prior	A bgms_indicator_prior object specifying the prior for edge inclusion. Allowed: bernoulli_prior() (default), beta_bernoulli_prior(), sbm_prior(). Legacy character strings ("Bernoulli", "Beta-Bernoulli", "Stochastic-Block") are accepted but deprecated. For continuous and mixed models with beta_bernoulli_prior() or sbm_prior(), the hyperparameter updates carry a normalizing-constant correction; see the note below. Default: bernoulli_prior(0.5).
precision_graph_prior	Character. How the precision prior composes with the edge prior under edge selection for continuous (GGM) data. "hierarchical" (default): the hierarchical specification \(p(\Gamma)\,p(K \mid \Gamma)\) with the precision prior normalized per graph, so that the graph marginal is exactly the edge prior. Each edge move then evaluates a normalizer ratio with a fast local approximation, built once at the start of the analysis; a trust gauge runs in-chain and audits that approximation, and its summary is returned as fit$zratio_diag (see summarize_zratio_gauge()). Set options(bgms.zratio_gauge_sweeps = 0L) before fitting to switch the gauge off. "joint": the un-normalized joint specification, whose graph marginal is the edge prior times a graph-dependent normalizing constant; with beta_bernoulli_prior() or sbm_prior() the hyperparameter updates then carry the normalizing-constant correction described in the note below. The hierarchical specification requires edge_selection = TRUE, a normal_prior() or cauchy_prior() interaction prior, and continuous data — all-continuous, or mixed with at least two continuous variables (the normalizer then lives on the continuous block, and only continuous-continuous edge moves are affected). When the choice has no effect (edge_selection = FALSE, or no continuous precision block) the argument is accepted rather than an error: the surface and the gauge are skipped, and the posterior is identical under either value. Setting the argument yourself in that case prints a note when verbose = TRUE. Default: "hierarchical".

NoteWhat a continuous fit builds before the chains start

Both specifications precompute something once per analysis, and both cache it on disk via tools::R_user_dir("bgms", "cache"), so a repeat fit of the same configuration skips the build.

Under the default precision_graph_prior = "hierarchical" it is the normalizer surface the edge moves evaluate their ratio against. It announces itself with the size it covers and the number of cores it uses, then reports its elapsed time — this is the pause before the chains start on a first hierarchical fit. At the full 80-variable range it costs about 24 seconds serial and about 7 on four cores; a fit on fewer variables sizes the build to its own variable count. A hierarchical fit needs no edge-prior correction table, whatever its edge prior.

Under precision_graph_prior = "joint" with beta_bernoulli_prior() or sbm_prior() it is instead the edge-prior correction table, which the hyperparameter updates read: under the determinant-tilted precision prior, the plain conjugate updates would target the wrong marginals for the inclusion probability and the block structure. In mixed models the tilt acts on the continuous precision block, so the correction reads continuous-continuous edges only; with fewer than two continuous variables the plain conjugate updates apply unchanged. The table is keyed on the interaction prior and takes of the order of minutes to build. It announces itself only when it actually builds; a cache hit is silent. With beta_bernoulli_prior() the sampled inclusion probability is returned per chain in fit$inclusion_parameter_samples.

NoteDeprecated arguments

The following arguments are deprecated as of bgms 0.2.0. They still work but emit a lifecycle warning. Use the prior-object equivalents in the table above.

Deprecated	Replacement
standardize	No replacement; see the note below the table
pairwise_scale = s	interaction_prior = cauchy_prior(scale = s)
main_alpha = a, main_beta = b	threshold_prior = beta_prime_prior(alpha = a, beta = b)
inclusion_probability = p	edge_prior = bernoulli_prior(p)
beta_bernoulli_alpha = a, beta_bernoulli_beta = b	edge_prior = beta_bernoulli_prior(a, b)
beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda	edge_prior = sbm_prior(alpha_between = …, beta_between = …, dirichlet_alpha = …, lambda = …)
interaction_scale = s	interaction_prior = cauchy_prior(scale = s) (deprecated since 0.1.6.0)
threshold_alpha = a, threshold_beta = b	threshold_prior = beta_prime_prior(alpha = a, beta = b) (deprecated since 0.1.6.0)
burnin = n	warmup = n (deprecated since 0.1.6.0)
save	No replacement; raw samples are always stored in raw_samples (deprecated since 0.1.6.0)

standardize scaled each pair’s interaction prior by the product of the two variables’ maximum scores; pairwise interactions are now on the association scale and share one prior scale, so the adjustment is gone. standardize = FALSE, the old default, is what the sampler does, so it warns and fits; standardize = TRUE errors and points you at setting the scale directly through interaction_prior.

Sampler settings

Argument	Description
iter	Integer. Number of post-warmup iterations (per chain). Default: 2000.
warmup	Integer. Number of warmup iterations before collecting samples. Default: 2000. Short warmups trigger progressive warnings (NUTS only); see Warmup Schedule for the adaptation stages. With update_method = "gibbs" and edge selection, the first 15% of the warmup runs the full model (all edges included) so the precision matrix settles, and edge selection is active for the remaining 85%; both windows scale with the warmup budget.
update_method	Character. Specifies how the MCMC sampler updates the model parameters. "nuts" (default): the No-U-Turn Sampler, a gradient-based sampler available for all variable types, including under edge selection; continuous precision blocks use a free-element Cholesky parameterization that keeps the precision matrix positive-definite and encodes excluded edges through per-column null-space coordinates. "adaptive-metropolis": componentwise adaptive Metropolis–Hastings with Robbins–Monro proposal adaptation. "gibbs": a Gibbs sampler for the Gaussian graphical model, with a conjugate row-block draw of the precision matrix; available only for all-continuous data with a Normal or Cauchy interaction prior and a Gamma scale prior on the precision diagonal. Edge selection is supported for both slabs; the graph is updated by a full-conditional birth/death between-step (a Cauchy slab uses its scale-mixture representation).
target_accept	Numeric between 0 and 1. Target acceptance rate for the sampler. Defaults are set automatically if not supplied: 0.44 for adaptive Metropolis and 0.80 for NUTS. Under update_method = "nuts" this governs the step-size adaptation; the edge-move Metropolis proposals then use a fixed 0.44 target regardless of this argument. Ignored for "gibbs", which draws from exact conditionals and has no acceptance target.
nuts_max_depth	Integer. Maximum tree depth in NUTS. Must be positive. Default: 10.
learn_mass_matrix	Logical. If TRUE, adapt a diagonal mass matrix during warmup (NUTS only). If FALSE, use the identity matrix. Default: TRUE.
chains	Integer. Number of parallel chains to run. Default: 4.
cores	Integer. Number of CPU cores for parallel execution. Sampling uses min(cores, chains); some computations outside of sampling (such as building the edge-selection prior correction table) use all cores. Default: parallel::detectCores().
display_progress	Character. Controls progress reporting during sampling. Options: "per-chain" (separate bar per chain), "total" (single combined bar), or "none" (no progress). Default: "per-chain".
seed	Optional integer. Random seed for reproducibility. Must be a single non-negative integer. On Windows with RcppParallel >= 6.0.0, a fixed seed reproduces a run only at the same cores setting: under the oneTBB 2022 scheduler, runs at different core counts are statistically equivalent but not bit-for-bit identical.
verbose	Logical. If TRUE, prints informational messages during data processing (e.g., missing data handling, variable recoding). Defaults to getOption("bgms.verbose", TRUE). Set options(bgms.verbose = FALSE) to suppress messages globally.
progress_callback	Optional R function with signature function(completed, total), called at regular intervals during sampling with the number of completed and total iterations across all chains. Useful for external front-ends (e.g., JASP) that supply their own progress reporting. Default: NULL (no callback).

Value

An S7 object of class bgms with posterior summaries, posterior mean matrices, and raw MCMC draws. Supports print(), summary(), coef(), predict(), and simulate().

Main components:

• posterior_summary_main — Data frame with posterior summaries (mean, sd, MCSE, ESS, Rhat) for main-effect parameters. For OMRF models these are category thresholds; for mixed MRF models these are discrete thresholds and continuous means. NULL for GGM models.
• posterior_summary_quadratic — Data frame with posterior summaries for the residual variance parameters (GGM and mixed MRF). NULL for OMRF models.
• posterior_summary_pairwise — Data frame with posterior summaries for partial association parameters.
• posterior_summary_indicator — Data frame with posterior summaries for edge inclusion indicators (if edge_selection = TRUE).
• posterior_mean_main — Posterior mean of main-effect parameters. NULL for GGM models. For OMRF: a matrix (p x max_categories) of category thresholds. For mixed MRF: a list with $discrete (threshold matrix) and $continuous (q x 1 matrix of means).
• posterior_mean_pairwise — Symmetric matrix of posterior mean partial associations (zero diagonal). For continuous variables these are unstandardized partial correlations; for discrete variables these are half the log adjacent-category odds ratio. Use extract_precision(), extract_partial_correlations(), or extract_log_odds() to convert to interpretable scales.
• posterior_mean_residual_variance — Named numeric vector of posterior mean residual variances \(1/\Theta_{ii}\). Present for GGM and mixed MRF models; NULL for OMRF.
• posterior_mean_indicator — Symmetric matrix of posterior mean inclusion probabilities (if edge selection was enabled).
• raw_samples — A list of raw MCMC draws per chain with the following sublists:
  • main — List of main effect samples.
  • pairwise — List of pairwise effect samples.
  • indicator — List of indicator samples (if edge selection enabled).
  • rb_inclusion, rb_counts — Rao-Blackwellized inclusion draws (aligned with the indicator samples) and the per-edge odds accumulators behind them (if edge selection enabled); see extract_posterior_inclusion_probabilities() and extract_inclusion_bf().
  • allocations — List of cluster allocations (if SBM prior used).
  • nchains — Number of chains.
  • niter — Number of post-warmup iterations per chain.
  • parameter_names — Named lists of parameter labels.
• arguments — A list of function call arguments and metadata.

Additional summaries when edge_prior = sbm_prior(...) (see Sekulovski et al., 2025):

• posterior_summary_pairwise_allocations — Data frame with posterior summaries (mean, sd, MCSE, ESS, Rhat) for the pairwise cluster co-occurrence of the nodes. This serves to indicate whether the estimated posterior allocations, co-clustering matrix and posterior cluster probabilities have converged.
• posterior_mean_coclustering_matrix — Symmetric matrix of pairwise proportions of occurrence of every variable. This matrix can be plotted to visually inspect the estimated number of clusters and visually inspect nodes that tend to switch clusters.
• posterior_mean_allocations — A vector with the posterior mean of the cluster allocations of the nodes.
• posterior_mode_allocations — A vector with the posterior mode of the cluster allocations of the nodes.
• posterior_num_blocks — A data frame with the estimated posterior probabilities for all the possible numbers of clusters, computed under the shifted-Poisson prior on the number of clusters (see Edge Clustering).

Fields tied to specific settings:

• inclusion_parameter_samples — Per-chain draws of the sampled inclusion probability, present for continuous and mixed models with edge_prior = beta_bernoulli_prior(...).
• zratio_diag — Per-chain trust-gauge summary, present on any fit that runs the deployed hierarchical path: a continuous or mixed model with edge selection at the default precision_graph_prior = "hierarchical". It is NULL under "joint", where there is no approximation to audit, on models with no continuous precision block, and when the gauge is switched off with options(bgms.zratio_gauge_sweeps = 0L). See summarize_zratio_gauge().
• nuts_diag — NUTS diagnostics (tree depth, divergences, acceptance probabilities, energy, E-BFMI, and a warmup-completeness check), present when update_method = "nuts". Its summary includes mean_accept_prob, warmup_incomplete, and an overall has_issues flag that triggers a printed warning after fitting.

Details

For full explanations, see the guide. Model-specific assumptions and parameterization are described in Ordinal MRF, Gaussian Graphical Model, and Mixed MRF. Prior specification and structure learning are discussed in Prior Basics, Edge Selection, and Edge Clustering. To interpret outputs and verify chain quality, see Model Output and MCMC Diagnostics.

Examples

Fit an ordinal MRF to the Wenchuan PTSD data (the model below was run in advance; fitting takes a few minutes):

library(bgms)
fit = bgm(Wenchuan, iter = 1e4, warmup = 5e3, seed = 123)

Only the run length departs from the defaults; the longer run resolves the inclusion probabilities precisely enough to read as evidence (see From estimation to evidence).

summary(fit)

Posterior summaries from Bayesian estimation:

Category thresholds: 
                 mean  mcse    sd     n_eff  Rhat
intrusion (2)   0.445 0.002 0.249 12914.547 1.001
intrusion (3)  -2.075 0.005 0.389  6574.586 1.001
intrusion (4)  -5.356 0.009 0.657  5515.388 1.001
intrusion (5) -10.496 0.014 1.040  5784.383 1.001
dreams (2)     -0.817 0.002 0.212 10813.006 1.000
dreams (3)     -4.386 0.005 0.416  6887.781 1.000
... (use `summary(fit)$main` to see full output)

Pairwise interactions:
                     mean  mcse    sd     n_eff share_incl  Rhat
intrusion-dreams    0.365 0.000 0.036 22742.840      0.000 1.000
intrusion-flash     0.173 0.000 0.035 16205.811      0.000 1.000
intrusion-upset     0.062 0.001 0.052  1751.586      0.928 1.005
intrusion-physior   0.014 0.001 0.031  3481.262      0.954 1.001
intrusion-avoidth  -0.001 0.000 0.006 45609.250      0.162 1.000
intrusion-avoidact  0.000 0.000 0.005 38820.463      0.011 1.000
... (use `summary(fit)$pairwise` to see full output)

Inclusion probabilities:
                    mean  mcse    sd     n_eff  Rhat n0->1 n1->0
intrusion-dreams   1.000       0.000                     0     0
intrusion-flash    1.000     0 0.001 39460.069     1     0     0
intrusion-upset    0.664 0.013 0.460  1300.913 1.005   659   660
intrusion-physior  0.217 0.008 0.388  2477.845 1.001  1090  1090
intrusion-avoidth  0.034 0.001 0.109 21794.876     1  1018  1018
intrusion-avoidact 0.029 0.001 0.094 23071.446     1   913   913
... (use `summary(fit)$indicator` to see full output)
Note: NA values are suppressed in the print table; they occur for indicators
that were not updated or whose draws are constant, so ESS/Rhat are undefined.
`summary(fit)$indicator` still contains all computed values.

Use `summary(fit)$<component>` to access full results.
Use `extract_log_odds(fit)` for log odds ratios.
See the `easybgm` package for other summary and plotting tools.

The threshold and pairwise tables report, in that column order, the posterior mean, posterior standard deviation, Monte Carlo standard error (mcse), effective sample size, and R-hat for each parameter; the pairwise table adds share_incl, the share of that Monte Carlo error owed to inclusion rather than slab uncertainty. In the pairwise table, an estimate of 0.000 (e.g. intrusion-avoidact) is an edge that edge selection excluded in essentially every draw. The inclusion table (further down the printed output) gives each edge’s posterior inclusion probability in the mean column. See Model Output for a walkthrough of every component and MCMC Diagnostics for how to act on n_eff and Rhat.

See also

Getting Started, Methods, Extractor Functions, and bgmCompare().

References

Marsman, M., van den Bergh, D., & Haslbeck, J. M. B. (2025). Bayesian analysis of the ordinal Markov random field. Psychometrika, 90(1), 146–182. https://doi.org/10.1017/psy.2024.4

Sekulovski, N., Arena, G., Haslbeck, J. M. B., Huth, K. B. S., Friel, N., & Marsman, M. (2025). A stochastic block prior for clustering in graphical models. Manuscript Submitted for Publication. https://osf.io/preprints/psyarxiv/29p3m_v3