The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Model Comparison with bgmCompare

Introduction

The function bgmCompare() extends bgm() to independent-sample designs. It estimates whether edge weights and category thresholds differ across groups in an ordinal Markov random field (MRF).

Posterior inclusion probabilities indicate how plausible it is that a group difference exists in a given parameter. These can be converted to Bayes factors for hypothesis testing.

Boredom dataset

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

library(bgms)

?Boredom
data_french = Boredom[Boredom$language == "fr", -1]
data_french = data_french[, 1:5]
data_english = Boredom[Boredom$language != "fr", -1]
data_english = data_english[, 1:5]

Fitting a model

fit = bgmCompare(x = data_french, y = data_english, seed = 1234)

Note: During fitting, progress bars are shown in interactive sessions. In this vignette, they are suppressed for clarity. Sampling can take a while; the progress bars usually help track progress.

Posterior summaries

The summary shows both baseline effects and group differences:

summary(fit)
#> Posterior summaries from Bayesian grouped MRF estimation (bgmCompare):
#> 
#> Category thresholds:
#>        parameter    mean  mcse    sd    n_eff  Rhat
#> 1 loose_ends (1)  -0.938 0.002 0.096 3415.827 1.001
#> 2 loose_ends (2)  -2.518 0.003 0.139 2069.362 1.004
#> 3 loose_ends (3)  -3.796 0.005 0.183 1543.211 1.005
#> 4 loose_ends (4)  -5.082 0.007 0.243 1236.538 1.006
#> 5 loose_ends (5)  -7.609 0.010 0.335 1164.760 1.007
#> 6 loose_ends (6) -10.106 0.014 0.452 1092.310 1.006
#> ... (use `summary(fit)$main` to see full output)
#> 
#> Pairwise interactions:
#>                parameter  mean mcse    sd    n_eff  Rhat
#> 1   loose_ends-entertain 0.170    0 0.013 2245.500 1.002
#> 2  loose_ends-repetitive 0.058    0 0.012 2059.968 1.007
#> 3 loose_ends-stimulation 0.126    0 0.012 2648.651 1.002
#> 4   loose_ends-motivated 0.141    0 0.013 2490.105 1.002
#> 5   entertain-repetitive 0.067    0 0.012 2712.597 1.000
#> 6  entertain-stimulation 0.109    0 0.012 2591.243 1.001
#> ... (use `summary(fit)$pairwise` to see full output)
#> 
#> Inclusion probabilities:
#>                          parameter  mean    sd  mcse n0->0 n0->1 n1->0 n1->1
#>                  loose_ends (main) 0.009 0.094 0.003  3947    16    16    20
#>    loose_ends-entertain (pairwise) 0.020 0.139 0.002  3845    75    75     4
#>   loose_ends-repetitive (pairwise) 0.507 0.500 0.018  1640   331   330  1698
#>  loose_ends-stimulation (pairwise) 0.042 0.200 0.003  3680   152   152    15
#>    loose_ends-motivated (pairwise) 0.018 0.132 0.002  3860    68    68     3
#>                   entertain (main) 0.001 0.039 0.001  3988     5     5     1
#>     n_eff  Rhat
#>  1156.240 1.027
#>  3755.694 1.001
#>   792.304 1.002
#>  3617.928 1.012
#>  3805.320 1.004
#>  2864.511 1.081
#> ... (use `summary(fit)$indicator` to see full output)
#> Note: NA values are suppressed in the print table. They occur when an indicator
#> was constant (all 0 or all 1) across all iterations, so sd/mcse/n_eff/Rhat
#> are undefined; `summary(fit)$indicator` still contains the NA values.
#> 
#> Group differences (main effects):
#>              parameter   mean    sd  mcse    n_eff  Rhat
#>  loose_ends (diff1; 1)  0.000 0.002 0.000   47.557 1.000
#>  loose_ends (diff1; 2)  0.004 0.044 0.001 1082.785 1.001
#>  loose_ends (diff1; 3)  0.004 0.043 0.001  902.638 1.001
#>  loose_ends (diff1; 4)  0.004 0.042 0.001  923.443 1.001
#>  loose_ends (diff1; 5)  0.001 0.010 0.001  113.681 1.000
#>  loose_ends (diff1; 6) -0.003 0.031 0.001  630.010 1.000
#> ... (use `summary(fit)$main_diff` to see full output)
#> Note: NA values are suppressed in the print table. They occur here when an
#> indicator was zero across all iterations, so mcse/n_eff/Rhat are undefined;
#> `summary(fit)$main_diff` still contains the NA values.
#> 
#> Group differences (pairwise effects):
#>                       parameter  mean    sd  mcse    n_eff  Rhat
#>    loose_ends-entertain (diff1) 0.000 0.000 0.000   92.057 1.001
#>   loose_ends-repetitive (diff1) 0.021 0.022 0.001  837.278 1.001
#>  loose_ends-stimulation (diff1) 0.001 0.005 0.000 2789.970 1.001
#>    loose_ends-motivated (diff1) 0.000 0.002 0.000 1517.163 1.000
#>    entertain-repetitive (diff1) 0.004 0.011 0.000 1174.633 1.001
#>   entertain-stimulation (diff1) 0.002 0.007 0.000 2567.839 1.001
#> ... (use `summary(fit)$pairwise_diff` to see full output)
#> Note: NA values are suppressed in the print table. They occur here when an
#> indicator was zero across all iterations, so mcse/n_eff/Rhat are undefined;
#> `summary(fit)$pairwise_diff` still contains the NA values.
#> 
#> 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
#> loose_ends(c1)   -0.9376669 -8.550637e-05
#> loose_ends(c2)   -2.5179548  4.207447e-03
#> loose_ends(c3)   -3.7960694  4.109076e-03
#> loose_ends(c4)   -5.0822845  4.002782e-03
#> loose_ends(c5)   -7.6089748  9.141143e-04
#> loose_ends(c6)  -10.1063777 -2.905428e-03
#> entertain(c1)    -0.8692630 -4.662515e-04
#> entertain(c2)    -2.2458877 -1.496695e-04
#> entertain(c3)    -3.8177841  2.483365e-04
#> entertain(c4)    -5.1697251 -3.493918e-04
#> entertain(c5)    -7.0455504 -2.961615e-05
#> entertain(c6)    -9.5780234  8.034366e-04
#> repetitive(c1)   -0.1379151 -8.806321e-03
#> repetitive(c2)   -0.6611929 -1.585117e-02
#> repetitive(c3)   -1.0866300 -5.037304e-03
#> repetitive(c4)   -1.8627494  6.470289e-03
#> repetitive(c5)   -3.2512960  2.020118e-02
#> repetitive(c6)   -5.0299108  2.038924e-02
#> stimulation(c1)  -0.5457673 -3.167586e-03
#> stimulation(c2)  -1.7894587 -1.004249e-03
#> stimulation(c3)  -2.5304682 -1.332632e-03
#> stimulation(c4)  -3.6190432 -2.226982e-03
#> stimulation(c5)  -5.0975023 -1.200563e-03
#> stimulation(c6)  -7.0661001 -3.494762e-03
#> motivated(c1)    -0.5585283 -2.017340e-03
#> motivated(c2)    -1.8120817 -1.070582e-03
#> motivated(c3)    -3.2996542  1.534816e-03
#> motivated(c4)    -4.7945890  2.748786e-03
#> motivated(c5)    -6.7933940 -6.137974e-04
#> motivated(c6)    -9.1663037  1.980850e-03
#> 
#> $pairwise_effects_raw
#>                          baseline         diff1
#> loose_ends-entertain   0.16991838  2.642296e-05
#> loose_ends-repetitive  0.05811138  2.137907e-02
#> loose_ends-stimulation 0.12622492  9.635558e-04
#> loose_ends-motivated   0.14097684 -2.398129e-04
#> entertain-repetitive   0.06659031  4.431347e-03
#> entertain-stimulation  0.10937128  1.801550e-03
#> entertain-motivated    0.08666756  3.824078e-03
#> repetitive-stimulation 0.05631314  6.388122e-04
#> repetitive-motivated   0.13750770  1.397853e-02
#> stimulation-motivated  0.10798919  1.154900e-04
#> 
#> $main_effects_groups
#>                      group1      group2
#> loose_ends(c1)   -0.9376242  -0.9377097
#> loose_ends(c2)   -2.5200585  -2.5158511
#> loose_ends(c3)   -3.7981240  -3.7940149
#> loose_ends(c4)   -5.0842859  -5.0802831
#> loose_ends(c5)   -7.6094318  -7.6085177
#> loose_ends(c6)  -10.1049250 -10.1078305
#> entertain(c1)    -0.8690299  -0.8694961
#> entertain(c2)    -2.2458129  -2.2459625
#> entertain(c3)    -3.8179083  -3.8176599
#> entertain(c4)    -5.1695504  -5.1698998
#> entertain(c5)    -7.0455356  -7.0455652
#> entertain(c6)    -9.5784251  -9.5776217
#> repetitive(c1)   -0.1335119  -0.1423182
#> repetitive(c2)   -0.6532674  -0.6691185
#> repetitive(c3)   -1.0841114  -1.0891487
#> repetitive(c4)   -1.8659846  -1.8595143
#> repetitive(c5)   -3.2613966  -3.2411954
#> repetitive(c6)   -5.0401054  -5.0197162
#> stimulation(c1)  -0.5441836  -0.5473511
#> stimulation(c2)  -1.7889565  -1.7899608
#> stimulation(c3)  -2.5298019  -2.5311345
#> stimulation(c4)  -3.6179297  -3.6201567
#> stimulation(c5)  -5.0969021  -5.0981026
#> stimulation(c6)  -7.0643527  -7.0678475
#> motivated(c1)    -0.5575196  -0.5595369
#> motivated(c2)    -1.8115464  -1.8126170
#> motivated(c3)    -3.3004216  -3.2988868
#> motivated(c4)    -4.7959634  -4.7932146
#> motivated(c5)    -6.7930871  -6.7937009
#> motivated(c6)    -9.1672942  -9.1653133
#> 
#> $pairwise_effects_groups
#>                            group1     group2
#> loose_ends-entertain   0.16990517 0.16993159
#> loose_ends-repetitive  0.04742184 0.06880091
#> loose_ends-stimulation 0.12574314 0.12670670
#> loose_ends-motivated   0.14109674 0.14085693
#> entertain-repetitive   0.06437464 0.06880598
#> entertain-stimulation  0.10847051 0.11027206
#> entertain-motivated    0.08475552 0.08857960
#> repetitive-stimulation 0.05599374 0.05663255
#> repetitive-motivated   0.13051843 0.14449696
#> stimulation-motivated  0.10793144 0.10804693
#> 
#> $indicators
#>             loose_ends entertain repetitive stimulation motivated
#> loose_ends     0.00900   0.01975    0.50725     0.04175   0.01775
#> entertain      0.01975   0.00150    0.13725     0.06650   0.11425
#> repetitive     0.50725   0.13725    0.03675     0.03175   0.36100
#> stimulation    0.04175   0.06650    0.03175     0.00575   0.01525
#> motivated      0.01775   0.11425    0.36100     0.01525   0.00700

Visualizing group networks

We can use the output to plot the network for the French sample:

library(qgraph)

french_network = matrix(0, 5, 5)
french_network[lower.tri(french_network)] = coef(fit)$pairwise_effects_groups[, 1]
french_network = french_network + t(french_network)
colnames(french_network) = colnames(data_french)
rownames(french_network) = colnames(data_french)

qgraph(french_network,
       theme = "TeamFortress",
       maximum = 1,
       fade = FALSE,
       color = c("#f0ae0e"), vsize = 10, repulsion = .9,
       label.cex = 1, label.scale = "FALSE",
       labels = colnames(data_french))

Next steps

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.