bgms provides Bayesian analyses of graphical models, all of which are Markov random fields (MRFs). In an MRF, two variables are connected by an edge when they remain associated after accounting for all other variables in the model. The variable_type argument selects the model family: the ordinal Markov random field, with binary variables as a special case; the Blume–Capel model; the Gaussian graphical model for continuous variables; or the mixed Markov random field, which joins discrete and continuous variables in one model.

The graph is treated as unknown rather than fixed. The analysis returns a posterior distribution over graphs and, for every pair of variables, an inclusion Bayes factor: the factor by which the data update the odds that the edge is present. The reporting convention assigns each pair one of three verdicts: evidence of presence, evidence of absence, or undecided. The verdicts() function reports these classifications.

This vignette is the front door: enough to run an analysis and read what comes back, and no more. The inclusion Bayes factor and its three-state interpretation are developed in Huth et al. (2023) and Sekulovski et al. (2024). Extended teaching material, including model background, Markov Chain Monte Carlo (MCMC) output, and worked analyses, is available on the package website, https://bayesian-graphical-modelling-lab.github.io/bgms-docs/. The other vignettes develop the individual components in greater detail.

The models

There are two entry points. bgm() fits one Markov random field to one sample. bgmCompare() fits several groups at once and tests where their graphs differ; it takes binary and ordinal data.

Which member bgm() fits follows from variable_type, not from a separate argument.

"ordinal" (the default) fits the ordinal Markov random field: one threshold per category per variable, and one pairwise interaction per pair (Marsman et al., 2025). Binary variables are the two-category special case and are handled here.
"blume-capel" fits the Blume–Capel model, which replaces a variable’s free thresholds with a linear and a quadratic term around a baseline_category. It is the economical choice when a variable has many categories.
"continuous" fits the Gaussian graphical model: the pairwise parameters are read off the precision matrix, and an edge is a non-zero partial association.
A vector with one entry per variable mixes the members. A vector holding both discrete and "continuous" entries fits the mixed Markov random field, which carries discrete and continuous variables in one graph.

These are members of one family rather than separate modelling frameworks, and everything downstream – edge selection, the summaries, the verdicts, the plots – is the same whichever member is fitted.

A worked example

The Wenchuan dataset holds responses from survivors of the 2008 Wenchuan earthquake on 17 posttraumatic stress items (McNally et al., 2015). Nine items are enough to show the whole workflow.

library(bgms)
data = Wenchuan[, 1:9]
fit = bgm(data, seed = 1234)

That call is the whole model specification: bgm() reads the variable types, applies its default priors (below), turns edge selection on, and runs four chains of NUTS. Setting seed makes the run reproducible. The fit built for this vignette uses two chains rather than four to keep the build short; nothing else about it differs.

What the fit contains

summary() prints the posterior in three blocks, each with its own Monte Carlo diagnostics beside the estimates, so the numbers and the question of whether to trust them are never in separate places.

summary(fit)
#> Posterior summaries from Bayesian estimation:
#> 
#> Category thresholds: 
#>                 mean  mcse    sd    n_eff  Rhat
#> intrusion (2)  0.463 0.005 0.229 2434.348 1.000
#> intrusion (3) -1.857 0.009 0.334 1531.636 1.001
#> intrusion (4) -4.749 0.015 0.545 1345.727 1.000
#> intrusion (5) -9.341 0.022 0.869 1497.287 1.001
#> dreams (2)    -0.613 0.005 0.202 1788.399 1.001
#> dreams (3)    -3.864 0.010 0.385 1369.451 1.001
#> ... (use `summary(fit)$main` to see full output)
#> 
#> Pairwise interactions:
#>                     mean  mcse    sd    n_eff share_incl  Rhat
#> intrusion-dreams   0.332 0.001 0.034 3048.269      0.000 1.000
#> intrusion-flash    0.153 0.001 0.031 3185.816      0.000 1.000
#> intrusion-upset    0.128 0.002 0.037  262.216      0.035 1.008
#> intrusion-physior  0.038 0.003 0.045  202.706      0.965 1.013
#> intrusion-avoidth  0.000 0.000 0.006 2225.366      0.013 1.000
#> intrusion-avoidact 0.000 0.000 0.005 3770.473      0.132 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.000     4000     1     0     0
#> intrusion-upset    0.991 0.003 0.093  801.618 1.008     5     5
#> intrusion-physior  0.471 0.039 0.484  153.872 1.011    92    91
#> intrusion-avoidth  0.030 0.002 0.095 1737.468     1    87    87
#> intrusion-avoidact 0.030 0.002 0.103 2513.991 1.001    80    80
#> ... (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 first block holds the category thresholds, the second the pairwise interactions, and the third the posterior inclusion probability of each edge. The mcse, n_eff, and Rhat columns describe the sampling rather than the data: how precisely each quantity was estimated, and whether the chains agree. The diagnostics vignette explains how to read them, including why the usual Rhat < 1.01 rule does not transfer to the binary edge indicators.

