Introduction

The function bgmCompare() extends bgm() to independent-sample designs. It estimates whether pairwise interactions and category thresholds differ across groups in an ordinal Markov random field (Marsman et al., 2025), the member of the family that bgm() fits to binary and ordinal data.

Each difference is put under selection, so it carries a posterior inclusion probability and an inclusion Bayes factor of its own – the factor by which the data shift the odds that the groups differ on that parameter (Huth et al., 2023; Sekulovski et al., 2024). verdicts() reads those Bayes factors against a threshold and returns evidence of a difference, evidence of no difference, or undecided.

ADHD dataset

We illustrate with a subset from the ADHD dataset included in bgms.

library(bgms)

data_adhd = ADHD[ADHD$group == 1, -1]
data_adhd = data_adhd[, 1:5]
data_no_adhd = ADHD[ADHD$group == 0, -1]
data_no_adhd = data_no_adhd[, 1:5]

Fitting a model

fit = bgmCompare(x = data_adhd, y = data_no_adhd, seed = 1234)

The baseline pairwise prior is the same at both entry points: interaction_prior = normal_prior(1) in bgm() and in bgmCompare(). A bgmCompare() fit set beside two separate bgm() fits at their stated defaults therefore differs in the model for the groups, not in the prior on the baseline interactions. The group differences are priced separately, by difference_family ("Normal" by default) and difference_scale; these govern both the pairwise-interaction differences and the threshold differences, and interaction_prior does not reach them.

Groups often differ in which answer categories they actually use. bgmCompare() models the union of the categories observed across the groups, so a category one group uses is kept even when the other group never uses it; only a category value that no group uses at all is dropped. When a kept category is empty in some group, that group’s data say nothing about its threshold there, so the reported difference for that group-by-category combination comes from the prior and will be large and very uncertain. bgmCompare() warns and names each such case, and extract_arguments(fit)$category_support holds the per-group category counts. The pairwise (edge) parameters are unaffected. See ?bgmCompare for the details, including why Blume–Capel variables are exempt.

Posterior summaries

The summary shows both baseline effects and group differences:

summary(fit)
#> Posterior summaries from Bayesian grouped MRF estimation (bgmCompare):
#> 
#> groups: 1 = x (n = 146), 2 = y (n = 209)
#> 
#> Category thresholds:
#>      parameter   mean  mcse    sd    n_eff  Rhat
#> 1    avoid (1) -2.810 0.006 0.363 3920.182 1.000
#> 2 closeatt (1) -2.408 0.006 0.353 3695.700 1.000
#> 3 distract (1) -0.797 0.005 0.268 2806.912 1.000
#> 4   forget (1) -1.801 0.005 0.298 3134.105 1.001
#> 5 instruct (1) -2.649 0.007 0.352 2286.136 1.001
#> 
#> Pairwise interactions:
#>           parameter   mean  mcse    sd    n_eff  Rhat
#> 1    avoid-closeatt  0.500 0.006 0.219 1503.666 1.000
#> 2    avoid-distract  0.880 0.003 0.170 3507.710 1.000
#> 3      avoid-forget  0.274 0.004 0.182 1773.694 1.000
#> 4    avoid-instruct  0.242 0.004 0.209 2965.034 1.001
#> 5 closeatt-distract -0.079 0.004 0.195 3047.824 1.000
#> 6   closeatt-forget  0.098 0.002 0.144 4424.713 1.001
#> ... (use `summary(fit)$pairwise` to see full output)
#> 
#> Inclusion probabilities:
#>                  parameter  mean  mcse    sd    n_eff  Rhat n0->1 n1->0
#>               avoid (main)                                      0     0
#>  avoid-closeatt (pairwise)  0.72 0.014 0.404  799.092     1   301   300
#>  avoid-distract (pairwise) 0.256 0.005 0.298 4003.781     1   774   775
#>    avoid-forget (pairwise) 0.772 0.013 0.375  875.315 1.002   276   276
#>  avoid-instruct (pairwise) 0.999     0 0.024 2724.857 1.001     1     1
#>            closeatt (main)                                      0     0
#> ... (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.
#> 
#> Group differences (main effects):
#>            parameter   mean  mcse    sd    n_eff share_incl Rhat
#>     avoid (diff1; 1) -1.894 0.010 0.552 2884.173          0    1
#>  closeatt (diff1; 1) -2.151 0.011 0.575 2580.451          0    1
#>  distract (diff1; 1) -1.821 0.010 0.487 2189.023          0    1
#>    forget (diff1; 1) -2.137 0.011 0.505 2240.655          0    1
#>  instruct (diff1; 1) -1.583 0.016 0.648 1667.178          0    1
#> 
#> Group differences (pairwise effects):
#>                  parameter   mean  mcse    sd    n_eff share_incl  Rhat
#>     avoid-closeatt (diff1)  0.532 0.012 0.458 1524.633      0.822 1.000
#>     avoid-distract (diff1)  0.033 0.002 0.165 4386.372      0.057 1.000
#>       avoid-forget (diff1)  0.543 0.010 0.408 1672.416      0.801 1.000
#>     avoid-instruct (diff1) -1.476 0.007 0.406 2993.427      0.009 1.000
#>  closeatt-distract (diff1) -0.260 0.007 0.373 2528.300      0.647 1.000
#>    closeatt-forget (diff1)  0.052 0.003 0.202 3995.748      0.159 1.001
#> ... (use `summary(fit)$pairwise_diff` to see full output)
#> 
#> Use `summary(fit)$<component>` to access full results.
#> See the `easybgm` package for other summary and plotting tools.

