Package 'bgms'

Bayesian edge selection or Bayesian estimation for a Markov random field model for binary and/or ordinal variables.

Description

The function bgm explores the joint pseudoposterior distribution of parameters and possibly edge indicators for a Markov Random Field model for mixed binary and ordinal variables.

Usage

bgm(
  x,
  variable_type = "ordinal",
  reference_category,
  iter = 10000,
  burnin = 500,
  interaction_scale = 2.5,
  threshold_alpha = 0.5,
  threshold_beta = 0.5,
  edge_selection = TRUE,
  edge_prior = c("Bernoulli", "Beta-Bernoulli", "Stochastic-Block"),
  inclusion_probability = 0.5,
  beta_bernoulli_alpha = 1,
  beta_bernoulli_beta = 1,
  dirichlet_alpha = 1,
  lambda = 1,
  na_action = c("listwise", "impute"),
  save = FALSE,
  display_progress = TRUE
)

Arguments

x	A data frame or matrix with n rows and p columns containing binary and ordinal variables for n independent observations and p variables in the network. Regular binary and ordinal variables are recoded as non-negative integers (0, 1, ..., m) if not already done. Unobserved categories are collapsed into other categories after recoding (i.e., if category 1 is unobserved, the data are recoded from (0, 2) to (0, 1)). Blume-Capel ordinal variables are also coded as non-negative integers if not already done. However, since “distance” to the reference category plays an important role in this model, unobserved categories are not collapsed after recoding.
variable_type	What kind of variables are there in x? Can be a single character string specifying the variable type of all p variables at once or a vector of character strings of length p specifying the type for each variable in x separately. Currently, bgm supports “ordinal” and “blume-capel”. Binary variables are automatically treated as “ordinal’’. Defaults to variable_type = "ordinal".
reference_category	The reference category in the Blume-Capel model. Should be an integer within the range of integer scores observed for the “blume-capel” variable. Can be a single number specifying the reference category for all Blume-Capel variables at once, or a vector of length p where the i-th element contains the reference category for variable i if it is Blume-Capel, and bgm ignores its elements for other variable types. The value of the reference category is also recoded when bgm recodes the corresponding observations. Only required if there is at least one variable of type “blume-capel”.
iter	How many iterations should the Gibbs sampler run? The default of 1e4 is for illustrative purposes. For stable estimates, it is recommended to run the Gibbs sampler for at least 1e5 iterations.
burnin	The number of iterations of the Gibbs sampler before saving its output. 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. When edge_selection = TRUE, the bgm function will perform 2 * burnin iterations, first burnin iterations without edge selection, then burnin iterations with edge selection. This helps ensure that the Markov chain used for estimation starts with good parameter values and that the adaptive MH proposals are properly calibrated.
interaction_scale	The scale of the Cauchy distribution that is used as a prior for the pairwise interaction parameters. Defaults to 2.5.
threshold_alpha, threshold_beta	The shape parameters of the beta-prime prior density for the threshold parameters. Must be positive values. If the two values are equal, the prior density is symmetric about zero. If threshold_beta is greater than threshold_alpha, the distribution is skewed to the left, and if threshold_beta is less than threshold_alpha, it is skewed to the right. Smaller values tend to lead to more diffuse prior distributions.
edge_selection	Should the function perform Bayesian edge selection on the edges of the MRF in addition to estimating its parameters (edge_selection = TRUE), or should it just estimate the parameters (edge_selection = FALSE)? The default is edge_selection = TRUE.
edge_prior	The inclusion or exclusion of individual edges in the network is modeled with binary indicator variables that capture the structure of the network. The argument edge_prior is used to set a prior distribution for the edge indicator variables, i.e., the structure of the network. Currently, three options are implemented: The Bernoulli model edge_prior = "Bernoulli" assumes that the probability that an edge between two variables is included is equal to inclusion_probability and independent of other edges or variables. When inclusion_probability = 0.5, this means that each possible network structure is given the same prior weight. The Beta-Bernoulli model edge_prior = "Beta-Bernoulli" assumes a beta prior for the unknown inclusion probability with shape parameters beta_bernoulli_alpha and beta_bernoulli_beta. If beta_bernoulli_alpha = 1 and beta_bernoulli_beta = 1, this means that networks with the same complexity (number of edges) get the same prior weight. The Stochastic Block model edge_prior = "Stochastic-Block" assumes that nodes can be organized into blocks or clusters. In principle, the assignment of nodes to such clusters is unknown, and the model as implemented here considers all possible options (i.e., specifies a Dirichlet prior on the probability of allocations as described by Geng et al. 2019). This model is advantageous when nodes are expected to fall into distinct clusters. The inclusion probabilities for the edges are defined at the level of the clusters, with a beta prior for the unknown inclusion probability with shape parameters beta_bernoulli_alpha and beta_bernoulli_beta, and a Dirichlet prior on the cluster assignment probabilities with a common concentration parameter dirichlet_alpha and a zero-truncated Poisson prior on the number of clusters with a rate parameter lambda, indicating the expected number of clusters. The default is edge_prior = "Bernoulli".
inclusion_probability	The prior edge inclusion probability for the Bernoulli model. Can be a single probability, or a matrix of p rows and p columns specifying an inclusion probability for each edge pair. The default is inclusion_probability = 0.5.
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 beta_bernoulli_alpha = 1 and beta_bernoulli_beta = 1.
dirichlet_alpha	The shape of the Dirichlet prior on the node-to-block allocation probabilities for the Stochastic Block model.
lambda	The rate parameter of the zero-truncated Poisson prior on the number of cluster for the Stochastic Block model.
na_action	How do you want the function to handle missing data? If na_action = "listwise", listwise deletion is used. If na_action = "impute", missing data are imputed iteratively during the MCMC procedure. Since imputation of missing data can have a negative impact on the convergence speed of the MCMC procedure, it is recommended to run the MCMC for more iterations. Also, since the numerical routines that search for the mode of the posterior do not have an imputation option, the bgm function will automatically switch to interaction_prior = "Cauchy" and adaptive = TRUE.
save	Should the function collect and return all samples from the Gibbs sampler (save = TRUE)? Or should it only return the (model-averaged) posterior means (save = FALSE)? Defaults to FALSE.
display_progress	Should the function show a progress bar (display_progress = TRUE)? Or not (display_progress = FALSE)? The default is TRUE.