A pairwise interaction is on the association scale: it is the coefficient omega that enters every conditional distribution as 2 * omega * x, which for an ordinal model makes it half the log odds ratio between adjacent response categories. This changed in 0.2.0.0 – earlier versions stored twice this number – so a value taken from a 0.1.6.3 run is not comparable without the factor of two. summary(fit)$pairwise and its siblings return each block in full, and coef(fit) returns the posterior means on their own.

Which edges the evidence settles

verdicts() reads each edge’s inclusion Bayes factor against a threshold and its reciprocal, and returns one of three answers.

verdicts(fit)
#> Edge verdicts at an inclusion Bayes factor of 10 (and 0.1 for absence):
#> presence: log BF > 2.30; absence: log BF < -2.30
#> 
#>   presence 11 | undecided 15 | absence 10   (36 indicators)
#> 
#>           parameter   pip  log_bf   verdict fragile
#>    intrusion-dreams 1.000 279.606  presence   FALSE
#>     intrusion-flash 1.000  17.175  presence   FALSE
#>     intrusion-upset 0.991   4.654  presence   FALSE
#>   intrusion-physior 0.471  -0.118 undecided   FALSE
#>   intrusion-avoidth 0.030  -3.469   absence   FALSE
#>  intrusion-avoidact 0.030  -3.463   absence   FALSE
#>   intrusion-amnesia 0.036  -3.283   absence   FALSE
#>   intrusion-lossint 0.201  -1.382 undecided   FALSE
#>        dreams-flash 1.000 138.442  presence   FALSE
#>        dreams-upset 0.516   0.066 undecided   FALSE
#> ... (26 more rows)
#> 
#> 6 verdicts are Monte-Carlo fragile: a verdict boundary lies within two standard
#> errors of the evidence, so the verdict could change on a rerun. Consider a
#> longer run.

The default threshold is 10, so an edge is called presence when the data multiply its inclusion odds by more than 10, absence when they divide them by more than 10, and undecided in between. That third category is not a failure of the analysis. It is the honest answer when the data separate neither hypothesis from the other, and being able to report it is the point of treating the graph as unknown.

The evidence is displayed as log_bf, the natural logarithm of the Bayes factor: zero is even odds, positive favours presence, negative favours absence, and the log scale stays finite where the Bayes factor itself would overflow. extract_inclusion_bf() returns the Bayes factors themselves, or their natural logarithms with log = TRUE.

The fragile column is the part that cannot be read off the Bayes factor. A verdict is fragile when a threshold sits within two standard errors of the estimated evidence, which is the regime in which Monte Carlo noise, rather than the data, decides the answer. A fragile verdict is not a wrong verdict; it is one the run was too short to settle, and the remedy is more iterations. The known-truth calibration study behind the flag is described in ?verdicts.

The edge evidence plot

plot() draws the three verdicts as three panels on one shared layout: the pairs the data support, the pairs the data rule out, and the pairs the data cannot decide. A single drawing of the graph would have to collapse the last two into the same blank.

plot(fit)

Each panel is titled with what it holds, how many pairs are in it, and the rule that put them there. Only the first panel is weighted – line width is the posterior mean association and colour carries its sign – because that is where the effect sizes are; the other two are drawn at uniform width, since for those pairs the classification is the result. The layout is computed once from all pairs, so a node sits in the same place in every panel. The threshold is the one verdicts() uses, so the figure and the table cannot say different things.

Three panels side by side want a wide device: open one at roughly width = 13, height = 5 before plotting, or the labels crowd. Drawing needs the qgraph package, which bgms suggests rather than depends on; without it plot() stops and points at verdicts() for the same information as a table. plot(fit, type = "centrality") draws posterior strength centrality instead, and plot_edge_posterior() opens a single pair, showing how its posterior mass divides between the edge being absent and the edge being present.