You can extract posterior means and inclusion probabilities:

coef(fit)
#> $main_effects_raw
#>                baseline     diff1
#> avoid(c1)    -2.8096192 -1.894262
#> closeatt(c1) -2.4076644 -2.151033
#> distract(c1) -0.7974013 -1.821006
#> forget(c1)   -1.8006281 -2.136781
#> instruct(c1) -2.6488650 -1.582641
#> 
#> $pairwise_effects_raw
#>                      baseline       diff1
#> avoid-closeatt     0.49976998  0.53224721
#> avoid-distract     0.88006538  0.03259238
#> avoid-forget       0.27412877  0.54326548
#> avoid-instruct     0.24241175 -1.47642958
#> closeatt-distract -0.07912528 -0.26014631
#> closeatt-forget    0.09814889  0.05185378
#> closeatt-instruct  0.80534315  0.19141717
#> distract-forget    0.26455030  0.02286185
#> distract-instruct  0.68680669  0.43372707
#> forget-instruct    0.57667810  0.25842143
#> 
#> $main_effects_groups
#>                  group1    group2
#> avoid(c1)    -1.8624881 -3.756750
#> closeatt(c1) -1.3321478 -3.483181
#> distract(c1)  0.1131017 -1.707904
#> forget(c1)   -0.7322375 -2.869019
#> instruct(c1) -1.8575444 -3.440186
#> 
#> $pairwise_effects_groups
#>                        group1     group2
#> avoid-closeatt    0.233646382  0.7658936
#> avoid-distract    0.863769190  0.8963616
#> avoid-forget      0.002496037  0.5457615
#> avoid-instruct    0.980626542 -0.4958030
#> closeatt-distract 0.050947872 -0.2091984
#> closeatt-forget   0.072221999  0.1240758
#> closeatt-instruct 0.709634570  0.9010517
#> distract-forget   0.253119371  0.2759812
#> distract-instruct 0.469943153  0.9036702
#> forget-instruct   0.447467383  0.7058888
#> 
#> $indicators
#>            avoid closeatt distract  forget instruct
#> avoid    1.00000  0.71475  0.25775 0.77200   0.9995
#> closeatt 0.71475  1.00000  0.49050 0.28325   0.4295
#> distract 0.25775  0.49050  1.00000 0.23450   0.7205
#> forget   0.77200  0.28325  0.23450 1.00000   0.5380
#> instruct 0.99950  0.42950  0.72050 0.53800   1.0000

Difference verdicts are scale-contingent. bgmCompare() prices group differences on the association scale through difference_scale, and the calibration of that default is under study, so a verdict close to a decision threshold can move with the scale.

Visualizing the groups

plot() on a bgmCompare() fit draws the differences by default, split by what the data settle about each pair. type = "groups" draws each group’s own graph instead, on one shared layout, so the two pictures can be read side by side: line width is that group’s posterior mean association and colour carries its sign.

plot(fit, type = "groups")

Drawing needs the qgraph package, which bgms suggests rather than depends on. plot(fit) gives the difference display and plot(fit, type = "centrality") the posterior strength centrality of one group; see ?plot.bgmCompare.

Next steps

For a one-sample analysis, see the Getting Started vignette.
For diagnostics and convergence checks, see the Diagnostics vignette.
For additional analysis tools and more advanced plotting options, consider using the easybgm package, which integrates smoothly with bgms objects.

References

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

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