Type: | Package |
Title: | Bayesian Analysis of Networks of Binary and/or Ordinal Variables |
Version: | 0.1.4.2 |
Date: | 2024-12-03 |
Maintainer: | Maarten Marsman <m.marsman@uva.nl> |
Description: | Bayesian variable selection methods for analyzing the structure of a Markov Random Field model for a network of binary and/or ordinal variables. Details of the implemented methods can be found in: Marsman, van den Bergh, and Haslbeck (in press) <doi:10.31234/osf.io/ukwrf>. |
License: | GPL-2 | GPL-3 [expanded from: GPL (≥ 2)] |
URL: | https://maartenmarsman.github.io/bgms/ |
BugReports: | https://github.com/MaartenMarsman/bgms/issues |
Imports: | Rcpp (≥ 1.0.7), Rdpack, methods |
RdMacros: | Rdpack |
LinkingTo: | Rcpp, RcppProgress |
RoxygenNote: | 7.3.2 |
Depends: | R (≥ 2.10) |
LazyData: | true |
Encoding: | UTF-8 |
Suggests: | knitr, qgraph, rmarkdown, testthat (≥ 3.0.0) |
VignetteBuilder: | knitr |
Config/testthat/edition: | 3 |
Config/Needs/website: | tidyverse/tidytemplate |
NeedsCompilation: | yes |
Packaged: | 2024-12-05 09:44:17 UTC; maartenmarsman |
Author: | Maarten Marsman |
Repository: | CRAN |
Date/Publication: | 2024-12-05 10:20:06 UTC |
bgms: Bayesian Analysis of Networks of Binary and/or Ordinal Variables
Description
The R
package bgms provides tools for Bayesian analysis of
graphical models describing networks of variables. The package uses Markov
chain Monte Carlo methods combined with a pseudolikelihood approach to
estimate the posterior distribution of model parameters.
Gibbs variable selection (George and McCulloch 1993) is used to model the underlying network structure of the graphical model. By imposing a discrete spike and slab prior on the pairwise interactions, it is possible to shrink the interactions to exactly zero. The Gibbs sampler embeds a Metropolis approach for mixtures of mutually singular distributions (Gottardo and Raftery 2008) to account for the discontinuity at zero. The goal is to provide these tools for Markov Random Field (MRF) models for a wide range of variable types in the bgms package, and it currently provides them for analyzing networks of binary and/or ordinal variables (Marsman et al. in press).
While the goal is to provide the above tools for Markov Random Field (MRF) models for a wide range of variable types in the bgms package, it currently provides tools for analyzing networks of binary and/or ordinal variables (Marsman et al. in press).
MRFs are a special class of graphical models whose graph structure reflects the conditional associations between their variables, making them useful for testing for conditional independence or dependence. For example, the inclusion Bayes factor tests for conditional independence or dependence of a pair of variables in the network by comparing the predictive adequacy of models that include the edge between these variables and models that exclude the edge. (Huth et al. 2023; Sekulovski et al. in press).
The bgms package offers several tools for analyzing the structure of the MRF:
Simulate response data from the MRF using the Gibbs sampler.
Simulate
mrfSampler
.
Estimate the posterior distribution of the MRF's parameters and possibly its network structure using Gibbs variable selection.
Bayesian estimation or Bayesian edge selection with
bgm
.
Author(s)
Maintainer: Maarten Marsman m.marsman@uva.nl (ORCID)
Other contributors:
Karoline Huth (ORCID) [contributor]
Nikola Sekulovski (ORCID) [contributor]
Don van den Bergh (ORCID) [contributor]
References
George EI, McCulloch RE (1993).
“Variable selection via Gibbs sampling.”
Journal of the American Statistical Association, 88, 881–889.
doi:10.1080/01621459.1993.10476353.
Gottardo R, Raftery AE (2008).
“Markov Chain Monte Carlo With Mixtures of Mutually Singular Distributions.”
Journal of Computational and Graphical Statistics, 17, 949–975.
doi:10.1198/106186008X386102.
Huth K, de Ron J, Goudriaan AE, Luigjes K, Mohammadi R, van Holst RJ, Wagenmakers E, Marsman M (2023).
“Bayesian analysis of cross-sectional networks: A tutorial in R and JASP.”
Advances in Methods and Practices in Psychological Science, 6, 1–18doi.
Marsman M, van den Bergh D, Haslbeck JMB (in press).
“Bayesian analysis of the ordinal Markov random field.”
Psychometrika.
Sekulovski N, Keetelaar S, Huth K, Wagenmakers E, van Bork R, van den Bergh D, Marsman M (in press).
“Testing conditional independence in psychometric networks: An analysis of three bayesian methods.”
Multivariate Behavioral Research.
doi:10.1080/00273171.2024.2345915.
See Also
Useful links:
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.
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,
na.action = c("listwise", "impute"),
save = FALSE,
display_progress = TRUE
)
Arguments
x |
A data frame or matrix with |
variable_type |
What kind of variables are there in |
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
|
iter |
How many iterations should the Gibbs sampler run? The default of
|
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 |
interaction_scale |
The scale of the Cauchy distribution that is used as
a prior for the pairwise interaction parameters. Defaults to |
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
|
edge_selection |
Should the function perform Bayesian edge selection on
the edges of the MRF in addition to estimating its parameters
( |
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 |
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 |
dirichlet_alpha |
The shape of the Dirichlet prior on the node-to-block allocation parameters for the Stochastic Block model. |
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
( |
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 withp
rows andp
columns, containing the posterior inclusion probabilities of individual edges. -
interactions
: A matrix withp
rows andp
columns, containing model-averaged posterior means of the pairwise associations. -
thresholds
: A matrix withp
rows andmax(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.
If save = TRUE
, the result is a list of class “bgms” containing:
-
indicator
: A matrix withiter
rows andp * (p - 1) / 2
columns, containing the edge inclusion indicators from every iteration of the Gibbs sampler. -
interactions
: A matrix withiter
rows andp * (p - 1) / 2
columns, containing parameter states from every iteration of the Gibbs sampler for the pairwise associations. -
thresholds
: A matrix withiter
rows andsum(m)
columns, containing parameter states from every iteration of the Gibbs sampler for the category thresholds.
Column averages of these matrices provide the model-averaged posterior means.
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 the Markov random field model for binary and/or ordinal variables in two independent samples.
Description
The bgmCompare
function estimates the pseudoposterior distribution of
the parameters of a Markov Random Field model for mixed binary and ordinal
variables, and the differences in pairwise interactions and category thresholds
between two groups. The groups are assumed to be two independent samples.
Usage
bgmCompare(
x,
y,
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,
display_progress = TRUE
)
Arguments
x |
A data frame or matrix with |
y |
A data frame or matrix with |
difference_selection |
Logical. If |
main_difference_model |
A string specifying options for how the bgmCompare function should handle the comparison of threshold parameters when the observed categories in the samples do not match. The "Collapse" option tells bgmCompare to collapse the two categories into one (for the data set where both categories were observed). The "Constrain" option sets the difference between the category thresholds in the two data sets to zero if the category is not observed in one of the two data sets. The "Free" option tells bgmCompare to estimate a separate set of thresholds in the two samples and to not model their differences. |
variable_type |
A string or vector specifying the type of variables
in |
reference_category |
The reference category in the Blume-Capel model.
Should be an integer within the range of integer values observed for the "blume-capel"
variable. Can be a single number that sets the reference category for all
Blume-Capel variables at once, or a vector of length |
pairwise_difference_scale |
The scale of the Cauchy distribution that is
used as the prior for the pairwise difference parameters. Defaults to
|
main_difference_scale |
The scale of the Cauchy distribution that
is used as the prior for the threshold difference parameters. Defaults to
|
pairwise_difference_prior |
A character string that specifies the model to use for the inclusion probability of pairwise differences. Options are "Bernoulli" or "Beta-Bernoulli". Default is "Bernoulli". |
main_difference_prior |
A character string that specifies the model to use for the inclusion probability of threshold differences. Options are "Bernoulli" or "Beta-Bernoulli". Default is "Bernoulli". |
pairwise_difference_probability |
The inclusion probability for a
pairwise difference in the Bernoulli model. Can be a single probability or a
matrix of |
main_difference_probability |
The inclusion probability for a
threshold difference in the Bernoulli model. Defaults to 0.5, implying no
prior preference. Can be a single probability or a vector of length |
pairwise_beta_bernoulli_alpha |
The alpha parameter of the beta distribution for the Beta-Bernoulli model for group differences in pairwise interactions. Default is 1. |
pairwise_beta_bernoulli_beta |
The beta parameter of the beta distribution for the Beta-Bernoulli model for group differences in pairwise interactions. Default is 1. |
main_beta_bernoulli_alpha |
The alpha parameter of the beta distribution for the Beta-Bernoulli model for group differences in category thresholds. Default is 1. |
main_beta_bernoulli_beta |
The beta parameter of the beta distribution for the Beta-Bernoulli model for group differences in category thresholds. Default is 1. |
interaction_scale |
The scale of the Cauchy distribution that is used as
a prior for the nuisance pairwise interaction parameters. Defaults to
|
threshold_alpha , threshold_beta |
The shape parameters of the beta-prime
prior density for the nuisance threshold parameters. Must be positive values.
If the two values are equal, the prior density is symmetric about zero. If
|
iter |
The function uses a Gibbs sampler to sample from the posterior
distribution of the model parameters and indicator variables. How many
iterations should this Gibbs sampler run? The default of |
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 |
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
( |
Details
In the first group, the pairwise interactions between the variables i
and j
are modeled as
\sigma_{\text{ij}} = \theta_{\text{ij}} + \delta_{\text{ij}} / 2,
and in the second group as
\sigma_{\text{ij}} = \theta_{\text{ij}} - \delta_{\text{ij}} / 2,
The pairwise interaction parameter \theta_{\text{ij}}
denotes an overall effect that is considered nuisance, and attention is focused
on the pairwise difference parameter \delta_{\text{ij}}
,
which reflects the difference in the pairwise interaction between the two groups.
The bgmCompare
function supports two types of ordinal variables, which
can be mixed. The default ordinal variable introduces a threshold parameter
for each category except the lowest category. For this variable type, the threshold parameter for
variable i
, category c
, is modeled as
\mu_{\text{ic}} = \tau_{\text{ic}} + \epsilon_{\text{ic}} / 2,
in the first group and in the second group as
\mu_{\text{ic}} = \tau_{\text{ic}} - \epsilon_{\text{ic}} / 2,
The category threshold parameter \tau_{\text{ic}}
denotes
an overall effect that is considered nuisance, and attention is focused on the
threshold difference parameter \epsilon_{\text{ic}}
,
which reflects the difference in threshold of for variable i
, category
c
between the two groups.
The Blume-Capel ordinal variable assumes that there is a specific reference category, such as “neutral” in a Likert scale, and responses are scored according to their distance from this reference category. In the first group, the threshold parameters are modelled as
\mu_{\text{ic}} = (\tau_{\text{i1}} + \epsilon_{\text{i1}} / 2) \times \text{c} + (\tau_{\text{i2}} + \epsilon_{\text{i2}} / 2) \times (\text{c} - \text{r})^2,
and in the second groups as
\mu_{\text{ic}} = (\tau_{\text{i1}} - \epsilon_{\text{i1}} / 2) \times \text{c} + (\tau_{\text{i2}} - \epsilon_{\text{i2}} / 2) \times (\text{c} - \text{r})^2.
The linear and quadratic category threshold parameters
\tau_{\text{i1}}
and \tau_{\text{i2}}
denote overall effects that are considered nuisance, and attention is focused
on the two threshold difference parameters
\epsilon_{\text{i1}}
and
\epsilon_{\text{i2}}
, which reflect the differences
in the quadratic model for the variable i
between the two groups.
Bayesian variable selection is used to model the presence or absence of the
difference parameters \delta
and \epsilon
, which
allow us to assess parameter differences between the two groups. Independent
spike and slab priors are specified for these difference parameters. The spike
and slab priors use binary indicator variables to select the difference parameters,
assigning them a diffuse Cauchy prior with an optional scaling parameter if
selected, or setting the difference parameter to zero if not selected.
The function offers two models for the probabilistic inclusion of parameter differences:
-
Bernoulli Model: This model assigns a fixed probability of selecting a parameter difference, treating them as independent events. A probability of 0.5 indicates no preference, giving equal prior weight to all configurations.
-
Beta-Bernoulli Model: Introduces a beta distribution prior for the inclusion probability that models the complexity of the configuration of the difference indicators. When the alpha and beta shape parameters of the beta distribution are set to 1, the model assigns the same prior weight to the number of differences present (i.e., a configuration with two differences or with four differences is a priori equally likely).
Inclusion probabilities can be specified for pairwise interactions with
pairwise_difference_probability
and for category thresholds with
threshold_difference_probability
.
The pairwise interaction parameters \theta
, the category
threshold parameters \tau
, and, in the not yet implemented paired-samples designs,
the between-sample interactions \omega
are considered
nuisance parameters that are common to all models. The pairwise interaction
parameters \theta
and the between-sample interactions
\omega
are assigned a diffuse Cauchy prior with an optional
scaling parameter. The exponent of the category threshold parameters
\tau
are assigned beta-prime distribution with optional scale
values.
Value
If save = FALSE
(the default), the result is a list of class
“bgmCompare” containing the following matrices:
-
indicator
: A matrix withp
rows andp
columns containing the posterior inclusion probabilities of the differences in pairwise interactions on the off-diagonal and the posterior inclusion probabilities of the differences in category thresholds on the diagonal. -
difference_pairwise
: A matrix withp
rows andp
columns, containing model-averaged posterior means of the differences in pairwise interactions. -
difference_threshold
: A matrix withp
rows andmax(m)
columns, containing model-averaged posterior means of the differences in category thresholds. -
interactions
: A matrix withp
rows andp
columns, containing posterior means of the nuisance pairwise interactions. -
thresholds
: A matrix withp
rows andmax(m)
columns containing the posterior means of the nuisance 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.
If save = TRUE
, the result is a list of class “bgmCompare”
containing the following matrices:
-
indicator_pairwise
: A matrix withiter
rows andp * (p - 1) / 2
columns containing the inclusion indicators for the differences in pairwise interactions from each iteration of the Gibbs sampler. -
difference_pairwise
: A matrix withiter
rows andp * (p - 1) / 2
columns, containing parameter states for the differences in pairwise interactions from each iteration of the Gibbs sampler. -
indicator_threshold
: A matrix withiter
rows andsum(m)
columns, containing the inclusion indicators for the differences in category thresholds from each iteration of the Gibbs sampler. -
difference_threshold
: A matrix withiter
rows andsum(m)
columns, containing the parameter states for the differences in category thresholds from each iteration of the Gibbs sampler. -
interactions
: A matrix withiter
rows andp * (p - 1) / 2
columns, containing parameter states for the nuisance pairwise interactions in each iteration of the Gibbs sampler. -
thresholds
: A matrix withiter
rows andsum(m)
columns, containing parameter states for the nuisance category thresholds in each iteration of the Gibbs sampler.
Column averages of these matrices provide the model-averaged posterior means.
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.
Extractor Functions.
Description
Extractor Functions.
Usage
extract_arguments(bgms_object)
## S3 method for class 'bgms'
extract_arguments(bgms_object)
## S3 method for class 'bgmCompare'
extract_arguments(bgms_object)
extract_indicators(bgms_object)
## S3 method for class 'bgms'
extract_indicators(bgms_object)
## S3 method for class 'bgmCompare'
extract_indicators(bgms_object)
extract_posterior_inclusion_probabilities(bgms_object)
## S3 method for class 'bgms'
extract_posterior_inclusion_probabilities(bgms_object)
## S3 method for class 'bgmCompare'
extract_posterior_inclusion_probabilities(bgms_object)
extract_indicator_priors(bgms_object)
## S3 method for class 'bgms'
extract_indicator_priors(bgms_object)
## S3 method for class 'bgmCompare'
extract_indicator_priors(bgms_object)
extract_pairwise_interactions(bgms_object)
extract_category_thresholds(bgms_object)
extract_pairwise_difference.bgmCompare(bgms_object)
extract_main_difference.bgmCompare(bgms_object)
extract_edge_indicators(bgms_object)
extract_pairwise_thresholds(bgms_object)
Arguments
bgms_object |
A fit object created by the bgms package or specifically by the bgm function. |
Details
Extract results from bgm objects in a safe way. Mainly intended for developers of packages that build on top of the bgms package.
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 |
interactions |
A symmetric |
thresholds |
A |
variable_type |
What kind of variables are simulated? Can be a single
character string specifying the variable type of all |
reference_category |
An integer vector of length |
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 |
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 |
... |
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 |
... |
Ignored. |