Getting Started with bgms

Introduction

The bgms package implements Bayesian methods for analyzing graphical models of binary and ordinal variables. It estimates main effects (category thresholds) and pairwise interactions in an ordinal Markov random field (MRF), with optional Bayesian edge selection via spike–and–slab priors. The package provides two main entry points:

This vignette walks through the basic workflow: fitting a model, summarizing posterior output, and visualizing results.

Wenchuan dataset

The dataset Wenchuan contains responses from survivors of the 2008 Wenchuan earthquake on posttraumatic stress items. Here, we analyze a subset of the first five items as a demonstration.

library(bgms)

# Analyse a subset of the Wenchuan dataset
?Wenchuan
data = Wenchuan[, 1:5]
head(data)
#>      intrusion dreams flash upset physior
#> [1,]         2      2     2     2       3
#> [2,]         2      2     2     3       3
#> [3,]         2      4     4     4       3
#> [4,]         2      1     2     2       1
#> [5,]         2      2     2     2       2
#> [6,]         4      3     2     2       2

Fitting a model

The main entry point is bgm() for single-group models and bgmCompare() for multiple-group comparisons.

fit = bgm(data, seed = 1234)

Posterior summaries

summary(fit)
#> Posterior summaries from Bayesian estimation:
#> 
#> Category thresholds:
#>                 mean  mcse    sd   n_eff  Rhat
#> intrusion (1)  0.502 0.008 0.243 854.956 1.000
#> intrusion (2) -1.848 0.018 0.363 400.710 1.003
#> intrusion (3) -4.748 0.034 0.583 296.020 1.002
#> intrusion (4) -9.358 0.052 0.940 322.463 1.001
#> dreams (1)    -0.595 0.006 0.191 861.969 1.006
#> dreams (2)    -3.799 0.013 0.348 711.950 1.006
#> ... (use `summary(fit)$main` to see full output)
#> 
#> Pairwise interactions:
#>                    mean    sd  mcse    n_eff  Rhat
#> intrusion-dreams  0.630 0.002 0.068  885.215 1.002
#> intrusion-flash   0.339 0.002 0.065 1105.927 1.002
#> intrusion-upset   0.180 0.083 0.009   88.374 1.007
#> intrusion-physior 0.195 0.072 0.006  161.937 1.000
#> dreams-flash      0.500 0.001 0.061 1726.881 1.000
#> dreams-upset      0.231 0.002 0.055  613.238 1.003
#> ... (use `summary(fit)$pairwise` 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` still contains the NA values.
#> 
#> Inclusion probabilities:
#>                    mean    sd  mcse n0->0 n0->1 n1->0 n1->1  n_eff  Rhat
#> intrusion-dreams  1.000 0.000           0     0     0  1999             
#> intrusion-flash   1.000 0.000           0     0     0  1999             
#> intrusion-upset   0.886 0.318 0.042   217    11    11  1760 55.981 1.071
#> intrusion-physior 0.941 0.236 0.026   109     9     9  1872  84.48 1.047
#> dreams-flash      1.000 0.000           0     0     0  1999             
#> dreams-upset      1.000 0.000           0     0     0  1999             
#> ... (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.
#> 
#> Use `summary(fit)$<component>` to access full results.
#> See the `easybgm` package for other summary and plotting tools.

You can also access posterior means or inclusion probabilities directly:

coef(fit)
#> $main
#>              cat (1)   cat (2)   cat (3)    cat (4)
#> intrusion  0.5024739 -1.848388 -4.748210  -9.358041
#> dreams    -0.5948872 -3.799160 -7.126397 -11.571395
#> flash     -0.1090829 -2.579439 -5.389731  -9.712605
#> upset      0.4213477 -1.292089 -3.358005  -7.007588
#> physior   -0.6090363 -3.170423 -6.228609 -10.577326
#> 
#> $pairwise
#>           intrusion      dreams      flash      upset     physior
#> intrusion 0.0000000 0.629828116 0.33879204 0.18013791 0.194800560
#> dreams    0.6298281 0.000000000 0.49962958 0.23067829 0.006216451
#> flash     0.3387920 0.499629583 0.00000000 0.01100187 0.306863313
#> upset     0.1801379 0.230678286 0.01100187 0.00000000 0.714164491
#> physior   0.1948006 0.006216451 0.30686331 0.71416449 0.000000000
#> 
#> $indicator
#>           intrusion dreams  flash  upset physior
#> intrusion     0.000  1.000 1.0000 0.8860   0.941
#> dreams        1.000  0.000 1.0000 1.0000   0.062
#> flash         1.000  1.000 0.0000 0.0985   1.000
#> upset         0.886  1.000 0.0985 0.0000   1.000
#> physior       0.941  0.062 1.0000 1.0000   0.000

Network plot

To visualize the network structure, we threshold the posterior inclusion probabilities at 0.5 and plot the resulting adjacency matrix.

library(qgraph)

median_probability_network = coef(fit)$pairwise
median_probability_network[coef(fit)$indicator < 0.5] = 0.0

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

Next steps