# Getting Started with bgms

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.

# Comparing groups

`bgmCompare()` estimates group differences in the category thresholds and in the pairwise interactions, and puts those differences under selection, so each one carries an inclusion Bayes factor of its own.

    fit_compare = bgmCompare(
      x = ADHD[ADHD$group == 1, 2:6],
      y = ADHD[ADHD$group == 0, 2:6],
      seed = 1234
    )
    verdicts(fit_compare)
    plot(fit_compare)

`verdicts()` and `plot()` mean the same things here as above, with presence reading as “the groups differ on this pair” and absence as evidence that they do not – a conclusion a group comparison can otherwise rarely state. One caveat belongs on the front door: difference verdicts are priced by `difference_scale`, whose calibration under the association-scale parameterization is still under study, so a difference verdict close to a threshold should be read as scale-contingent. The Model Comparison vignette works through a full example.

# Where to go next

- Checking your fitted model – `verdicts()` in depth, `calibration_check()`, and building a posterior predictive display on `simulate()`.
- Diagnostics and Spike-and-Slab Summaries – R-hat, effective sample size, the NUTS diagnostics, and the hierarchical prior trust gauge.
- Checking Prior Sensitivity – whether the verdicts depend on the slab scale.
- Model Comparison with bgmCompare – group differences on real data.
- The easybgm package builds on bgms fits and offers a wider set of summary and plotting workflows around them.
- `NEWS.md` lists what changed in 0.2.0.0, including the change of scale for the pairwise effects and the arguments that were removed.

# References

Geng, J., Bhattacharya, A., & Pati, D. (2019). Probabilistic community detection with unknown number of communities. Journal of the American Statistical Association, 114, 893–905. https://doi.org/10.1080/01621459.2018.1458618
Huth, K., de Ron, J., Goudriaan, A. E., Luigjes, K., Mohammadi, R., van Holst, R. J., Wagenmakers, E.-J., & Marsman, M. (2023). Bayesian analysis of cross-sectional networks: A tutorial in R and JASP. Advances in Methods and Practices in Psychological Science. https://doi.org/10.1177/25152459231193334
Marsman, M., Bergh, D. van den, & Haslbeck, J. M. B. (2025). Bayesian analysis of the ordinal Markov random field. Psychometrika, 90, 146–182. https://doi.org/10.1017/psy.2024.4
McNally, R. J., Robinaugh, D. J., Wu, G. W. Y., Wang, L., Deserno, M. K., & Borsboom, D. (2015). Mental disorders as causal systems: A network approach to posttraumatic stress disorder. Clinical Psychological Science, 3, 836–849. https://doi.org/10.1177/2167702614553230
Sekulovski, N., Keetelaar, S., Huth, K. B. S., Wagenmakers, E.-J., van Bork, R., van den Bergh, D., & Marsman, M. (2024). Testing conditional independence in psychometric networks: An analysis of three Bayesian methods. Multivariate Behavioral Research, 59, 913–933. https://doi.org/10.1080/00273171.2024.2345915