# Model Comparison with bgmCompare

# 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