bgm {bgms}  R Documentation 
The function bgm
explores the joint pseudoposterior distribution of
structures and parameters in a Markov Random Field for mixed binary and
ordinal variables.
bgm(
x,
iter = 10000,
burnin = 1000,
interaction_prior = c("Cauchy", "UnitInfo"),
cauchy_scale = 2.5,
edge_prior = c("Bernoulli", "BetaBernoulli"),
inclusion_probability = 0.5,
beta_bernoulli_alpha = 1,
beta_bernoulli_beta = 1,
threshold_alpha = 0.5,
threshold_beta = 0.5,
adaptive = FALSE,
na.action = c("listwise", "impute"),
save = FALSE,
display_progress = TRUE
)
x 
A data frame or matrix with 
iter 
The number of iterations of the Gibbs sampler. The default of

burnin 
The number of iterations of the Gibbs sampler before its output is saved. Since it may take some time for the Gibbs sampler to converge to the posterior distribution, it is recommended not to set this number too low. 
interaction_prior 
The type of prior distribution that is used for the
interaction effects. Currently, two prior densities are implemented: The
Cauchy prior ( 
cauchy_scale 
The scale of the Cauchy prior for interactions. Defaults
to 
edge_prior 
The prior distribution for the edges or structure of the
network. Two prior distributions are currently implemented: The Bernoulli
model 
inclusion_probability 
The prior edge inclusion probability for the
Bernoulli model. Can be a single probability, or a matrix of 
beta_bernoulli_alpha , beta_bernoulli_beta 
The two shape parameters of
the Beta prior density for the Bernoulli inclusion probability. Must be
positive numbers. Defaults to 
threshold_alpha , threshold_beta 
The shape parameters of the betaprime
prior density for the threshold parameters. Must be positive values. If the
two values are equal, the prior density is symmetric about zero. If

adaptive 
Should the function use an adaptive Metropolis algorithm to
update interaction parameters within the model? If 
na.action 
How do you want the function to handle missing data? If

save 
Should the function collect and return all samples from the Gibbs
sampler ( 
display_progress 
Should the function show a progress bar
( 
A discrete spike and slab prior distribution is stipulated on the pairwise interactions. By formulating it as a mixture of mutually singular distributions, the function can use a combination of MetropolisHastings and Gibbs sampling to create a Markov chain that has the joint posterior distribution as invariant. Current options for the slab distribution are the unitinformation prior or a Cauchy with an optional scaling parameter. A Betaprime distribution is used for the exponent of the category parameters. A uniform prior is used for edge inclusion variables (i.e., the prior probability that the edge is included is 0.5).
If save = FALSE
(the default), the result is a list of class
“bgms” containing the following matrices:
gamma
: A matrix with p
rows and p
columns,
containing posterior inclusion probabilities of individual edges.
interactions
: A matrix with p
rows and p
columns,
containing modelaveraged posterior means of the pairwise associations.
thresholds
: A matrix with p
rows and max(m)
columns, containing modelaveraged category thresholds.
If save = TRUE
, the result is a list of class “bgms” containing:
gamma
: A matrix with iter
rows and
p * (p  1) / 2
columns, containing the edge inclusion indicators from
every iteration of the Gibbs sampler.
interactions
: A matrix with iter
rows and
p * (p  1) / 2
columns, containing parameter states from every
iteration of the Gibbs sampler for the pairwise associations.
thresholds
: A matrix with iter
rows and
sum(m)
columns, containing parameter states from every iteration of
the Gibbs sampler for the category thresholds.
Column averages of these matrices provide the modelaveraged posterior means.
In addition to the analysis results, the bgm output lists some of the arguments of its call. This is useful for postprocessing the results.
#Store user par() settings
op < par(no.readonly = TRUE)
##Analyse the Wenchuan dataset
# Here, we use 1e4 iterations, for an actual analysis please use at least
# 1e5 iterations.
fit = bgm(x = Wenchuan)
#
# INCLUSION  EDGE WEIGHT PLOT
#
par(mar = c(6, 5, 1, 1))
plot(x = fit$interactions[lower.tri(fit$interactions)],
y = fit$gamma[lower.tri(fit$gamma)], ylim = c(0, 1),
xlab = "", ylab = "", axes = FALSE, pch = 21, bg = "gray", cex = 1.3)
abline(h = 0, lty = 2, col = "gray")
abline(h = 1, lty = 2, col = "gray")
abline(h = .5, lty = 2, col = "gray")
mtext("Posterior Mode Edge Weight", side = 1, line = 3, cex = 1.7)
mtext("Posterior Inclusion Probability", side = 2, line = 3, cex = 1.7)
axis(1)
axis(2, las = 1)
#
# EVIDENCE  EDGE WEIGHT PLOT
#
#The bgms package currently assumes that the prior odds are 1:
prior.odds = 1
posterior.inclusion = fit$gamma[lower.tri(fit$gamma)]
posterior.odds = posterior.inclusion / (1  posterior.inclusion)
log.bayesfactor = log(posterior.odds / prior.odds)
log.bayesfactor[log.bayesfactor > 5] = 5
par(mar = c(5, 5, 1, 1) + 0.1)
plot(fit$interactions[lower.tri(fit$interactions)], log.bayesfactor, pch = 21, bg = "#bfbfbf",
cex = 1.3, axes = FALSE, xlab = "", ylab = "", ylim = c(5, 5.5),
xlim = c(0.5, 1.5))
axis(1)
axis(2, las = 1)
abline(h = log(1/10), lwd = 2, col = "#bfbfbf")
abline(h = log(10), lwd = 2, col = "#bfbfbf")
text(x = 1, y = log(1 / 10), labels = "Evidence for Exclusion", pos = 1,
cex = 1.7)
text(x = 1, y = log(10), labels = "Evidence for Inclusion", pos = 3, cex = 1.7)
text(x = 1, y = 0, labels = "Absence of Evidence", cex = 1.7)
mtext("LogInclusion Bayes Factor", side = 2, line = 3, cex = 1.5, las = 0)
mtext("Posterior Mean Interactions ", side = 1, line = 3.7, cex = 1.5, las = 0)
#
# THE LOCAL MEDIAN PROBABILITY NETWORK
#
tmp = fit$interactions[lower.tri(fit$interactions)]
tmp[posterior.inclusion < 0.5] = 0
median.prob.model = matrix(0, nrow = ncol(Wenchuan), ncol = ncol(Wenchuan))
median.prob.model[lower.tri(median.prob.model)] = tmp
median.prob.model = median.prob.model + t(median.prob.model)
rownames(median.prob.model) = colnames(Wenchuan)
colnames(median.prob.model) = colnames(Wenchuan)
library(qgraph)
qgraph(median.prob.model,
theme = "TeamFortress",
maximum = .5,
fade = FALSE,
color = c("#f0ae0e"), vsize = 10, repulsion = .9,
label.cex = 1.1, label.scale = "FALSE",
labels = colnames(Wenchuan))
#Restore user par() settings
par(op)