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.799 0.007 0.364 2903.156    1
#> 2 closeatt (1) -2.403 0.006 0.343 3738.953    1
#> 3 distract (1) -0.793 0.005 0.258 2717.506    1
#> 4   forget (1) -1.795 0.005 0.303 3081.048    1
#> 5 instruct (1) -2.645 0.007 0.357 2893.943    1
#> 
#> Pairwise interactions:
#>           parameter   mean  mcse    sd    n_eff Rhat
#> 1    avoid-closeatt  0.504 0.005 0.215 1681.491    1
#> 2    avoid-distract  0.876 0.003 0.171 3399.601    1
#> 3      avoid-forget  0.272 0.004 0.178 1820.401    1
#> 4    avoid-instruct  0.238 0.003 0.209 3723.882    1
#> 5 closeatt-distract -0.083 0.004 0.196 3008.148    1
#> 6   closeatt-forget  0.097 0.002 0.142 4193.970    1
#> ... (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.727 0.014 0.396  808.518     1   307   307
#>  avoid-distract (pairwise) 0.249 0.005 0.298 4252.905     1   765   764
#>    avoid-forget (pairwise) 0.785 0.013 0.367   811.04 1.004   262   261
#>  avoid-instruct (pairwise)     1     0 0.007 3119.822     1     0     0
#>            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.895 0.010 0.544 3169.142          0 1.000
#>  closeatt (diff1; 1) -2.157 0.011 0.572 2727.445          0 1.000
#>  distract (diff1; 1) -1.830 0.010 0.486 2320.306          0 1.000
#>    forget (diff1; 1) -2.174 0.010 0.506 2578.157          0 1.000
#>  instruct (diff1; 1) -1.602 0.015 0.635 1835.366          0 1.001
#> 
#> Group differences (pairwise effects):
#>                  parameter   mean  mcse    sd    n_eff share_incl  Rhat
#>     avoid-closeatt (diff1)  0.536 0.012 0.458 1563.007      0.788 1.000
#>     avoid-distract (diff1)  0.028 0.002 0.161 4201.686      0.043 1.000
#>       avoid-forget (diff1)  0.550 0.010 0.403 1553.870      0.787 1.003
#>     avoid-instruct (diff1) -1.476 0.007 0.402 3257.192      0.001 1.000
#>  closeatt-distract (diff1) -0.254 0.009 0.368 1815.316      0.476 1.002
#>    closeatt-forget (diff1)  0.050 0.003 0.200 3915.938      0.145 1.000
#> ... (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.7986291 -1.895160
#> closeatt(c1) -2.4030223 -2.157442
#> distract(c1) -0.7932048 -1.830019
#> forget(c1)   -1.7947925 -2.173816
#> instruct(c1) -2.6445981 -1.602313
#> 
#> $pairwise_effects_raw
#>                      baseline       diff1
#> avoid-closeatt     0.50398981  0.53556359
#> avoid-distract     0.87609885  0.02817147
#> avoid-forget       0.27211734  0.55013711
#> avoid-instruct     0.23775527 -1.47642687
#> closeatt-distract -0.08310233 -0.25388949
#> closeatt-forget    0.09711549  0.04969847
#> closeatt-instruct  0.80698877  0.19239735
#> distract-forget    0.26610463  0.02789478
#> distract-instruct  0.68556872  0.42890877
#> forget-instruct    0.57825160  0.27139793
#> 
#> $main_effects_groups
#>                  group1    group2
#> avoid(c1)    -1.8510493 -3.746209
#> closeatt(c1) -1.3243011 -3.481744
#> distract(c1)  0.1218046 -1.708214
#> forget(c1)   -0.7078845 -2.881700
#> instruct(c1) -1.8434416 -3.445755
#> 
#> $pairwise_effects_groups
#>                         group1     group2
#> avoid-closeatt     0.236208011  0.7717716
#> avoid-distract     0.862013114  0.8901846
#> avoid-forget      -0.002951213  0.5471859
#> avoid-instruct     0.975968704 -0.5004582
#> closeatt-distract  0.043842414 -0.2100471
#> closeatt-forget    0.072266259  0.1219647
#> closeatt-instruct  0.710790094  0.9031874
#> distract-forget    0.252157243  0.2800520
#> distract-instruct  0.471114331  0.9000231
#> forget-instruct    0.442552636  0.7139506
#> 
#> $indicators
#>            avoid closeatt distract  forget instruct
#> avoid    1.00000  0.72575  0.24900 0.78175  1.00000
#> closeatt 0.72575  1.00000  0.47950 0.27875  0.44650
#> distract 0.24900  0.47950  1.00000 0.23175  0.71625
#> forget   0.78175  0.27875  0.23175 1.00000  0.55525
#> instruct 1.00000  0.44650  0.71625 0.55525  1.00000

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