Details

Currently, bgm supports two types of ordinal variables. The regular, default, ordinal variable type has no restrictions on its distribution. Every response category except the first receives its own threshold parameter. The Blume-Capel ordinal variable assumes that there is a specific reference category, such as the “neutral” in a Likert scale, and responses are scored in terms of their distance to this reference category. Specifically, the Blume-Capel model specifies the following quadratic model for the threshold parameters:

$\mu_{\text{c}} = \alpha \times \text{c} + \beta \times (\text{c} - \text{r})^2,$

where $\mu_{\text{c}}$ is the threshold for category c. The parameter $\alpha$ models a linear trend across categories, such that $\alpha > 0$ leads to an increasing number of observations in higher response categories and $\alpha <0$ leads to a decreasing number of observations in higher response categories. The parameter $\beta$ models the response style in terms of an offset with respect to the reference category $r$; if $\beta<0$ there is a preference to respond in the reference category (i.e., the model introduces a penalty for responding in a category further away from the reference_category category r), while if $\beta > 0$ there is preference to score in the extreme categories further away from the reference_category category.

The Bayesian estimation procedure (edge_selection = FALSE) simply estimates the threshold and pairwise interaction parameters of the ordinal MRF, while the Bayesian edge selection procedure (edge_selection = TRUE) also models the probability that individual edges should be included or excluded from the model. Bayesian edge selection imposes a discrete spike and slab prior distribution on the pairwise interactions. By formulating it as a mixture of mutually singular distributions, the function can use a combination of Metropolis-Hastings and Gibbs sampling to create a Markov chain that has the joint posterior distribution as an invariant. The current option for the slab distribution is a Cauchy with an optional scaling parameter. The slab distribution is also used as the prior for the interaction parameters for Bayesian estimation. A beta-prime distribution is used for the exponent of the category parameters. For Bayesian edge selection, two prior distributions are implemented for the edge inclusion variables (i.e., the prior probability that an edge is included); the Bernoulli prior and the Beta-Bernoulli prior.

Value

If save = FALSE (the default), the result is a list of class “bgms” containing the following matrices with model-averaged quantities:

• indicator: A matrix with p rows and p columns, containing the posterior inclusion probabilities of individual edges.

• interactions: A matrix with p rows and p columns, containing model-averaged posterior means of the pairwise associations.

