Large Gaussian Graphical Models
Gaussian graphical models often contain hundreds of variables, making edge selection computationally demanding. A graph with 200 variables has 19,900 candidate edges, and repeatedly updating the graph and its precision matrix can therefore become costly. For such analyses, bgm() provides a faster route, selected by two arguments:
fit = bgm(y,
variable_type = "continuous",
precision_graph_prior = "joint",
update_method = "gibbs"
)update_method = "gibbs" replaces the default No-U-Turn sampler with a conjugate sampler that redraws an entire column of the precision matrix at a time. precision_graph_prior = "joint" normalizes the prior over the precision matrix and graph once, globally, rather than separately for each graph. This puts every term in the edge move in closed form.
Neither option is the default, and the second changes the model being fitted. Under the joint specification, the graph prior is the edge prior reweighted toward sparsity. Its inclusion Bayes factors are therefore not interchangeable with those from the default specification. Gaussian Graphical Model describes both arguments and this trade-off.
This page compares the faster specification with two CRAN packages designed for the same task, Bayesian edge selection in Gaussian graphical models: ssgraph [version 1.16; Mohammadi (2025)], which samples the precision matrix under the continuous spike-and-slab prior of Wang (2015), and modelSelection [version 1.0.7; Rossell et al. (2026)], which performs Bayesian model selection over the entries of the precision matrix one column at a time, implementing the approach of Sulem et al. (2025).
Everything below uses simulated data: Erdős–Rényi graphs with average degree 3 and p = 50 to 200. Fits were run one at a time on the same machine (Apple M5 Pro, R 4.6.0, serial reference BLAS), with each package evaluated at its defaults and at settings recommended in its own documentation. The full protocol, the rationale for each setting quoted from the package help pages, and the complete tables are available in the repository.
Time to a converged ranking
Default run lengths are conventions, and here they differ by orders of magnitude: bgms runs 4 × 4000 full sweeps, ssgraph runs 5000, and modelSelection runs 1000 iterations, each visiting about sqrt(p) of the columns. Comparing default timings therefore compares conventions, not samplers. The relevant question is how much wall-clock time passes before the edge ranking stops changing. By that measure, the fast route earns its name.
We measured convergence on ten replicate datasets with p = 100 at each of three sample sizes: n = 1000, 500, and 200. We ran one chain on one thread per package, so chain parallelism does not confound the comparison.
Sample size controls the difficulty of the ranking problem. The generator’s per-edge partial correlations are 0.14 by construction, making n = 1000 an easy case and n = 500 and 200 substantially harder ones, as the converged AUCs below show. T is the time, including burn-in, required for the inclusion-probability AUC to come within 0.01 of its converged value. This is the convergence-cost metric of Vogels et al. (2024). Each grid point is a separate complete run using the package’s own burn-in policy. The table reports the median [min, max] over the ten replicates:
| n | bgms | ssgraph, g.prior 0.5 | modelSelection |
|---|---|---|---|
| 1000 | 0.50 s [0.45, 0.56] | 13.3 s [3.3, 13.5] | 1.02 s [1.00, 1.05] |
| 500 | 0.56 s [0.52, 0.62] | 65.4 s [13.6, 65.9] | 0.94 s [0.89, 0.98] |
| 200 | 0.65 s [0.55, 1.04] | 189 s [63, 254] | 0.53 s [0.44, 0.60] |
Each package is scored on the inclusion probabilities it reports to its user. bgms reaches its converged ranking in well under a second at every sample size. Relative to ssgraph with the matched 0.5 edge prior, the per-replicate median advantage is 25× at n = 1000 [range 6, 30] and grows to 235× at n = 200 [range 61, 424]; ssgraph is slower in all 30 replicate comparisons. The direction is consistent, but the magnitude at any one sample size is reliable only to within a factor of two or three, so we report medians with ranges rather than point estimates.
One thing these times are not: a recommended run length. The Rao-Blackwellized probabilities that bgm() reports average a conditional inclusion probability, so they stabilize before the indicator chain has fully explored — at n = 200 the median convergence point is 8 sweeps, where nearly half the pairs have not yet changed state. The times measure how quickly the ranking settles, not how long a chain should run. Scored on raw inclusion frequencies instead, which cannot converge until indicators move, bgms takes 0.5 to 7 seconds and still leads ssgraph by a median 18 to 25 times; those tables are in the repository.
The speed difference has a mechanical source. ssgraph’s spike is continuous, so an absent edge retains a small non-zero precision entry, and every column update draws a dense (p − 1)-variate conditional however sparse the sampled graph is. bgms uses a point-mass spike: an absent edge is exactly zero, and a column update works only through the edges included in the current graph. A point-mass sampler therefore benefits computationally from the sparsity it infers — telling bgms to expect a sparse graph (a 2/(p−1) edge prior instead of the Bernoulli(0.5) default) cuts a complete p = 200 fit from 394 to 122 seconds, while the same change moves ssgraph’s runtime by three percent.
modelSelection converges within about a second everywhere, and that speed is real. It reached the floor of the tested grid at every sample size in all ten replicates. Its reported time is therefore an upper bound: the short default run is converged for ranking, not undersampled. Its trade-offs are a different posterior (a 1/p edge prior, on the correlation scale), somewhat lower recall at a matched threshold (152 true edges against 168 for bgms under a comparable sparse prior), and no exact reproducibility: it exposes no seed, and its counts vary by about two percent across identical calls. bgms and ssgraph runs are reproducible bit for bit.
The converged rankings are effectively tied. bgms and ssgraph have identical median AUCs at every sample size: 1.00, 0.98, and 0.87 as the signal weakens. modelSelection is slightly lower at the two harder sample sizes, with median AUCs of 0.97 and 0.84, but the ranges overlap. Across the full benchmark, the inclusion-probability AUC exceeds 0.99 in every setting but one. The packages agree on the final ranking; where they differ is how quickly they reach it and where their default edge priors — 0.5 for bgms, 0.2 for ssgraph, 1/p for modelSelection — cut that ranking into selected and discarded edges. The prior, rather than the sampler, is what moves the false-positive count.
Several limits apply. The study fixes p = 100 and the average degree at 3. Because the advantage of the point-mass specification depends on sparsity, denser graphs may reduce it. The generator also assigns uniform edge strengths by construction, so the ranking problem lacks the heavy tail of genuinely weak edges that a G-Wishart draw would produce. Convergence of the AUC does not imply convergence of inclusion Bayes factors or tail quantities.
Finally, modelSelection is not part of the prior-matched comparison. Its posterior uses a 1/p edge prior and operates on the correlation scale, so it is not matched to the posterior used by bgms and ssgraph.
Rebuilding the tables
Every number on this page is traceable to a seeded runner script. The scripts, the rationale for each package setting, and the complete tables at all three sizes — including complete-run wall clocks and edge recovery at each package’s default and documented settings — are available in benchmarks/ggm-comparison/ in this site’s repository. The comparison can therefore be rerun whenever one of the packages releases a new version. Timings are facts about a machine; the ratios will travel further than the seconds.