Priors

Edge selection rests on a spike-and-slab prior for each pairwise interaction: a point mass at zero for the edge being absent, and a slab for the values it can take when the edge is present. Priors are supplied as objects built by constructor functions rather than as loose numbers, so a prior always states its family as well as its scale.

fit = bgm(data,
  interaction_prior = cauchy_prior(scale = 2.5),
  edge_prior = beta_bernoulli_prior(alpha = 1, beta = 1)
)

The parameter priors are normal_prior(), cauchy_prior(), and beta_prime_prior(); the precision-diagonal priors are exponential_prior() and gamma_prior(); the edge priors are bernoulli_prior(), beta_bernoulli_prior(), and sbm_prior(). The last of these is a stochastic-block prior, under which the edges are free to cluster (Geng et al., 2019).

The defaults bgm() runs at are:

Argument	Default	Governs
`interaction_prior`	`normal_prior(scale = 1)`	pairwise interactions (the slab)
`threshold_prior`	`beta_prime_prior(0.5, 0.5)`	category thresholds
`means_prior`	`normal_prior(scale = 1)`	continuous means (mixed models)
`precision_scale_prior`	`exponential_prior(eta = 1)`	precision diagonal (continuous and mixed models)
`edge_prior`	`bernoulli_prior(0.5)`	edge inclusion

bgmCompare() prices its baseline pairwise interactions with the same normal_prior(scale = 1), so the two entry points no longer ship different priors under one name; its group differences are priced separately, by difference_family ("Normal" by default) and difference_scale.

Two consequences are worth stating plainly. The default slab is Normal, where 0.1.6.3 used a Cauchy at scale 2.5, and the two are different models; interaction_prior = cauchy_prior(scale = 2.5) restores the old one, on the new coordinate. And because the edge prior defaults to a fifty-fifty bernoulli_prior(0.5), the inclusion Bayes factor is exactly the posterior inclusion odds, which is why the pip and log_bf columns of verdicts() track each other so closely.

The old scalar arguments – pairwise_scale, inclusion_probability, beta_bernoulli_alpha and the rest – still work, but they warn and are translated into the corresponding object. Note that pairwise_scale = s means cauchy_prior(scale = s), preserving the family it had in 0.1.6.3, so a deprecated call is not the same model as a call at the new default.

The scale of the interaction prior is the assumption a reader is most likely to disagree with. prior_sensitivity_check() traces every edge’s inclusion Bayes factor across a range of slab scales and reports which verdicts depend on the choice; the Checking Prior Sensitivity vignette walks through its report.

The Gaussian and mixed members

A continuous fit is the same call with a different variable_type.

# 200 draws from a five-variable Gaussian graphical model whose precision
# matrix is a chain, 1-2-3-4-5.
set.seed(1234)
K = diag(5)
K[cbind(1:4, 2:5)] = -0.4
K[cbind(2:5, 1:4)] = -0.4
continuous_data = matrix(rnorm(200 * 5), 200, 5) %*% chol(solve(K))

fit_ggm = bgm(continuous_data, variable_type = "continuous", seed = 1234)
verdicts(fit_ggm)

The pairwise effects of a Gaussian graphical model are unstandardized partial associations, on the same association scale as everywhere else – half the off-diagonal precision entry, -0.5 * K_ij. They are not partial correlations. extract_partial_correlations() converts them, and extract_precision() returns the precision matrix itself. A mixed fit reports the same quantities for its continuous block.

One default is specific to continuous data. precision_graph_prior chooses how the prior on the precision matrix composes with the prior on the graph, and it defaults to "hierarchical", under which the graph marginal is exactly the edge prior that was asked for. Obtaining that requires a normalizing constant at every edge move, which bgms approximates rather than computes exactly; a trust gauge then audits the approximation against the exact calculation after sampling, stores the audit in fit$zratio_diag, and warns when the shortcut is changing edge decisions or leaning the inclusion probabilities. ?bgm documents the approximation and its measured accuracy, and the diagnostics vignette explains how to read the gauge and what to do about a flag.

Missing values are handled by listwise deletion by default; na_action = "impute" integrates over them during sampling instead, for all members of the family.

Getting Started with bgms