• thresholds: A matrix with p rows and max(m) columns, containing model-averaged category thresholds. In the case of “blume-capel” variables, the first entry is the parameter for the linear effect and the second entry is the parameter for the quadratic effect, which models the offset to the reference category.

• In the case of edge_prior = "Stochastic-Block", two additional elements are returned:

  • A vector allocations with the estimated cluster assignments of the nodes, calculated using a method proposed by (Dahl2009) and also used by (GengEtAl_2019).

  • A matrix components with the estimated posterior probability of the number of components (clusters) in the network. These probabilities are calculated based on Equation 3.7 in (miller2018mixture), which computes the conditional probability of the number of components given the number of clusters. The number of clusters is derived from the cardinality of the sampled allocations vector for each iteration of the MCMC sampler (see save = TRUE).

If save = TRUE, the result is a list of class “bgms” containing:

• indicator: 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.

• In the case of edge_prior = "Stochastic-Block" a matrix allocations with the cluster assignments of the nodes from each iteration is returned. This matrix can be used to calculate the posterior probability of the number of clusters by utilizing the summary_SBM function.

Column averages of these matrices provide the model-averaged posterior means. Except for the allocations matrix, for which the summary_SBM needs to be utilized.

In addition to the analysis results, the bgm output lists some of the arguments of its call. This is useful for post-processing the results.

References

Geng J, Bhattacharya A, Pati D (2019). “Probabilistic community detection with unknown number of communities.” Journal of the American Statistical Association, 114, 893–905. doi:10.1080/01621459.2018.1458618.

Examples

