library(bgms)
fit = bgm(Wenchuan, iter = 1e4, warmup = 5e3, seed = 123)Edge Selection
Bayesian edge selection determines which edges are present in the network by placing spike-and-slab priors on the partial association parameters. When edge_selection = TRUE (the default), bgm() estimates both the edge structure and the effect sizes simultaneously.
Spike-and-slab priors
For each edge \((i, j)\), an indicator variable \(\gamma_{ij} \in \{0, 1\}\) governs whether the interaction is included in the model:
- When \(\gamma_{ij} = 0\) (spike), the partial association \(\omega_{ij}\) is set to exactly zero — the edge is absent.
- When \(\gamma_{ij} = 1\) (slab), the partial association receives the slab prior set via
interaction_prior(defaultnormal_prior(scale = 1)) — the edge is present and the effect size is estimated.
In this formulation of the spike-and-slab prior (Gottardo & Raftery, 2008; van den Bergh et al., 2026), each edge is either exactly zero or follows a continuous distribution — with no overlap between the two components. An excluded edge is genuinely absent from the model, not merely shrunk toward zero, which is what makes exact statements like “the data support excluding this edge” possible.
The prior probability of inclusion is controlled by the edge_prior argument. See Prior Basics for the available options (Bernoulli, Beta-Bernoulli, Stochastic-Block).
Posterior inclusion probabilities
The posterior inclusion probability for edge \((i, j)\) is the posterior mean of \(\gamma_{ij}\). By default bgms estimates it by averaging the Rao-Blackwellized inclusion draws rather than the 0/1 indicators themselves, which targets the same quantity with less Monte Carlo noise (see Model Output). It ranges from 0 to 1:
- Values near 1.0: strong evidence that the edge is present.
- Values near 0.0: strong evidence that the edge is absent.
- Values near 0.5: the data are uninformative — there is not enough evidence to decide whether the edge is present or absent.
These are reported by coef(fit)$indicator and in the summary() output.
Bayes factors for edges
When the prior inclusion probability is \(\frac{1}{2}\), the posterior inclusion probability can be directly transformed into a Bayes factor (Kass & Raftery, 1995):
\[ \text{BF}_{10} = \frac{p(\gamma_{ij} = 1 \mid \mathbf{x})}{1 - p(\gamma_{ij} = 1 \mid \mathbf{x})} \]
For other prior inclusion probabilities, the Bayes factor divides the posterior inclusion odds by the prior inclusion odds — see The Bayesian Approach for the general identity.
In practice, the recommended tool is extract_inclusion_bf(), which returns the inclusion Bayes factor for every edge with the prior odds already divided out. It is computed from Rao-Blackwellized odds accumulators recorded during sampling, so it remains finite even for edges whose posterior inclusion probability saturates at 0 or 1 — exactly the edges where the transformation above overflows to \(\infty\) or collapses to 0. Pass log = TRUE to get the natural logarithm instead, which stays finite for evidence too extreme to represent as a Bayes factor.
The reciprocal \(1 / \text{BF}_{10}\) quantifies the evidence for edge absence, i.e., conditional independence (Sekulovski et al., 2024). Values of \(\text{BF}_{10}\) close to one indicate that the data do not discriminate between the two hypotheses — absence of evidence rather than evidence of absence.
When posterior odds equal the Bayes factor
The shortcut “posterior odds = Bayes factor” holds when the prior inclusion odds equal one, and that takes two ingredients. First, the edge prior must put probability \(\frac{1}{2}\) on inclusion. Under bernoulli_prior(0.5) it does; under, say, bernoulli_prior(0.1) the prior odds are \(1/9\), and the posterior odds understate the evidence for inclusion by a factor of nine, whatever the data type. Second, the model must use the probability you wrote down. For discrete (ordinal) data it always does. For continuous (GGM) data — and for the continuous-continuous edges of a mixed model — it does under the default precision_graph_prior = "hierarchical", which normalizes the precision prior per graph so that the prior over structures is exactly the edge prior (see Prior Basics).
The second ingredient fails if you opt into precision_graph_prior = "joint". The edge prior there combines with the un-normalized precision-matrix prior, and the effective prior inclusion probability of a continuous edge differs from the nominal edge-prior value. On a joint fit, compute the Bayes factor from the general identity — posterior inclusion odds divided by prior inclusion odds — with the prior probabilities from extract_prior_inclusion_probabilities():
post = extract_posterior_inclusion_probabilities(fit)
prior = extract_prior_inclusion_probabilities(fit)
BF10 = (post / (1 - post)) / (prior / (1 - prior))Both extractors return matrices of matching shape, so the odds combine element-wise.
In every case, extract_inclusion_bf() computes the Bayes factor from that same identity with the prior odds the fit actually used, so it needs no adjustment: not for an edge prior away from \(\frac{1}{2}\), and not for either precision_graph_prior specification.
Other approaches to testing conditional independence
The inclusion Bayes factor is not the only Bayesian approach to testing conditional independence. Two common alternatives can also be computed from bgms output when edge_selection = FALSE. Sekulovski et al. (2024) provide a detailed comparison of all three methods.
Credible intervals. With edge_selection = FALSE, posterior credible intervals for each partial association \(\omega_{ij}\) can be constructed from the MCMC samples (available in fit$raw_samples). If the interval excludes zero, this is sometimes taken as evidence against the absence of an edge.
Savage–Dickey Bayes factor. With edge_selection = FALSE, the posterior samples from the full model can be used to compute a Savage–Dickey density ratio for each edge, testing whether \(\omega_{ij} = 0\). This approach assumes that the rest of the network is fully connected.
As an alternative to these two, the inclusion Bayes factor is what bgms reports by default. When edge_selection = TRUE, the edge indicator samples are used to compute an inclusion Bayes factor for each edge, testing whether \(\omega_{ij} = 0\). Unlike the Savage–Dickey approach, this method does not assume that the rest of the network is fully connected, but weighs each possible structure by its posterior plausibility. It thus tests each edge while accounting for structural uncertainty. See The Bayesian Approach.
Evidence categories
The following classification, based on Kass & Raftery (1995), is a common convention for interpreting Bayes factors:
| Bayes factor | Evidence |
|---|---|
| 1–3 | Not worth more than a bare mention |
| 3–10 | Substantial |
| 10–100 | Strong |
| > 100 | Decisive |
The same scale applies for evidence of absence using \(1 / \text{BF}_{10}\). These thresholds are conventions, not absolute standards — the Bayes factor is a continuous measure of evidence and does not require strict cutoffs.
Evidence in practice
As an illustration from psychology, Huth et al. (2026) reanalyzed 293 networks from 126 studies and found that the network structure is highly uncertain: for the majority of edges, the data do not clearly decide whether the edge should be present or absent. Only about one in five edges was supported by strong evidence; for the rest, the evidence was weak or inconclusive. The exact proportions depend on the evidence thresholds used, but the overall picture is clear — uncertainty about individual edges is the norm, not the exception.
This is precisely why a Bayesian approach matters. Without it, a missing edge in an estimated network is ambiguous: did the method find evidence that the edge is absent, or was there simply not enough data to detect it? The Bayes factor makes this distinction explicit. When the evidence is inconclusive, it says so — and that information is just as valuable as a decisive result, because it prevents conclusions that the data cannot support.
A worked example
To see these quantities on real output, this snippet reads three edges off the Wenchuan fit that From estimation to evidence runs end to end (run in advance here):
That page owns the analysis, including why this fit runs longer than the default; here we only read the quantities this page has introduced.
The summary(fit)$indicator table holds the posterior inclusion probability of every edge in its mean column (the Rao-Blackwellized estimate; see Model Output). With the default prior inclusion probability of 0.5, the inclusion Bayes factor is the posterior odds:
ind = summary(fit)$indicator
edges = ind[c("intrusion-dreams", "intrusion-upset", "intrusion-avoidact"),
"mean", drop = FALSE]
edges$BF10 = edges$mean / (1 - edges$mean)
round(edges, 3) mean BF10
intrusion-dreams 1.000 Inf
intrusion-upset 0.664 1.976
intrusion-avoidact 0.029 0.030Three different conclusions, one from each row:
intrusion-dreams: the inclusion probability rounds to 1 — every sweep left this edge included — so the posterior-odds transformation overflows to infinity.intrusion-upset: an inclusion probability of about 0.66 gives \(\text{BF}_{10} \approx 2.0\) — the data lean toward including this edge, but the evidence is weak.intrusion-avoidact: an inclusion probability of about 0.03 gives \(\text{BF}_{01} \approx 33\) — strong evidence that this edge is absent (evidence of absence).
For the saturated edge, extract_inclusion_bf() recovers a finite answer from the odds accumulators:
bf = extract_inclusion_bf(fit)
signif(bf["intrusion", c("dreams", "upset", "avoidact")], 3) dreams upset avoidact
2.12e+119 1.98e+00 3.04e-02 intrusion-dreams gets a Bayes factor of about \(2 \times 10^{119}\): astronomically strong evidence, but now an actual number. The other two entries reproduce the posterior-odds calculation above (\(2.0\) for intrusion-upset, and \(0.030 \approx 1/33\) for intrusion-avoidact).
Evidence that strong runs into the limits of double precision: a Bayes factor beyond about \(1.8 \times 10^{308}\) comes back as Inf. Reading the same three edges on the log scale keeps every entry finite and comparable:
round(extract_inclusion_bf(fit, log = TRUE)["intrusion",
c("dreams", "upset", "avoidact")], 2) dreams upset avoidact
274.76 0.68 -3.49 Median probability model
A common summary of the posterior network is the median probability model (Barbieri & Berger, 2004): retain all edges with posterior inclusion probability above 0.5. This corresponds to including edges for which there is more evidence for presence than absence.
network = coef(fit)$pairwise
network[coef(fit)$indicator < 0.5] = 0Use this as a summary of what the posterior favors, not as a picture of what the data decided. Thresholding at 0.5 collapses the three-way distinction this page has just drawn: an edge kept at an inclusion probability of 0.51 is undecided, an edge dropped at 0.49 is undecided too, and neither blank in the resulting drawing means the data support conditional independence. When the network is to be displayed, draw the three-panel edge evidence plot instead, which keeps the undecided pairs visible as their own category; Getting Started shows how.
Prior guidance for edge selection
The default edge_prior = bernoulli_prior(0.5) is a standard noninformative choice. Considerations for changing it:
- Sparse networks: lower the inclusion probability (e.g.,
bernoulli_prior(0.25)) to encode a prior expectation that most edges are absent. - Dense networks: raise the inclusion probability if prior knowledge suggests many edges.
- Learned inclusion probability: use
edge_prior = beta_bernoulli_prior(1, 1)to let the data inform the overall sparsity level. - Structured sparsity: use
edge_prior = sbm_prior(...)when edges are expected to cluster (see Edge Clustering).
The interaction_prior controls the slab width. A wider slab makes it harder for the Bayes factor to favor inclusion when the true effect is small, because the prior spreads probability over a wide range. The default is normal_prior(scale = 1).
To see how much your edge verdicts depend on that width, run prior_sensitivity_check() on the fit. It recovers each edge’s inclusion Bayes factor across a range of slab scales around the one you used, and separates edges whose verdict genuinely changes from edges that merely wobble between reruns of the same prior.
References
Barbieri, M. M., & Berger, J. O. (2004). Optimal predictive model selection. The Annals of Statistics, 32(3), 870–897. https://doi.org/10.1214/009053604000000238
Gottardo, R., & Raftery, A. E. (2008). Markov chain Monte Carlo with mixtures of mutually singular distributions. Journal of Computational and Graphical Statistics, 17(4), 949–975. https://doi.org/10.1198/106186008X386102
Huth, K. B. S., Haslbeck, J. M. B., Keetelaar, S., van Holst, R. J., & Marsman, M. (2026). Statistical evidence in psychological networks. Nature Human Behaviour, 10, 333–346. https://doi.org/10.1038/s41562-025-02314-2
Kass, R. E., & Raftery, A. E. (1995). Bayes factors. Journal of the American Statistical Association, 90(430), 773–795. https://doi.org/10.2307/2291091
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
van den Bergh, D., Clyde, M. A., Raftery, A. E., & Marsman, M. (2026). Reversible jump MCMC with no regrets: Bayesian variable selection using mixtures of mutually singular distributions. Manuscript in Preparation.