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.
ADHD dataset
We illustrate with a subset from the ADHD dataset
included in bgms.
library(bgms)
?ADHD
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]
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 avoid (1) -2.676 0.010 0.389 1411.000 1.000
#> 2 closeatt (1) -2.249 0.011 0.376 1198.435 1.002
#> 3 distract (1) -0.486 0.013 0.337 637.870 1.001
#> 4 forget (1) -1.580 0.010 0.330 1070.957 1.000
#> 5 instruct (1) -2.416 0.014 0.401 797.064 1.008
#>
#> Pairwise interactions:
#> parameter mean mcse sd n_eff Rhat
#> 1 avoid-closeatt 0.993 0.017 0.459 723.538 1.000
#> 2 avoid-distract 1.704 0.009 0.365 1598.070 1.003
#> 3 avoid-forget 0.486 0.014 0.381 709.989 1.009
#> 4 avoid-instruct 0.387 0.014 0.464 1029.911 1.001
#> 5 closeatt-distract -0.257 0.011 0.393 1238.404 1.001
#> 6 closeatt-forget 0.146 0.008 0.312 1513.112 1.003
#> ... (use `summary(fit)$pairwise` to see full output)
#>
#> Inclusion probabilities:
#> parameter mean sd mcse n0->0 n0->1 n1->0 n1->1 n_eff
#> avoid (main) 1.000 0.000 0 0 0 1999
#> avoid-closeatt (pairwise) 0.821 0.383 0.015 209 149 148 1493 678.075
#> avoid-distract (pairwise) 0.391 0.488 0.012 779 438 438 344 1703.716
#> avoid-forget (pairwise) 0.819 0.385 0.017 244 118 118 1519 496.957
#> avoid-instruct (pairwise) 1.000 0.022 0 0 1 1 1997 2002.003
#> closeatt (main) 1.000 0.000 0 0 0 1999
#> Rhat
#>
#> 1
#> 1
#> 1.031
#> 1.291
#>
#> ... (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
#> avoid (diff1; 1) -2.523 0.728 1.000
#> closeatt (diff1; 1) -3.008 0.726 1.002
#> distract (diff1; 1) -2.552 0.689 1.001
#> forget (diff1; 1) -2.831 0.654 1.007
#> instruct (diff1; 1) -2.374 0.884 1.003
#> 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
#> avoid-closeatt (diff1) 1.296 0.911 0.031 890.715 1.000
#> avoid-distract (diff1) 0.219 0.362 0.011 1018.055 1.000
#> avoid-forget (diff1) 1.235 0.831 0.032 676.485 1.006
#> avoid-instruct (diff1) -2.805 0.970 0.036 740.334 1.004
#> closeatt-distract (diff1) -0.174 0.348 0.012 795.120 1.000
#> closeatt-forget (diff1) 0.169 0.333 0.012 715.853 1.000
#> ... (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
#> avoid(c1) -2.6758947 -2.523160
#> closeatt(c1) -2.2489431 -3.007718
#> distract(c1) -0.4859767 -2.551852
#> forget(c1) -1.5803152 -2.830772
#> instruct(c1) -2.4164664 -2.373887
#>
#> $pairwise_effects_raw
#> baseline diff1
#> avoid-closeatt 0.9927303 1.2964029
#> avoid-distract 1.7035207 0.2190339
#> avoid-forget 0.4858398 1.2349851
#> avoid-instruct 0.3870157 -2.8050136
#> closeatt-distract -0.2573117 -0.1735373
#> closeatt-forget 0.1456844 0.1691829
#> closeatt-instruct 1.5675969 0.6163591
#> distract-forget 0.4044984 0.2301470
#> distract-instruct 1.2586856 1.2586632
#> forget-instruct 1.1314787 0.8426717
#>
#> $main_effects_groups
#> group1 group2
#> avoid(c1) -1.4143149 -3.937475
#> closeatt(c1) -0.7450839 -3.752802
#> distract(c1) 0.7899495 -1.761903
#> forget(c1) -0.1649295 -2.995701
#> instruct(c1) -1.2295228 -3.603410
#>
#> $pairwise_effects_groups
#> group1 group2
#> avoid-closeatt 0.34452892 1.6409318
#> avoid-distract 1.59400382 1.8130377
#> avoid-forget -0.13165272 1.1033324
#> avoid-instruct 1.78952252 -1.0154911
#> closeatt-distract -0.17054303 -0.3440803
#> closeatt-forget 0.06109296 0.2302758
#> closeatt-instruct 1.25941733 1.8757764
#> distract-forget 0.28942488 0.5195718
#> distract-instruct 0.62935398 1.8880172
#> forget-instruct 0.71014281 1.5528145
#>
#> $indicators
#> avoid closeatt distract forget instruct
#> avoid 1.0000 0.8210 0.3915 0.8190 0.9995
#> closeatt 0.8210 1.0000 0.3745 0.3915 0.6030
#> distract 0.3915 0.3745 1.0000 0.3970 0.8440
#> forget 0.8190 0.3915 0.3970 1.0000 0.7485
#> instruct 0.9995 0.6030 0.8440 0.7485 1.0000
Visualizing group networks
We can use the output to plot the network for the ADHD group:
library(qgraph)
adhd_network = matrix(0, 5, 5)
adhd_network[lower.tri(adhd_network)] = coef(fit)$pairwise_effects_groups[, 1]
adhd_network = adhd_network + t(adhd_network)
colnames(adhd_network) = colnames(data_adhd)
rownames(adhd_network) = colnames(data_adhd)
qgraph(adhd_network,
theme = "TeamFortress",
maximum = 1,
fade = FALSE,
color = c("#f0ae0e"), vsize = 10, repulsion = .9,
label.cex = 1, label.scale = "FALSE",
labels = colnames(data_adhd)
)