#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$indicator[lower.tri(fit$indicator)], 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
 #------------------------------------------------------------------------------|

 #For the default choice of the structure prior, the prior odds equal one:
 prior.odds = 1
 posterior.inclusion = fit$indicator[lower.tri(fit$indicator)]
 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("Log-Inclusion 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 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)

---

Bayesian Variable Selection or Bayesian Estimation for Differences in Markov Random Fields

Description

The 'bgmCompare' function estimates the pseudoposterior distribution of the parameters of a Markov Random Field (MRF) model for mixed binary and ordinal variables, as well as differences in pairwise interactions and category thresholds across groups. Groups are assumed to be 'G' independent samples.

Usage

bgmCompare(
  x,
  y,
  g,
  difference_selection = TRUE,
  main_difference_model = c("Free", "Collapse", "Constrain"),
  variable_type = "ordinal",
  reference_category,
  pairwise_difference_scale = 1,
  main_difference_scale = 1,
  pairwise_difference_prior = c("Bernoulli", "Beta-Bernoulli"),
  main_difference_prior = c("Bernoulli", "Beta-Bernoulli"),
  pairwise_difference_probability = 0.5,
  main_difference_probability = 0.5,
  pairwise_beta_bernoulli_alpha = 1,
  pairwise_beta_bernoulli_beta = 1,
  main_beta_bernoulli_alpha = 1,
  main_beta_bernoulli_beta = 1,
  interaction_scale = 2.5,
  threshold_alpha = 0.5,
  threshold_beta = 0.5,
  iter = 10000,
  burnin = 500,
  na_action = c("listwise", "impute"),
  save = FALSE,
  save_main = FALSE,
  save_pairwise = FALSE,
  save_indicator = FALSE,
  display_progress = TRUE
)

Arguments

x	Data frame or matrix with binary and ordinal responses. Regular ordinal variables should be coded as integers starting from 0. Missing categories are collapsed for regular ordinal variables but retained for Blume-Capel variables.
y	A data frame or matrix similar to 'x', used for two-group designs. 'x' contains Group 1 data, and 'y' contains Group 2 data. Ignored for multi-group designs.
g	Group membership vector for rows in 'x'. Required for multi-group designs; ignored if 'y' is provided.
difference_selection	Logical. Enables modeling of inclusion/exclusion of parameter differences ('TRUE') or estimation of all differences ('FALSE'). Default: 'TRUE'.
main_difference_model	Character. Specifies how to handle threshold differences when categories are unmatched. Options: '"Collapse"', '"Free"'. The option "Collapse" collapses categories unobserved in one or more groups. The option "Free" option estimates thresholds separately for each group and does not model their difference. Default: '"Free"'.
variable_type	Character or vector. Specifies the type of variables in 'x' ('"ordinal"' or '"blume-capel"'). Default: '"ordinal"'.
reference_category	Integer or vector. Reference category for Blume-Capel variables. Required if there is at least one Blume-Capel variable.
pairwise_difference_scale	Double. Scale parameter for the Cauchy prior on pairwise differences. Default: '1'.
main_difference_scale	Double. Scale parameter for the Cauchy prior on threshold differences. Default: '1'.
pairwise_difference_prior, main_difference_prior	Character. Specifies the inclusion probability model ('"Bernoulli"' or '"Beta-Bernoulli"'). Default: '"Bernoulli"'.
pairwise_difference_probability	A numeric value or a $p \times p$ matrix specifying the prior inclusion probability of a pairwise difference in the Bernoulli model. A single value applies the same probability to all pairs, while a matrix allows for edge-specific probabilities. Defaults to 0.5 (equal prior probability for inclusion and exclusion).
main_difference_probability	A numeric value or a length-$p$ vector specifying the prior inclusion probability of a threshold difference in the Bernoulli model. A single value applies the same probability to all variables, while a vector allows for variable-specific probabilities. Defaults to 0.5, indicating no prior preference.
pairwise_beta_bernoulli_alpha, pairwise_beta_bernoulli_beta	Double. Shape parameters for the Beta-Bernoulli prior on pairwise differences.
main_beta_bernoulli_alpha, main_beta_bernoulli_beta	Double. Shape parameters for the Beta-Bernoulli prior on threshold differences.
interaction_scale	Double. Scale of the Cauchy prior for nuisance pairwise interactions. Default: '2.5'.
threshold_alpha, threshold_beta	Double. Shape parameters for the beta-prime prior on nuisance threshold parameters.
iter, burnin	Integer. Number of Gibbs iterations ('iter') and burn-in iterations ('burnin'). Defaults: 'iter = 1e4', 'burnin = 5e2'.
na_action	Character. Specifies handling of missing data. '"listwise"' deletes rows with missing values; '"impute"' imputes values during Gibbs sampling. Default: '"listwise"'.
save	Logical. If true, sampled states for all parameters are returned. Deprecated.
save_main, save_pairwise, save_indicator	Logical. Enable saving sampled states for 'main_effects', 'pairwise_effects', and 'indicator', respectively. Default: 'FALSE'.
display_progress	Logical. Show progress bar during computation. Default: 'TRUE'.

Details

### Pairwise Interactions Pairwise interactions between variables 'i' and 'j' are modeled as:

$\boldsymbol{\theta}_{ij} = \phi_{ij} + \boldsymbol{\delta}_{ij},$

where: - $\boldsymbol{\theta}_{ij}$ is the vector of pairwise interaction parameters of length 'G'. - $\phi_{ij}$ is the overall pairwise interaction (nuisance parameter). - $\boldsymbol{\delta}_{ij}$ represents group-specific differences constrained to sum to zero for identification.

### Ordinal Variables The function supports two types of ordinal variables: 1. **Regular ordinal variables**: Introduce a threshold parameter for each category except the lowest, modeled as:

$\boldsymbol{\mu}_{ic} = \tau_{ic} + \boldsymbol{\epsilon}_{ic},$

where: - $\tau_{ic}$ denotes an overall effect (nuisance parameter). - $\boldsymbol{\epsilon}_{ic}$ represents group-specific differences constrained to sum to zero.

2. **Blume-Capel ordinal variables**: Assume a specific reference category and score responses based on distance to it:

$\boldsymbol{\mu}_{ic} = (\tau_{i1} + \boldsymbol{\epsilon}_{i1}) \cdot c + (\tau_{i2} + \boldsymbol{\epsilon}_{i2}) \cdot (c - r)^2,$

where: - 'r' is the reference category. - $\tau_{i1}$ and $\tau_{i2}$ are nuisance parameters. - $\boldsymbol{\epsilon}_{i1}$ and $\boldsymbol{\epsilon}_{i2}$ represent group-specific differences.

### Variable Selection Bayesian variable selection enables testing of parameter differences or equivalence across groups. Independent spike-and-slab priors are applied to difference parameters: - **Bernoulli Model**: Assigns a fixed probability to parameter inclusion. - **Beta-Bernoulli Model**: Incorporates a beta prior to model inclusion probabilities.

### Saving Options Users can store sampled states for parameters ('main_effects', 'pairwise_effects', 'indicator') during Gibbs sampling. Enabling these options ('save_main', 'save_pairwise', 'save_indicator') increases output size and memory usage, so use them judiciously.

Value

A list containing the posterior means and, optionally, sampled states based on the 'save_*' options. The returned components include: - 'posterior_mean_main', 'posterior_mean_pairwise', and 'posterior_mean_indicator' for posterior means. - If saving options are enabled: - 'raw_samples_main' for sampled states of 'main_effects'. - 'raw_samples_pairwise' for sampled states of 'pairwise_effects'. - 'raw_samples_indicator' for sampled states of the inclusion indicators.

In addition to the results of the analysis, the output lists some of the arguments of its call. This is useful for post-processing the results.

---

Sample observations from the ordinal MRF

Description

This function samples states from the ordinal MRF using a Gibbs sampler. The Gibbs sampler is initiated with random values from the response options, after which it proceeds by simulating states for each variable from a logistic model using the other variable states as predictor variables.

Usage

mrfSampler(
  no_states,
  no_variables,
  no_categories,
  interactions,
  thresholds,
  variable_type = "ordinal",
  reference_category,
  iter = 1000
)

Arguments

no_states	The number of states of the ordinal MRF to be generated.
no_variables	The number of variables in the ordinal MRF.
no_categories	Either a positive integer or a vector of positive integers of length no_variables. The number of response categories on top of the base category: no_categories = 1 generates binary states.
interactions	A symmetric no_variables by no_variables matrix of pairwise interactions. Only its off-diagonal elements are used.
thresholds	A no_variables by max(no_categories) matrix of category thresholds. The elements in row i indicate the thresholds of variable i. If no_categories is a vector, only the first no_categories[i] elements are used in row i. If the Blume-Capel model is used for the category thresholds for variable i, then row i requires two values (details below); the first is $\alpha$, the linear contribution of the Blume-Capel model and the second is $\beta$, the quadratic contribution.
variable_type	What kind of variables are simulated? Can be a single character string specifying the variable type of all p variables at once or a vector of character strings of length p specifying the type for each variable separately. Currently, bgm supports “ordinal” and “blume-capel”. Binary variables are automatically treated as “ordinal’’. Defaults to variable_type = "ordinal".
reference_category	An integer vector of length no_variables specifying the reference_category category that is used for the Blume-Capel model (details below). Can be any integer value between 0 and no_categories (or no_categories[i]).
iter	The number of iterations used by the Gibbs sampler. The function provides the last state of the Gibbs sampler as output. By default set to 1e3.

Details

There are two modeling options for the category thresholds. The default option assumes that the category thresholds are free, except that the first threshold is set to zero for identification. The user then only needs to specify the thresholds for the remaining response categories. This option is useful for any type of ordinal variable and gives the user the most freedom in specifying their model.

The Blume-Capel option is specifically designed for ordinal variables that have a special type of reference_category category, such as the neutral category in a Likert scale. The Blume-Capel model specifies the following quadratic model for the threshold parameters:

$\mu_{\text{c}} = \alpha \times \text{c} + \beta \times (\text{c} - \text{r})^2,$

where $\mu_{\text{c}}$ is the threshold for category c (which now includes zero), $\alpha$ offers a linear trend across categories (increasing threshold values if $\alpha > 0$ and decreasing threshold values if $\alpha <0$), if $\beta < 0$, it offers an increasing penalty for responding in a category further away from the reference_category category r, while $\beta > 0$ suggests a preference for responding in the reference_category category.

Value

A no_states by no_variables matrix of simulated states of the ordinal MRF.

Examples

# Generate responses from a network of five binary and ordinal variables.
no_variables = 5
no_categories = sample(1:5, size = no_variables, replace = TRUE)

Interactions = matrix(0, nrow = no_variables, ncol = no_variables)
Interactions[2, 1] = Interactions[4, 1] = Interactions[3, 2] =
  Interactions[5, 2] = Interactions[5, 4] = .25
Interactions = Interactions + t(Interactions)
Thresholds = matrix(0, nrow = no_variables, ncol = max(no_categories))

x = mrfSampler(no_states = 1e3,
               no_variables = no_variables,
               no_categories = no_categories,
               interactions = Interactions,
               thresholds = Thresholds)

# Generate responses from a network of 2 ordinal and 3 Blume-Capel variables.
no_variables = 5
no_categories = 4

Interactions = matrix(0, nrow = no_variables, ncol = no_variables)
Interactions[2, 1] = Interactions[4, 1] = Interactions[3, 2] =
  Interactions[5, 2] = Interactions[5, 4] = .25
Interactions = Interactions + t(Interactions)

Thresholds = matrix(NA, no_variables, no_categories)
Thresholds[, 1] = -1
Thresholds[, 2] = -1
Thresholds[3, ] = sort(-abs(rnorm(4)), decreasing = TRUE)
Thresholds[5, ] = sort(-abs(rnorm(4)), decreasing = TRUE)

x = mrfSampler(no_states = 1e3,
               no_variables = no_variables,
               no_categories = no_categories,
               interactions = Interactions,
               thresholds = Thresholds,
               variable_type = c("b","b","o","b","o"),
               reference_category = 2)

---

Print method for bgms objects

Description

Used to prevent bgms output cluttering the console.

Usage

## S3 method for class 'bgmCompare'
print(x, ...)

Arguments

x	An object of class bgms.
...	Ignored.

---

Print method for bgms objects

Description

Used to prevent bgms output cluttering the console.

Usage

## S3 method for class 'bgms'
print(x, ...)

Arguments

x	An object of class bgms.
...	Ignored.

---

Function for summarizing the sampled cluster allocation vectors

Description

Th summarySBM function summarizes the sampled allocation vectors from each iteration of the Gibbs sampler from the output of the bgm function ran with edge_prior = "Stochastic-Block" and save = TRUE. It also estimates the posterior distribution of the number of clusters.

Usage

summarySBM(bgm_object, internal_call = FALSE)

Arguments

bgm_object	A fit object created by the bgm function.
internal_call	A logical value indicating whether the function is used within bgms for calculating the posterior probabilities of the number of clusters or by the user. This argument is always set to FALSE.

Value

Returns a list of two elements: components and allocations, containing the posterior probabilities for the number of components (clusters) and the estimated cluster allocation of the nodes using Dahl's method.

Examples

# fit a model with the SBM prior
  bgm_object = bgm(
    Wenchuan[, c(1:5)],
    edge_prior = "Stochastic-Block",
    save = TRUE)

  summarySBM(bgm_object)

---

Post-traumatic stress disorder symptoms of Wenchuan earthquake survivors

Description

A data set containing items measuring symptoms of posttraumatic stress disorder (PTSD) (McNally et al. 2015). Participants were 362 Chinese adults who survived the Wenchuan earthquake and lost at least one child in the disaster. PTSD symptoms were reported using the civilian version of the Posttraumatic Checklist, which consists of 17 items, each assessing one of the DSM-IV symptoms of PTSD. Participants rated each item on a five-point scale ranging from “not at all” to “extremely” to indicate how much the symptom bothered them in the past month.

Usage

data("Wenchuan")

Format

A matrix with 362 rows and 17 columns:

intrusion

Repeated, disturbing memories, thoughts, or images of a stressful experience from the past?

dreams

Repeated, disturbing dreams of a stressful experience from the past?

flash

Suddenly acting or feeling as if a stressful experience were happening again (as if you were reliving it)?

upset

Feeling very upset when something reminded you of a stressful experience from the past?

physior

Having physical reactions (e.g., heart pounding, trouble breathing, sweating) when something reminded you of a stressful experience from the past?

avoidth

Avoiding thinking about or talking about a stressful experience from the past or avoiding having feelings related to it?

avoidact

Avoiding activities or situations because they reminded you of a stressful experience from the past?

amnesia

Trouble remembering important parts of a stressful experience from the past?

lossint

Loss of interest in activities that you used to enjoy?

distant

Feeling distant or cut off from other people?

numb

Feeling emotionally numb or being unable to have loving feelings for those close to you?

future

Feeling as if your future will somehow be cut short?

sleep

Trouble falling or staying asleep?

anger

Feeling irritable or having angry outbursts?

concen

Having difficulty concentrating?

hyper

Being "super-alert" or watchful or on guard?

startle

Feeling jumpy or easily startled?

Source

http://psychosystems.org/wp-content/uploads/2014/10/Wenchuan.csv

References

McNally RJ, Robinaugh DJ, Wu GWY, Wang L, Deserno MK, Borsboom D (2015). “Mental disorders as causal systems: A network approach to posttraumatic stress disorder.” Clinical Psychological Science, 6, 836–849. doi:10.1177/2167702614553230.

---