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

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`.

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.

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.