Reading the output

You have a fitted model, you call summary(), and the console fills with tables. This page walks through that printout block by block on a real fit: what each number is, which numbers are the result, which are housekeeping, and, for every block, how to tell whether its numbers are okay. Misreading one column is what this page protects you from: the posterior means in the pairwise table are model-averaged, so treating one as “the effect size, given that the edge exists” overstates edges the data have not decided on. Keep the page open next to your console; the closing table says where to find the numbers that summary() does not print.

The fit

library(bgms)

fit = bgm(Wenchuan, iter = 1e4, warmup = 5e3, seed = 123)

The Wenchuan data ship with the package: 17 ordinal PTSD symptom items, with 344 of the 362 survivors complete under the default na_action = "listwise". Everything except the run length is at its default; see bgm() for the arguments and From estimation to evidence for why this fit runs longer than the default.

The output in full

summary(fit)

Posterior summaries from Bayesian estimation:

Category thresholds: 
                 mean  mcse    sd     n_eff  Rhat
intrusion (2)   0.445 0.002 0.249 12914.547 1.001
intrusion (3)  -2.075 0.005 0.389  6574.586 1.001
intrusion (4)  -5.356 0.009 0.657  5515.388 1.001
intrusion (5) -10.496 0.014 1.040  5784.383 1.001
dreams (2)     -0.817 0.002 0.212 10813.006 1.000
dreams (3)     -4.386 0.005 0.416  6887.781 1.000
... (use `summary(fit)$main` to see full output)

Pairwise interactions:
                     mean  mcse    sd     n_eff share_incl  Rhat
