BGGM

Fit Bayesian Gaussian graphical models. The methods are separated into two Bayesian approaches for inference: hypothesis testing and estimation. There are extensions for confirmatory hypothesis testing, comparing Gaussian graphical models, and node wise predictability. These methods were recently introduced in the Gaussian graphical model literature, including Williams (2019) <doi:10.31234/osf.io/x8dpr>, Williams and Mulder (2019) <doi:10.31234/osf.io/ypxd8>, Williams, Rast, Pericchi, and Mulder (2019) <doi:10.31234/osf.io/yt386>.

Releases

13

Needs compilation

Yes

Needs R

>= 4.0.0

Config/roxygen2/version

8.0.0

Config/testthat/edition

3

Date

2026-09-23

LazyData

true

RdMacros

Rdpack

RoxygenNote

8.0.0

Earlier: Psychometric Models and Methods (until 7 Jul 2021).

• Graphical model estimation

Since March 2026

• BFpack
• GGally
• ggplot2
• ggridges
• grDevices
• MASS
• methods
• mvnfast
• network
• Rcpp
• RcppArmadillo
• RcppDist
• RcppProgress
• Rdpack
• reshape
• rlang
• sna
• stats
• truncnorm
• utils

Installing it adds

133 packages, 105.8 MB of source

Of those, need compilation

71

• abind
• assortnet
• knitr
• mice
• networktools
• psych
• rmarkdown
• testthat

Lines by kind (28,958)

Lines by language

Languages are counted from file extensions.

Files

231

R functions

205

Compiled functions

94

Object systems

S3: 31 methods

Median

1

Mean

3.93

Highest

43

Median lines per function

21

R functions

205

Calls

154

Unconnected

79

Groups

7

Has tests

No

Test frameworks

testthat

Test lines per line of code

0

• Controlling for Variables
• Custom Network Comparisons
• Custom Network Statistics
• Graphical VAR
• In Tandem: Confirmatory and Exploratory Testing
• MCMC Diagnostics
• Network Plots
• Predictability: Binary, Ordinal, and Continuous
• Testing Sums
• Three Ways to Test the Same Hypothesis

README

No

Release notes file

Yes

Help pages

84

Lines of documentation

6,259

Commits

1,212

Open pull requests

0

Tagged releases

0

Latest tagged release

–

Issues closed in

22 days (median)

Pull requests closed in

0 days (median)

Pull requests merged

100%

Created

1 Dec 2025

Licence file reads as

GPL-2.0

Method, repository practices

Across rast-lab

New features

• `select.explore()` gains a `method = "BMA"` option: Bayesian model averaging as an alternative to the hard Bayes factor threshold (method = "BF_cut"). For each edge the posterior is a spike-and-slab mixture with mass P(H0 | Y) at zero and the posterior under the alternative otherwise; the reported network weights (pcor_mat_zero) are the median of this mixture, computed exactly under a normal approximation of the posterior of the Fisher-z partial correlation (truncated at 0 for one-sided hypotheses). The result is therefore deterministic and does not require stored posterior draws. The new prior.prob.H0 argument (default 0.5) sets the prior probability of the null. All alternatives are supported. For alternative = "exhaustive" the mixture has three states -- a spike at zero (H0), a positive slab (H+), and a negative slab (H-) -- mixed by the posterior hypothesis probabilities; pos_mat/neg_mat/null_mat then classify each edge by the sign of the model-averaged median (so an edge is null when the median is 0, which can happen even when H0 is not the most probable hypothesis). Thanks to Joris Mulder for the suggestion.
• `select.explore()` returns edge inclusion probabilities (`incl_prob`): a matrix of posterior edge inclusion probabilities q * BF / (q * BF + 1 - q) with prior inclusion probability q = 1 - prior.prob.H0 (BF_10 for "two.sided", BF_20 for "greater"/"less", and 1 - P(H0 | Y) for "exhaustive").
• `prior.prob.H0` now affects the selected graph under `method = "BF_cut"`, not only under method = "BMA". An edge is selected when its posterior inclusion probability exceeds BF_cut / (BF_cut + 1), and is called null when the posterior probability of the null hypothesis exceeds that cutoff. With the default prior.prob.H0 = 0.5 this is identical to the previous Bayes factor rule (BF > BF_cut) for "two.sided", "greater" and "less". For alternative = "exhaustive" the null keeps prior.prob.H0 and the two directional hypotheses split the remainder equally (0.5 / 0.25 / 0.25 by default); because BF_1u + BF_2u = 2 the inclusion probability then equals the one of "two.sided", so both alternatives select the same edges, with "exhaustive" labelling each selected edge positive or negative by the larger of the two directional probabilities. alternative = "exhaustive" with method = "BF_cut" now also returns pcor_mat_zero. Previously "exhaustive" used fixed equal 1/3 priors and thresholded each hypothesis separately, which made BF_cut = 3 correspond to a two-sided Bayes factor of 1.5.
• `truncnorm` added to `Imports`: required for the truncated-normal draws used in the one-sided BMA alternatives.
• `explore()` no longer samples the prior by default: the Bayes factors only need the prior sd of the Fisher-z partial correlations, which is now computed analytically (numerical integration of the marginal prior rho ~ 2 * Beta(delta/2, delta/2) - 1, independent of p) and stored as prior_sd_z. Draws from the joint prior (prior_samp) are returned only with the new argument store_prior_draws = TRUE (default FALSE), which halves time and memory use for large networks by default. select.explore(), ggm_compare_explore() and bggm_missing() use prior_sd_z; older explore objects still work in select(). Bayes factors change slightly, because the previously sampled prior sd was affected by the matrix-F approximation (e.g. about 0.69-0.70 instead of 0.684 for prior_sd = 0.5).
• `explore()` gains `store_post_draws` (default `TRUE`): all samplers used by explore() now keep running posterior summaries of the partial correlations (post_samp$pcor_mat, pcor_sd, z_mean, z_sd). With store_post_draws = FALSE the p x p x iter arrays post_samp$pcors and post_samp$fisher_z are not stored, so memory use no longer grows with the number of iterations. select() (all methods and alternatives) and summary() work from these summaries when the draws are not stored; functions that need the draws (posterior_samples(), convergence(), coef(), pcor_to_cor(), predict(), predictability(), posterior_predict(), constrained_posterior()) stop with an informative error, and bggm_missing() always stores the draws.
• `explore()` gains `burnin` (default 50) and `thin` (default 1): after burnin iterations, iter post-burn-in iterations are run and every thin-th draw is stored; the running posterior summaries use all post-burn-in iterations. explore() objects no longer store the burn-in draws: post_samp$pcors, fisher_z, beta and thresh now contain only the kept draws (previously iter + 50, including 50 burn-in draws). Downstream functions locate the post-burn-in draws with an internal helper, so objects created with earlier versions keep working. Objects from estimate(), confirm() and the other functions are unchanged.

