bgm.em {bgms} R Documentation

EM variable selection for a Markov Random Field model for ordinal variables.

Description

The function bgm.em selects promising edges for the ordinal MRF using the joint pseudolikelihood and a continuous spike and slab prior distribution stipulated on the MRF's interaction or association parameters.

Usage

bgm.em(
x,
precision = 0.975,

Value

A list containing:

• interactions: A matrix with p rows and p columns, containing the pairwise association estimates in the off-diagonal elements.

• gamma: A matrix with p rows and p columns, containing the expected values of edge inclusion variables (local posterior probabilities of edge inclusion).

• thresholds: A matrix with p rows and max(m) columns, containing the category thresholds for each node.

• theta (if hierarchical == TRUE): A numeric value representing the modal estimate of the prior inclusion probability.

Examples


#Store user par() settings

##Analyse the Wenchuan dataset
fit = bgm.em(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 = "#bfbfbf", cex = 1.3)
abline(h = 0, lty = 2, col = "#bfbfbf")
abline(h = 1, lty = 2, col = "#bfbfbf")
abline(h = .5, lty = 2, col = "#bfbfbf")
mtext("Posterior Inclusion Probability", side = 1, line = 3, cex = 1.7)
mtext("Posterior Mode Edge Weight", side = 2, line = 3, cex = 1.7)
axis(1)
axis(2, las = 1)

#------------------------------------------------------------------------------|
# THE LOCAL MEDIAN PROBABILITY NETWORK
#------------------------------------------------------------------------------|

library(qgraph) #For plotting the estimated network

posterior.inclusion = fit$gamma[lower.tri(fit$gamma)]
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)

qgraph(median.prob.model,
theme = "TeamFortress",
maximum = .5,