intrusion-dreams    0.365 0.000 0.036 22742.840      0.000 1.000
intrusion-flash     0.173 0.000 0.035 16205.811      0.000 1.000
intrusion-upset     0.062 0.001 0.052  1751.586      0.928 1.005
intrusion-physior   0.014 0.001 0.031  3481.262      0.954 1.001
intrusion-avoidth  -0.001 0.000 0.006 45609.250      0.162 1.000
intrusion-avoidact  0.000 0.000 0.005 38820.463      0.011 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.001 39460.069     1     0     0
intrusion-upset    0.664 0.013 0.460  1300.913 1.005   659   660
intrusion-physior  0.217 0.008 0.388  2477.845 1.001  1090  1090
intrusion-avoidth  0.034 0.001 0.109 21794.876     1  1018  1018
intrusion-avoidact 0.029 0.001 0.094 23071.446     1   913   913
... (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.

For this model type, an ordinal Markov random field with edge selection, the printout has three tables. Each row is a parameter, each column a posterior summary or a sampler diagnostic. The console shows the first six rows of each table; the full tables live in summary(fit)$main, summary(fit)$pairwise, and summary(fit)$indicator. The blanks in the last table are suppressed NA values, read below. The sections that follow take the blocks in print order.

Category thresholds

Each variable’s thresholds control how likely its response categories are, apart from any associations with other variables. The lowest observed category is the reference and gets no row; each other category gets one row, labeled with its value on the original scale, so intrusion (2) is the threshold of category 2 on the intrusion item. The increasingly negative values down each variable’s rows say that higher symptom categories are increasingly rare. See Ordinal MRF for the model behind them.

This first table also introduces the five columns that every block repeats:

mean is the posterior mean, the parameter estimate.
sd is the posterior standard deviation: your uncertainty about the parameter. More data would shrink it; a longer run would not.
mcse is the Monte Carlo standard error: how precisely the finite chain has pinned down the posterior mean. A longer run shrinks it. See MCMC Diagnostics.
n_eff is the effective sample size: how many independent draws the chain is worth. As a rule of thumb, 100 buys one significant digit of precision, and values below 100 also make the diagnostics themselves unreliable (Vehtari et al., 2021). See MCMC Diagnostics.
Rhat compares the four chains with each other and each chain’s halves with each other. Values above 1.01 mean the chains disagree and the run should be longer (Vehtari et al., 2021). See MCMC Diagnostics.

Is this block okay? Here every threshold has n_eff in the thousands and Rhat at most 1.001, comfortably clear of both guidelines. The thresholds are nuisance parameters: they have to be in the model, but you will rarely report one.

Pairwise interactions

This block holds the result: one row per pair of variables, with mean the posterior mean partial association, the strength of the direct connection that remains after conditioning on the other 15 items.

The means are model-averaged. With edge selection on, sweeps that exclude an edge set its weight to exactly zero, and those zeros enter the average. The 0.365 for intrusion-dreams comes from an edge included in every sweep. The 0.062 for intrusion-upset mixes zeros from the third of sweeps that excluded the edge with nonzero values from the rest, so it is smaller than the typical weight in the sweeps that included it. Read it as the estimate under structure uncertainty, not as “the effect, given that the edge exists”. See Model Output.

The diagnostics account for that mixture. The Monte Carlo error of a model-averaged weight has two sources, how well the chain resolved whether the edge is in, and how well it resolved the weight when it was in, and mcse and n_eff combine both. The extra column, share_incl, says how much of the squared Monte Carlo error comes from the inclusion side: near 1 the bottleneck is resolving inclusion, near 0 it is resolving the weight. The 0.928 for intrusion-upset says its error is almost all inclusion uncertainty, while the 0.000 for intrusion-dreams reflects an edge that never left the model, so only the weight contributes. share_incl is a bookkeeping column, not a verdict: it tells you what a longer run would be spending its iterations on.

Is this block okay? The same two guidelines apply: across all 136 pairs the smallest n_eff is 1646 and the largest Rhat is 1.005, so every model-averaged weight here is estimated precisely enough to report.

Inclusion probabilities

One row per pair again, and mean is now the posterior inclusion probability: how plausible it is, after seeing the data, that the edge is in the network. A value of 0.99 is strong evidence for a direct association; a value near 0 is evidence against one; 0.5 means the data have not decided, which is a finding, not a failure. It is estimated by averaging the Rao-Blackwellized inclusion draws rather than the raw 0/1 indicator draws, which targets the same probability with less Monte Carlo noise. For evidence claims, turn probabilities into inclusion Bayes factors with extract_inclusion_bf(); see Edge Selection for the thresholds and their reading.

mcse, sd, n_eff, and Rhat are the familiar columns, computed on the Rao-Blackwellized draws. The two new columns watch the sampler’s exploration of structures: n0->1 and n1->0 count how often the edge indicator switched state across the four chains. Read them beside n_eff: n_eff says the probability estimate is precise, the counts say the chain actually moved between structures, and only together do they say converged and explored. See MCMC Diagnostics.

The blanks are the note’s suppressed NA values. Nine edges here, intrusion-dreams among them, were included in every sweep of every chain: their inclusion probabilities are pinned at 1, and a chain that never moves gives the precision estimates nothing to work with. The package reports NA for mcse and n_eff on those rows (and for all but one of them Rhat too) rather than print reassuring-looking numbers. The mean stays, and for these saturated edges the evidence verdict rests on extract_inclusion_bf(), which remains finite where the probability rounds to 1.

Is this block okay? Guidelines as before, applied to the non-blank rows: here the smallest n_eff is 995 and every Rhat is at most 1.006. Blank rows are not failures; they are the edges the data decided most emphatically.

Where is X?

The printout closes by pointing beyond itself: full tables via summary(fit)$<component>, and other scales via the extractor functions. The numbers users look for that are not in summary():

You want	Where it lives
Inclusion Bayes factors	`extract_inclusion_bf(fit)`; finite even for saturated edges, and `log = TRUE` for evidence too strong to print
Associations on an interpretable scale	`extract_log_odds(fit)` for the discrete pairs of this fit; `extract_partial_correlations(fit)` for continuous and mixed fits (it returns `NULL` on an all-ordinal fit like this one)
Posterior means as matrices, for plotting	`coef(fit)`: `$pairwise`, `$indicator`, and `$main`; see Methods
Raw posterior draws, for intervals and custom summaries	`fit$raw_samples`, or the extractor functions such as `extract_pairwise_interactions(fit)`

What to report

For each association you interpret: the posterior mean partial association with its posterior uncertainty, as in “intrusion-dreams 0.37 (posterior SD 0.04)”, or a credible interval from the raw draws.
For each edge claim: the posterior inclusion probability together with the inclusion Bayes factor from extract_inclusion_bf(), against the evidence threshold you state.
One sampling-adequacy sentence: chains and iterations, the largest Rhat, and the smallest effective sample size across the summary tables, as in “four chains of 10,000 iterations; all Rhat at or below 1.01, smallest effective sample size 995”.
Not licensed by this output: an inclusion probability below 0.5 is not evidence that the edge is absent. Absence is a positive claim, and it needs the Bayes factor from extract_inclusion_bf() to clear your threshold in the other direction.

Where to next

Model Output for the full output structure, and MCMC Diagnostics for why the diagnostic guidelines are what they are.
Methods for the summary(), print(), and coef() contract, and Extractor Functions for everything the closing table points at.
Vehtari et al. (2021) for the ESS and Rhat guidelines, and Sekulovski et al. (2024) for testing edge presence and absence with inclusion Bayes factors, as the citable record.

A bgmCompare() fit prints more tables than this, adding group differences and their own inclusion probabilities; see Group Comparison until the comparing-groups example page ships.

References

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

Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: An improved $\widehat{R}$ for assessing convergence of MCMC (with discussion). Bayesian Analysis, 16(2), 667–718. https://doi.org/10.1214/20-BA1221