Bug fixes

• `explore(iter = )` is now the number of post-burn-in ITERATIONS, not the number of draws kept after thinning. Previously explore(iter = i, thin = t) ran burnin + i * t iterations and stored i draws, so raising thin silently multiplied the run time by t; it now runs burnin + i iterations and stores ceiling(i / thin) of them, returned in the new n_draws element. iter in the fitted object is the iteration count and n_draws the stored-draw count. With the default thin = 1 nothing changes.
• `seed = NULL` no longer re-seeds the RNG: explore(), confirm(), ggm_compare_confirm(), ggm_search() and var_estimate() called set.seed(seed) unconditionally before the guarded if (!is.null(seed)) set.seed(seed). With the default seed = NULL the first call was set.seed(NULL), which re-initialises the RNG from the current time and process ID, discarding any seed the user had set and making repeated calls on the same data irreproducible even inside set.seed(). The unguarded call has been removed, so with seed = NULL these functions now use the ambient RNG stream and are reproducible under set.seed(), as estimate() already was. Results of existing scripts that relied on the implicit re-seeding will change.
• Corrected the exhaustive posterior hypothesis probabilities: for alternative = "exhaustive", the positive/negative Bayes factors were referenced to the null model H0 instead of the unrestricted model Hu, so each carried an extra factor of the two-sided Bayes factor and double-counted the two-sided evidence. All three Bayes factors are now referenced to Hu per Eq. 9 of Williams & Mulder (2019). Affects post_prob and the derived null_mat/pos_mat/neg_mat under both method = "BF_cut" and method = "BMA".
• Corrected the prior sd in the `"greater"`, `"less"` and `"exhaustive"` selection branches of `select.explore()`: the "greater" branch averaged a hardcoded 3x3 prior mask (upper.tri(diag(3))), giving wrong Bayes factors whenever the number of variables was not 3, and the "less"/"exhaustive" branches averaged the (near-zero) matrix diagonal into the prior standard deviation. All branches now use the off-diagonal (edge) prior sd, matching the "two.sided" branch.
• Fixed a crash in `summary()` for `select.explore()` with `alternative = "less"`: the "less" summary branch produced an empty Relation column (mat_names[upper.tri(mat_names)] on an already-flattened vector), causing a "differing number of rows" error. The greater/less summaries now share one correct branch.
• `select.explore()` now uses all post-burn-in draws for the posterior mean/sd (previously the last 50 stored draws were dropped), and no longer returns NA on the diagonal of the exhaustive null_mat/pos_mat/neg_mat (where the posterior sd is 0).
• Fixed `bggm_missing()` dropping the wrong column with `mice` >= 3.17.0 (#2): bggm_missing() removed the .id column from mice::complete(action = "long") by position (column 2), which was correct only before mice 3.17.0. Since mice 3.17.0 places .imp/.id in the last two columns, this stripped a real data column and let .id leak into the model, corrupting the data passed to estimate()/ explore(). The column is now removed by name, which is robust to mice's column order. Reported by \@LilyTeesson.
• `bggm_missing()` now pools the posterior draws of the imputed data sets correctly: it stacked all draws of each fit, including their 50 burn-in draws, and set iter to iter * m + 50, so the methods that use draws 51:(iter + 50) included the burn-in draws of imputations 2 to m and indexed beyond the stored draws; summaries came from the first imputation only. The pooled object now contains the 50 burn-in draws of the first fit followed by the post-burn-in draws of all fits, with iter = iter * m. beta and thresh are pooled whenever present, and pcor_mat (and, for explore, the posterior summaries) are recomputed from the pooled draws. Also works for m = 1.
• `explore()` for large or n < p networks: the starting value is regularized (solve(cov(Y) + 0.1 I)), and the matrix-F prior now uses epsilon = min(0.01, 1 / (10 p)) instead of a fixed 0.01, so that nu = 1 / epsilon stays well above p - 1 (required for a proper prior, Williams & Mulder, 2020). Results change slightly for p > 10.
• C++ samplers no longer truncate the matrix-F degrees of freedom to integers (non-integer delta from prior_sd).
• Fixed `rref_ei` not found error: Added simple_rref() function in helpers.R to replace commented-out pracma::rref() call. This fixes a crash in create_matrices() when validating constraint matrices for hypothesis testing with multiple groups.
• Fixed `ggm_search()` crash: The C++ search function now handles edge cases where the adjacency matrix becomes all zeros or all ones, preventing "sample more elements than in x" errors.
• Fixed `qplot()` deprecation warning: Replaced deprecated qplot() with ggplot() + geom_density() in plot_prior().
• Fixed `ggm_search()` proposal-set shadowing bug and BIC creep: The C++ search() function had a variable-shadowing bug where zeros and nonzeros (pools of candidate edges) were re-declared instead of re-assigned after accepting a move, freezing them at the starting graph's configuration for the entire run. This caused systematic BIC drift under probabilistic acceptance and made reversing bad moves impossible. ggm_search() now implements a proper Metropolis-Hastings sampler with a birth-death Hastings correction as the default (probabilistic = TRUE); the greedy deterministic hill-climb is still available via probabilistic = FALSE. The burn_in parameter is restored and now actually applied to discard pre-burn-in samples before computing the Bayesian Model Averaging (BMA) solution (probabilistic search only). Note: the greedy hill-climb was found to be nearly non-functional on realistic test cases, typically accepting only a single edge flip out of thousands of attempts.

Performance

• Removed unused posterior arrays from the C++ samplers used by `explore()` (also used by estimate(), confirm() and the ggm_compare_*() functions) and by var_estimate(): Theta_mcmc, cors_mcmc and Sigma_mcmc (each p x p x iter) were allocated but never filled or returned in Theta_continuous, sample_prior, mv_continuous, mv_binary, mv_ordinal_albert, copula and var; copula also allocated an unused n x p x iter array of latent data. This reduces memory use substantially for large networks (e.g. mv_continuous, used with formula, allocated five p x p x iter arrays and now two).
• `missing_copula` (mixed data with `impute = TRUE`) no longer returns `post_samp$Y_collect`: this n x p x iter array was never filled (all zeros) and was not used anywhere in BGGM.

Maintenance

• Added rlang to Imports for proper use of .data pronoun in ggplot2 aesthetics.
• Added unit tests for multiple functions.
• Refactored `bggm_fast.cpp` for clarity and efficiency: Removed two functions unreachable from R (mv_ordinal_cowles, trunc_mvn; ~390 lines). Extracted cov_to_cor() and conditional_normal_params() helper functions, eliminating 17 duplicated blocks and a matrix inversion that was being computed twice per variable per MCMC iteration across 9 sampler functions. Removed dead commented-out debug code. No behavior change; file reduced from 3199 to ~2850 lines.
• Fixed the `ggm_search()` test suite and added regression tests: Corrected tests that were silently exercising non-existent parameters (bma = TRUE instead of bma_mean, a start = argument that doesn't exist, gamma = instead of prior_prob, prior_sd = which isn't a parameter at all) and a plot.ggm_search() test that could never fail since no such method exists. Rewrote with correct parameter names and assertions on the returned object's fields, and added regression tests tied directly to the bugs above (probabilistic search must accept a healthy fraction of proposals; bma_posterior() must not error on the new default's output). Added tests/testthat/test-search_cpp.R exercising the C++ search()/bic_fast()/hft_algorithm() functions directly, including a check for systematic BIC drift across a long run.