2.1 A psychometric approach to the ordinal Markov random field

Anderson and Vermunt (Reference Anderson and Vermunt2000) established the MRF in the psychometric literature as the distribution of manifest ordinal variables in contexts containing continuous latent variables. Such relationships between manifest and latent variables are central to psychometrics (e.g., Holland & Rosenbaum, Reference Holland and Rosenbaum1986), and the psychometric literature reports various approaches to establishing these relationships. One well-known method is the Dutch identity of Holland (Reference Holland1990), which uses the same underlying assumptions as the approach of Anderson and Vermunt (Reference Anderson and Vermunt2000; see Anderson & Yu, Reference Anderson and Yu2007, for a discussion). More recently, Marsman et al. (Reference Marsman, Sigurdardóttir, Bolsinova and Maris2019) proposed a related method that emerged from the connection between a binary variable MRF and item response theory (IRT) models in network psychometrics (e.g., Epskamp et al., Reference Epskamp, Maris, Waldorp, Borsboom, Irwing, Hughes and Booth2018a; Marsman et al., Reference Marsman, Borsboom, Kruis, Epskamp, van Bork, Waldorp and Maris2018; Marsman et al., Reference Marsman, Maris, Bechger and Glas2015). See Hessen (Reference Hessen2020) for yet another method. These approaches make the same assumptions about the distribution of unobserved latent variables and thus can all be used to obtain the results of Anderson and Vermunt (Reference Anderson and Vermunt2000). Here, we use the general setup of Marsman et al. (Reference Marsman, Sigurdardóttir, Bolsinova and Maris2019) because it provides a straightforward strategy for deriving the MRF from an existing IRT model.

Let $X_i$ denote an ordinal variable with $m+1$ response categories and realizations $x_i = 0\text {, }1\text {, }\dots \text {, }m$ , for $i = 1\text {, }2\text {, }\dots \text {, }p$ variables.Footnote ² We assume that the joint distribution for the full vector of p response variables $\mathbf {X}$ can be modeled using an item response theory (IRT) model, $P(\mathbf {X} \mid \boldsymbol {\theta })$ , in which a vector of continuous latent variables $\boldsymbol {\theta }$ of at most dimension $p-1$ fully accounts for the dependencies between the p response variables. In this context, Muraki (Reference Muraki1990) proposed the generalized partial credit model (GPCM) for ordinal responses. The multidimensional version of this GPCM is characterized by the following probability of observing category h of the ordinal response variable:

$$\begin{align*}P(X_i = h \mid \boldsymbol{\theta}) = \frac{\exp\left(h\,\boldsymbol{\alpha}_i^{\textsf T}\boldsymbol{\theta} + \beta_{ih}\right)}{1 + \sum_{u=1}^{m}\exp\left(u\,\boldsymbol{\alpha}_i^{\textsf T}\boldsymbol{\theta} + \beta_{iu}\right)}, \end{align*}$$

where $\beta _{ih}$ denotes a threshold parameter for category h of variable or item i, with $\beta _{i0}$ set to zero for identification, and $\boldsymbol {\alpha }_i$ denotes a vector of factor-loadings that relate the latent variables to the response probabilities for item i. If we let $f(\boldsymbol {\theta })$ denote the distribution for the latent variables, the full statistical model has the form

$$\begin{align*}P(\mathbf{X}) = P(X_1\text{, }\dots\text{, }X_p) = \int_{\mathbb{R}^{p-1}} \prod_{i=1}^p P(X_i\mid \boldsymbol{\theta})\, f(\boldsymbol{\theta})\,\textbf{d}\boldsymbol{\theta}. \end{align*}$$

We aim to find an expression for $P(\mathbf {X})$ , the marginal distribution of the response variables.

Marsman et al. (Reference Marsman, Sigurdardóttir, Bolsinova and Maris2019) showed that if an IRT model is in an exponential family form, we can define a latent variable distribution that allows us to analytically express the marginal distribution $P(\mathbf {X})$ . We define the exponential family form of the IRT model as follows:

$$\begin{align*}\prod_{i=1}^p P(X_i\mid \boldsymbol{\theta}) = \prod_{i=1}^p\frac{1}{\text{z}_i(\boldsymbol{\theta})}\, \text{b}_i(X_i)\,\exp\left(\textbf{S}_i(X_i)^{\textsf T}\boldsymbol{\theta}\right), \end{align*}$$

where $\textbf {S}_i(X_i)$ is a (possibly vector-valued) statistic that is sufficient for the (possibly vector-valued) latent variable $\boldsymbol {\theta }$ , $\text {b}_i(X_i)$ is a base measure that does not depend on the latent variable, and $\text {z}_i(\boldsymbol {\theta })$ is a normalizing constant that does not depend on the observed response variable and ensures that the probabilities add up to one. The normalizing constant is equal to

$$\begin{align*}\text{z}_i(\boldsymbol{\theta}) = \sum_{X_i = 0}^{m}\text{b}_i(X_i)\,\exp\left(\textbf{s}_i(X_i)^{\textsf T}\boldsymbol{\theta}\right). \end{align*}$$

We can express the multidimensional GPCM in this form. That is, we define the sufficient statistic as $\textbf {S}_i(X_i) = X_i\,\boldsymbol {\alpha }_i$ , the logarithm of the base measure as $\ln (\text {b}_i(X_i)) = \beta _{i\,X_i}$ , and the normalizing constant as

$$\begin{align*}\text{z}_i(\boldsymbol{\theta}) = 1 + \sum_{u=1}^{m}\exp\left(u\,\boldsymbol{\alpha}_i^{\textsf T}\boldsymbol{\theta} + \beta_{iu}\right). \end{align*}$$

We can now define a latent variable distribution as follows:

$$\begin{align*}f(\boldsymbol{\theta}) = \frac{1}{\text{Z}} \, \prod_{i=1}^p \text{z}_i(\boldsymbol{\theta})\, k(\boldsymbol{\theta}), \end{align*}$$

where $k(\boldsymbol {\theta })$ is a density function and $\text {Z}$ is a normalizing constant. The normalizing constant is equal to

(1) $$ \begin{align} \text{Z} = \int_{\mathbb{R}^{p-1}} \prod_{i=1}^p \text{z}_i(\boldsymbol{\theta})\, k(\boldsymbol{\theta})\,\textbf{d}\boldsymbol{\theta}. \end{align} $$

The proposed distribution $f(\boldsymbol {\theta })$ is valid for any density function $k(\boldsymbol {\theta })$ for which $0 < \text {Z} < \infty $ . Corollary 1 in Marsman et al. (Reference Marsman, Sigurdardóttir, Bolsinova and Maris2019) shows that if we assume that the density $k(\boldsymbol {\theta })$ is a multivariate normal distribution with a mean vector $\textbf {m} = \boldsymbol {0}$ and covariance matrix $\mathbf {V} = \mathbf {I}$ , the identity matrix, we find the following expression for the marginal distribution:

$$\begin{align*}P(\mathbf{X}\mid \boldsymbol{\alpha}\text{, }\boldsymbol{\beta}) = \frac{1}{\text{Z}} \exp\left(\sum_{i=1}^p\beta_{i\, X_i} + \sum_{i=1}^pX_i^2\boldsymbol{\alpha}_i^{\textsf T}\boldsymbol{\alpha}_i + 2\sum_{i=1}^{p-1}\sum_{j=i+1}^pX_iX_j\boldsymbol{\alpha}_i^{\textsf T}\boldsymbol{\alpha}_j\right). \end{align*}$$

From here, we recognize the quadratic form

(2) $$ \begin{align} P(\mathbf{X} \mid \boldsymbol{\mu}\text{, }\boldsymbol{\Sigma}) = \frac{1}{\text{Z}} \exp\left(\sum_{i=1}^p \sum_{h=1}^{m}\mathbb{I}(X_i = h)\,\mu_{ih} + 2\sum_{i=1}^{p-1}\sum_{j=i+1}^pX_iX_j\sigma_{ij}\right), \end{align} $$

where $\mu _{ih} = \beta _{ih} + h^2 \boldsymbol {\alpha }_i^{\textsf T}\boldsymbol {\alpha }_i$ and $\boldsymbol {\Sigma } = [\sigma _{ij}] = \boldsymbol {\alpha }_i^{\textsf T}\boldsymbol {\alpha }_j$ , and $\sigma _{ii} = 0$ , for $i=1\text {, }\dots \text {, }p$ . The marginal distribution $P(\mathbf {X} \mid \boldsymbol {\mu }\text {, }\boldsymbol {\Sigma })$ is the desired MRF and equal to Eq. (20) in Anderson and Vermunt (Reference Anderson and Vermunt2000; see also Theorem 2 in Hessen, Reference Hessen2020). The normalizing constant that we defined in Eq. (1) can now be restated as

(3) $$ \begin{align} \text{Z} = \sum_{\mathbf{X} \in {\Omega}_X} \exp\left(\sum_{i=1}^p \sum_{h=1}^{m}\mathbb{I}(X_i = h)\,\mu_{ih} + 2\sum_{i=1}^{p-1}\sum_{j=i+1}^pX_iX_j\sigma_{ij}\right), \end{align} $$

where ${\Omega}_{X} = \{0\text {, }\dots \text {, }m\}^{p}$ denotes the space of possible response patterns. Observe that when $m = 1$ , the MRF in Eq. (2) simplifies to the Ising model (Cox, Reference Cox1972; Ising, Reference Ising1925).

2.2 A Hammersley–Clifford approach to the ordinal Markov random field

Suggala et al. (Reference Suggala, Yang, Ravikumar, Precup and Teh2017) established the MRF in the statistical learning literature as a unique multivariate distribution, consistent with a specific full conditional distribution for ordinal variables. They extended the work of Yang et al. (Reference Yang, Allen, Liu, Ravikumar, Pereira, Burges, Bottou and Weinberger2012), who provided a method for constructing a multivariate graphical model from univariate full conditional distributions. These approaches follow the Hammersley–Clifford approach advocated by Besag (Reference Besag1974). Suggala et al. (Reference Suggala, Yang, Ravikumar, Precup and Teh2017) set out to use the method proposed by Yang et al. (Reference Yang, Allen, Liu, Ravikumar, Pereira, Burges, Bottou and Weinberger2012) to establish the multivariate distribution consistent with standard univariate models for ordinal variables based on cumulative, continuous, and consecutive ratio logits (Agresti, Reference Agresti2010, Reference Agresti2018). However, the method of Yang et al. (Reference Yang, Allen, Liu, Ravikumar, Pereira, Burges, Bottou and Weinberger2012) assumes that the conditional distributions are in the exponential family. Suggala et al. (Reference Suggala, Yang, Ravikumar, Precup and Teh2017) show that the cumulative and continuation-based logits do not belong to this class of models and consequently do not yield a multivariate distribution. On the other hand, the consecutive-based logit models belong to the exponential family and are also consistent with a multivariate distribution. These results are discussed next.

The Hammersley–Clifford theorem states that if a multivariate distribution $P(X_1\text {, }\dots \text {, }X_p)$ is consistent with certain full conditional distributions $P(X_i \mid X_1\text {, }\dots \text {, }X_{i-1}\text {, }X_{i+1}\text {, }\dots \text {, }X_p) = P(X_i \mid \mathbf {X}^{(i)}),$ then it is the only multivariate distribution consistent with these conditional distributions. Besag (Reference Besag1974) suggested using this idea to determine whether there is a multivariate distribution that is consistent with certain univariate full conditional distributions. Suggala et al. (Reference Suggala, Yang, Ravikumar, Precup and Teh2017) use this idea to analyze three popular models for ordinal variables: the cumulative ratio, the continuation ratio, and the consecutive ratio (or adjacent category) logits (Agresti, Reference Agresti2010, Reference Agresti2018). Theorem 3 in Suggala et al. (Reference Suggala, Yang, Ravikumar, Precup and Teh2017) shows that there is a multivariate distribution consistent with the consecutive ratio or adjacent categories logit,

$$\begin{align*}\log\left( \frac{P(X = h)}{P(X = h + 1)} \right) = \eta_h - \beta\text{, for }h = 0\text{, }\dots\text{, }m - 1. \end{align*}$$

We use this logit to define the full conditional of a node i given the remaining variables $\mathbf {X}^{(i)}$ , by setting $\eta _h = \eta _{ih}$ and $\beta = \beta _i(\mathbf {X}^{(i)})$ , such that the full conditional is of the form

$$\begin{align*}P(X_i = h \mid \mathbf{X}^{(i)}) = \frac{\exp\left(\sum_{q = 0}^{h} (\eta_{iq} - \beta_i(\mathbf{X}^{(i)}))\,\mathbb{I}(x_i \leq q)\right) }{1+\sum_{u=0}^{m-1}\exp\left(\sum_{q=0}^{m-1}(\eta_{iq} - \beta_i(\mathbf{X}^{(i)}))\,\mathbb{I}(u \leq q)\right)} \text{, for }h = 0\text{, }\dots\text{, }m - 1. \end{align*}$$

Theorem 3 in Suggala et al. (Reference Suggala, Yang, Ravikumar, Precup and Teh2017) shows that the multivariate distribution that is consistent with this full conditional is equal to

$$ \begin{align*} P(\mathbf{X} \mid \boldsymbol{\eta}\text{, }\boldsymbol{\Sigma}) &= \frac{1}{\text{Z}} \exp\left( \sum_{i=1}^p\sum_{q=0}^{m-1} \eta_{iq}\mathbb{I}(X_i\leq q) + \sum_{i=1}^{p-1}\sum_{j=i+1}^p\sum_{q=0}^{m-1}\sum_{u=0}^{m-1}\sigma_{ij}\mathbb{I}(X_i\leq q)\mathbb{I}(X_j\leq u) \right) \\&= \frac{1}{\text{Z}} \exp\left( \sum_{i=1}^p\sum_{q=0}^{m-1} \eta_{iq}\mathbb{I}(X_i\leq q) + \sum_{i=1}^{p-1}\sum_{j=i+1}^p\sigma_{ij}(m-X_i)(m-X_j) \right). \end{align*} $$

When we redefine $\mu _{ih}$ as $\sum _{q = h}^{m - 1}\eta _{iq}$ and $X_i^\ast $ as $m-X_i$ , this multivariate distribution is equal to the MRF in Eq. (2). Thus, the proposed MRF in Eq. (2) is a multivariate extension of the adjacent category model.

Suggala et al. (Reference Suggala, Yang, Ravikumar, Precup and Teh2017) showed that no multivariate distribution is consistent with either of the other two models for ordinal variables. The cumulative ratio or proportional odds is defined as

$$\begin{align*}\log\left( \frac{p(X \leq h)}{p(X> h)} \right) = \eta_h - \beta\text{, for }h = 0\text{, }\dots\text{, }m - 1. \end{align*}$$

This logit model is similar in nature to the probit model (Guo et al., Reference Guo, Levina, Michailidis and Zhu2015), in which the realizations of the ordinal variable correspond to adjacent intervals of an underlying latent variable. Theorem 1 in Suggala et al. (Reference Suggala, Yang, Ravikumar, Precup and Teh2017) shows that no multivariate distribution is consistent with having this logit model as full-conditional distribution. Furthermore, their Theorem 2 shows that there is also no multivariate distribution with full conditionals based on the continuation ratio:

$$\begin{align*}\log\left( \frac{p(X = h)}{p(X> h)} \right) = \eta_h - \beta\text{, for }h = 0\text{, }\dots\text{, }m - 1. \end{align*}$$

This implies that the proposed MRF is a rather unique multivariate model for ordinal variables.

2.3 Markov properties

The model in Eq. (2) is a Markov random field model in which the manifest variables are viewed as nodes of a graph or network whose structure is described by the matrix of pairwise interaction effects or edge weights $\boldsymbol {\Sigma }$ . If $\sigma _{ij} = 0$ , then the variables i and j do not directly interact in the network. In Appendix A, we show that the model satisfies a global Markov property. A convenient feature of this Markov property is that the structure of the graph—characterized by the parameters $\boldsymbol {\Sigma }$ —represents the conditional dependence relations between variables in the network. This means that if $\sigma _{ij} = 0$ , variables i and j are conditionally independent given the rest of the variables in the network. In other words, the remaining variables $\mathbf {X}^{(i\text {, }j)}$ fully account for a possible association between $X_i$ and $X_j$ .

One way to examine the properties of the model is to focus on its local properties, such as those implied by its full conditional distribution, the distribution of one or more variables given all other variables in the network. Appendix A shows the general form of conditional distributions for the ordinal MRF. Here we focus on the conditional distribution of one of the variables in the network:

(4) $$ \begin{align} p(X_i = x_i \mid \mathbf{x}^{(i)}) = \frac{\exp\left(\sum_{h=1}^{m}\mathbb{I}(x_i = h)\,\mu_{ih} + x_i\,\sum_{j\neq i}\sigma_{ij}\, x_j\right)}{1 + \sum_{u=1}^{m}\exp\left(\mu_{iu} + u\,\sum_{j\neq i}\sigma_{ij}\, x_j\right)}, \end{align} $$

which is similar to a logistic regression for ordinal outcomes using the remaining variables as predictors. Here we see that when the interaction between variables i and j is positive, responding in a higher category for variable i also tends to lead to responding in a higher category for variable j. And if there is no direct relationship between variables i and j (i.e., $\sigma _{ij}=0$ ), we see that variable j drops from the conditional distribution of variable i, a consequence of the Markov property of the model. In the absence of interactions, i.e., $\sigma _{ij}=0$ for all $j \neq i$ , the full conditional distribution simplifies to a categorical distribution (Agresti, Reference Agresti2018) with parameters

$$\begin{align*}p(X_i = h\mid \mathbf{x}^{(i)}) = \theta_{h} = \frac{\exp\left(\mu_{ih}\right)}{1 + \sum_{u=1}^{m}\exp\left(\mu_{iu} \right)}. \end{align*}$$

The category threshold parameters for the ordinal MRF thus express the tendency of ordinal response variables that cannot be explained by the remaining variables in the graph or network.

2.4 Parameter identification and interpretation

The model is in the exponential family and has sufficient statistics for each parameter. Here we use this to establish the identification of the model parameters. Let $P_{ij}(h\text {, }q) = P(X_i = h\text {, }X_j=q\text {, }\mathbf {X}^{(i\text {, }j)})$ denote the probability that variables i and j in the network take values h and q, respectively. Then we can see that the interaction parameter $\sigma _{ij}$ is related to the ratio of the adjacent category odds of the two ordinates:

$$ \begin{align*} \frac{P_{ij}(h\text{, }q)}{P_{ij}(h+1\text{, }q)} \bigg/ \frac{P_{ij}(h\text{, }q+1)}{P_{ij}(h+1\text{, }q+1)} &= \frac{P_{ij}(h\text{, }q)}{P_{ij}(h\text{, }q+1)} \bigg/ \frac{P_{ij}(h+1\text{, }q)}{P_{ij}(h+1\text{, }q+1)}\\ &= \exp\left(2 \sigma_{ij}\right). \end{align*} $$

Thus, there is a direct mapping between the manifest variable—the ratio of observed category odds—and the interaction parameter of the ordinal MRF. The interaction parameter thus indicates the change in the logarithm of the adjacent category odds ratio. If $\sigma _{ij}> 0$ , the odds of observing one response in a higher adjacent category are greater than observing the other response in a higher adjacent category. If $\sigma _{ij} < 0$ , the odds of observing one response in a higher neighboring category are smaller than observing the other response in a higher neighboring category.

In the derivations above, we already assumed that one of the category threshold parameters for each variable is zero. We inherited the identifying restriction that the threshold for the lowest category is zero, i.e., $\mu _{i0} =0$ , from the GPCM. Let $P_i(h) = P(X_i = h\text {, }\mathbf {X}^{(i)})$ denote the probability that the variable i takes the value h. Then we can establish a relationship between the difference in category thresholds and the adjacent category odds:

$$\begin{align*}\frac{P_{i}(h+1)}{P_{i}(h)} = \exp\left(\mu_{i(h+1)} - \mu_{ih} + 2\,\sum_{j\neq i}\sigma_{ij}X_j\right). \end{align*}$$

Note that this relation’s dependence on the interaction parameters is fine since we established their identification without using the category thresholds. We could fill in their log-odds expressions here.

Observe that we can only identify the difference in adjacent category threshold parameters. However, equipped with the restriction that the threshold for the lowest category is zero, i.e., $\mu _{i0} =0$ , we identify $\mu _{i1}$ from

$$\begin{align*}\frac{P_{i}(1)}{P_{i}(0)} = \exp\left(\mu_{i1} + 2\,\sum_{j\neq i}\sigma_{ij}X_j\right), \end{align*}$$

which we can then use to identify the other category thresholds. If there are no interactions between variables in the network (i.e., $\boldsymbol {\Sigma } = \boldsymbol {0}$ ), $\mu _{i(h+1)} - \mu _{ih}> 0$ implies that there is a higher probability of observing responses in category $h+1$ than in category h. Conversely, if $\mu _{i(h+1)} - \mu _{ih} < 0$ , there is a higher probability of observing responses in category h than in category $h+1$ . In the face of network interactions, the term $2\,\sum _{j\neq i}\sigma _{ij}X_j$ modifies these response tendencies. Note that this term is constant across neighboring categories. In this case, $\mu _{i(h+1)} - \mu _{ih}> 2\,\sum _{j\neq i}\sigma _{ij}X_j$ means that there is a higher probability of observing responses in category $h+1$ than in category h, and $\mu _{i(h+1)} - \mu _{ih} < 2\,\sum _{j\neq i}\sigma _{ij}X_j$ implies the opposite.

4.1 The computationally intractable likelihood

A key problem in the analysis of multivariate categorical models is that their normalizing constant tends to be computationally intractable. For example, for $p=10$ variables and $m=5$ ordinal categories, the normalizing constant $\text {Z}$ in Eq. (3) would consist of $m^p = 5^{10} = 9,765,625$ terms, which we would have to evaluate repeatedly in each iteration of the proposed Gibbs sampler. This is clearly computationally infeasible.

In network psychometrics, the problem of an intractable likelihood is usually solved by adopting the pseudolikelihood approximation proposed by Besag (Reference Besag1975), which replaces the intractable likelihood with a tractable one. This approximation underlies most existing frequentist methods for regularized parameter estimation of networks of discrete variables (e.g., Haslbeck & Waldorp, Reference Haslbeck and Waldorp2020; van Borkulo et al., Reference van Borkulo, Borsboom, Epskamp, Blanken, Boschloo, Schoevers and Waldorp2014), but more recently pseudolikelihoods have also been incorporated into Bayesian methods (Marsman, Huth, et al. Reference Marsman, Huth, Waldorp and Ntzoufras2022; Mohammadi et al., Reference Mohammadi, Schoonhoven, Vogels and Birbil2023; Pensar et al., Reference Pensar, Nyman, Niiranen and Corander2017). However, in other fields, especially in the area of social network analysis, several MCMC methods have been proposed that avoid the need to evaluate the full likelihood and its normalizing constant. Park and Haran (Reference Park and Haran2018) provide a recent review of available approaches and show that the double Metropolis hastings (DMH) algorithm of Liang (Reference Liang2010) performs best overall in terms of effective number of samples from the posterior per second. In this article, we consider the pseudolikelihood approach but discuss the DMH algorithm and its implementation for Bayesian edge selection in Appendix C.

4.1.1 The pseudolikelihood approach

The pseudolikelihood approach approximates the joint distribution of the vector variable $\mathbf {x}$ —i.e., the full MRF in Eq. (2)—with its respective full-conditional distributions (see Eq. (4)):

$$ \begin{align*} P(\mathbf{X} \mid \boldsymbol{\mu}\text{, }\boldsymbol{\Sigma}) \propto P^\ast(\mathbf{X} \mid \boldsymbol{\mu}\text{, }\boldsymbol{\Sigma}) &= \prod_{i=1}^p \frac{\exp\left(\sum_{h=1}^{m}\mathbb{I}(X_i = h)\,\mu_{ih} + X_i\,\sum_{j\neq i}\sigma_{ij}\, X_j\right)}{1 + \sum_{u=1}^{m}\exp\left(\mu_{iu} + u\,\sum_{j\neq i}\sigma_{ij}\, X_j\right)}\\&= \frac{\exp\left(\sum_{i=1}^p\sum_{h=1}^{m}\mathbb{I}(X_i = h)\,\mu_{ih} + 2\sum_{i=1}^{p-1}\sum_{j = i+1}^pX_iX_j\sigma_{ij}\right)}{\prod_{i=1}^p\left(1 + \sum_{u=1}^{m}\exp\left(\mu_{iu} + u\,\sum_{j\neq i}\sigma_{ij}\, X_j\right)\right)}. \end{align*} $$

Note that the pseudolikelihood is equivalent to the full likelihood in Eq. (2), except that it replaces the intractable normalization constant with a tractable one. For the ten-variable network we considered earlier, the pseudolikelihood normalization constant has only $m \times p = 5 \times 10 = 50$ terms instead of the $9,765,625$ terms in case of the full likelihood.

Although the pseudolikelihood is a crude approximation to the full likelihood, the maximum pseudolikelihood estimates are consistent (Arnold & Strauss, Reference Arnold and Strauss1991; Geys et al., Reference Geys, Molenberghs and Ryan2007), and this is also the case for the pseudoposterior distribution (Miller, Reference Miller2021). For exponential random graph models, a random graph model used in social network analysis, van Duijn et al. (Reference van Duijn, Gile and Handcock2009) showed that pseudolikelihood estimates are generally biased in scenarios where there is only one observation of the network. For the Ising model, a graphical model used in network psychometrics, Keetelaar et al. (Reference Marsman, Sekulovski and van den Bergh2024) showed that the bias in the pseudolikelihood estimates for the Ising model is comparable to that of the MLE in several scenarios typical of psychological applications. In these applications, we typically have many observations of the network. The overall good performance of the pseudolikelihood established by Keetelaar et al. (Reference Keetelaar, Sekulovski, Borsboom and Marsman2024) applies especially to the joint pseudolikelihood approximation, the approximation we consider here, and less so to the disjoint pseudolikelihood approximation, also known as nodewise regression. Estimates based on the disjoint pseudolikelihood approximation can be severely biased, especially when the sample size is small. To verify that Keetelaar et al.’s results for the Ising model extend to our ordinal MRF, we examine the bias of the pseudolikelihood estimates of the parameters of the ordinal MRF relative to their maximum likelihood estimates in Appendix B and confirm that the bias of the two estimators is indeed comparable.

But there is no free lunch for the pseudolikelihood approach. Although it is fast, consistent, and also accurate in settings encountered in psychological applications, its standard errors can be underestimated. This was demonstrated by Keetelaar et al. (Reference Keetelaar, Sekulovski, Borsboom and Marsman2024) for the Ising model. Despite this limitation, the pseudolikelihood approach can consistently estimate the unknown graph structure. Several studies on graph recovery using pseudolikelihood have established its consistency in high-dimensional settings where both n and p are allowed to grow (e.g., Barber & Drton, Reference Barber and Drton2015; Meinshausen & Bühlmann, Reference Meinshausen and Bühlmann2006; Ravikumar et al., Reference Ravikumar, Wainwright and Lafferty2010). Csiszár and Talata (Reference Csiszár and Talata2006) proposed an extension of BIC for pseudolikelihoods and showed that it can consistently reveal the Markov structure of the generating MRF in cases where the graphs grow in size but nodes have a finite number of neighbors. Pensar et al. (Reference Pensar, Nyman, Niiranen and Corander2017) used PIC to show that, under very general conditions, the marginal pseudolikelihood can consistently reveal the Markov cover of each node (the set of all neighbors of the node) and the global graph structure as the sample size increases (see also Vogels et al., Reference Vogels, Mohammadi, Schoonhoven and Birbilin press, for a practical illustration).

4.2 A metropolis within Gibbs approach to Bayesian edge selection

The Gibbs sampler is the standard solution for sampling values from an intractable posterior distribution. As indicated above, the discontinuity in the prior distribution for the interactions creates a serious complication for the Gibbs sampler. To explain why this is the case, suppose we use the following block updating scheme:

$$ \begin{align*} \boldsymbol{\mu} &\sim f(\boldsymbol{\mu} \mid \mathbf{X}\text{, }\boldsymbol{\Sigma}),\\ \boldsymbol{\Sigma} &\sim f(\boldsymbol{\Sigma} \mid \mathbf{X}\text{, }\boldsymbol{\mu}\text{, }\boldsymbol{\gamma}),\\ \boldsymbol{\gamma} &\sim f(\boldsymbol{\gamma} \mid \boldsymbol{\Sigma}). \end{align*} $$

If we use this scheme and simulate the value $\gamma _{ij} = 0$ , we have to set $\sigma _{ij} = 0$ . But then we cannot move to other states, since $P(\gamma _{ij} = 1 \mid \sigma _{ij} = 0) = 0$ , and so $\gamma _{ij}$ and $\sigma _{ij}$ would have to remain at zero. The Markov chain is then said to be reducible (e.g., Norris, Reference Norris1998).

As the Markov chain moves between two structures, we either exclude edges from the network and set their nonzero weights to zero, or include edges in the network and set the corresponding zero-value weights to nonzero values. As a result of these moves between models, the parameter dimensions change, and the MCMC method must account for these changes in order to generate a proper, irreducible Markov chain. Transdimensional MCMC methods have been proposed in the literature to account for these changes (e.g., Carlin & Chib, Reference Carlin and Chib1995; Green, Reference Green1995; Sisson, Reference Sisson2005). The most popular transdimensional MCMC method is the reversible jump algorithm (Fan & Sisson, Reference Fan, Sisson, Brooks, Gelman, Jones and Meng2011; Green, Reference Green1995), a Metropolis algorithm that uses auxiliary variables to adjust parameter dimensions as it moves between models. However, the reversible jump algorithm is difficult to tune and sensitive to how we adjust the parameter dimensions. Gottardo and Raftery (Reference Gottardo and Raftery2008) viewed the discrete spike and slab as a mixture of mutually singular (MoMS) distributions—probability distributions that have disjoint support (see also Petris & Tardella, Reference Petris and Tardella2003, for earlier ideas)—and showed that reversible jump is similar to a Metropolis algorithm that jointly updates the indicator variable and the focal parameter (here $\sigma _{ij}$ ). The advantage of this approach is that there is no need to specify a method for dimension matching as in reversible jump, and we can update the pairs $(\gamma _{ij}\text {, }\sigma _{ij})$ using the Metropolis algorithm. We embed the Metropolis pair update approach into the Gibbs sampler.

The proposed Gibbs sampler comprises two blocks of parameters to update

$$ \begin{align*} \boldsymbol{\mu} &\sim f(\boldsymbol{\mu} \mid \mathbf{X}\text{, }\boldsymbol{\Sigma}),\\ \boldsymbol{\gamma}\text{, }\boldsymbol{\Sigma} &\sim f(\boldsymbol{\gamma}\text{, }\boldsymbol{\Sigma} \mid \mathbf{X}\text{, }\boldsymbol{\mu}). \end{align*} $$

Although updating all parameters in a block at once is preferable because it minimizes autocorrelation, we choose to update the category thresholds one at a time and investigate two different strategies for updating the $(\gamma _{ij}\text {, }\sigma _{ij})$ pairs, also one pair at a time. For the category thresholds, we do this because we can formulate an efficient sampler for the individual parameters. For the $(\gamma _{ij}\text {, }\sigma _{ij})$ pairs, we do this because the transition kernels for multi-pair updates are difficult to design and implement.

Next, we discuss the proposed MoMS Gibbs sampler based on the pseudolikelihood approximation, called PL-MoMS, which is included in the R package bgms (Marsman et al., Reference Marsman, Sekulovski and van den Bergh2023) available from CRAN: https://CRAN.R-project.org/package=bgms. Appendix D discusses the proposed MoMS Gibbs sampler based on the DMH algorithm of Liang (Reference Liang2010), called DMH-MoMS, which is implemented in the R package dmhBGM (Marsman et al., Reference Marsman, Sekulovski and van den Bergh2024), which is available on Github: https://github.com/MaartenMarsman/dmhBGM.

4.2.1 The PL-MoMS Gibbs sampler

Block I: Updating the category thresholds. To update the category thresholds, we need to sample from conditional distributions of the form

$$\begin{align*}f\left(\mu_{ik} \mid \mathbf{X}\text{, }\boldsymbol{\Sigma}\text{, }\boldsymbol{\mu}_i^{(k)}\right) \propto \frac{\exp\left(\mu_{ik} \right)^{n_{ik}}}{\prod_{v=1}^n\left(\text{g}_v + \text{q}_v\exp\left(\mu_{ik}\right)\right)}\, \frac{\exp\left(\mu_{ik}\right)^\alpha}{\left(1+\exp\left(\mu_{ik}\right)\right)^{\alpha +\beta}}, \end{align*}$$

for $k = 1\text {, }\dots \text {, }m$ , and $i = 1\text {, }\dots \text {, }p$ . Here, we have used $n_{ik} = \sum _{v=1}^n\mathbb {I}(x_{vi} = k)$ , $\text {g}_v = 1 + \sum _{u \neq k}\exp \left (\mu _{iu} + u\,\sum _{j\neq i}\sigma _{ij}\, x_{vj}\right )$ , and $\text {q}_v = \exp \left (k\,\sum _{j\neq i}\sigma _{ij}\, x_{vj}\right )$ . These full conditionals are intractable, but note that their form resembles that of a generalized beta-prime distribution

$$\begin{align*}f(\mu_{ik}) \propto \frac{\left(\text{c}\,\exp\left(\mu_{ik}\right)\right)^{\text{a}}}{\left(1+\text{c}\,\exp\left(\mu_{ik}\right)\right)^{\text{a}+\text{b}}}. \end{align*}$$

We will use the generalized beta-prime as a proposal in an independence chain Metropolis algorithm (Tierney, Reference Tierney1994). Maris et al. (Reference Maris, Bechger and San Martin2015) suggested using the properties of the target distribution to set the parameters of this proposal. This idea has also been successfully used by Marsman, Huth, et al. (Reference Marsman, Huth, Waldorp and Ntzoufras2022) and Marsman and Huth (Reference Marsman and Huth2023) to estimate the threshold parameters of the Ising and Divide and Color model (Häggström, Reference Häggström2001), respectively.

The log of the target distribution is concave and has linear tails:

$$\begin{align*}\frac{\text{d}}{\text{d} \mu_{ik}}\, \ln\left\lbrace f\left(\mu_{ik} \mid \mathbf{X}\text{, }\boldsymbol{\Sigma}\right) \right\rbrace \longrightarrow \begin{cases} n_{ik} - n -\beta &\text{as }\mu_{ik} \rightarrow \infty,\\ n_{ik} + \alpha &\text{as }\mu_{ik} \rightarrow -\infty, \end{cases} \end{align*}$$

and the same holds for the proposal distribution:

$$\begin{align*}\frac{\text{d}}{\text{d} \mu_{ik}}\, \ln\left\lbrace f\left(\mu_{ik} \right) \right\rbrace \longrightarrow \begin{cases} -\text{b} &\text{as }\mu_{ik} \rightarrow \infty,\\ \text{a} &\text{as }\mu_{ik} \rightarrow -\infty. \end{cases} \end{align*}$$

We match the tails of the proposal distribution to that of the target distribution by setting $\text {a} = n_{ik} + \alpha $ and $\text {b} = \beta + n - n_{ik}$ . The last free parameter, $\text {c}$ , is used to ensure that the proposal closely matches the target distribution at the current state of the Markov chain. Specifically, $\text {c}$ is used to make the derivatives of the logarithms of the proposal and target distributions equal. If $\hat {\mu }$ is the current state of $\mu _{ik}$ in the Markov chain, then we equate the two derivatives and solve c, which yields

$$\begin{align*}\text{c} = \frac{\sum_{v=1}^n\frac{\text{q}_v}{\text{g}_v+\text{q}_v\,\exp(\hat{\mu})} + (\alpha+\beta) \, \frac{1}{1+ \,\exp(\hat{\mu})}}{\alpha+\beta+n - \sum_{v=1}^n\frac{\text{q}_v\,\exp(\hat{\mu})}{\text{g}_v+\text{q}_v\,\exp(\hat{\mu})} - (\alpha+\beta) \, \frac{\exp(\hat{\mu})}{1+ \,\exp(\hat{\mu})}}. \end{align*}$$

Now that we have the value for $\text {c}$ , we can sample a proposal from the generalized beta-prime distribution in the following way: we sample Y from a Beta $(\text {a}\text {, }\text {b})$ distribution and set the proposed value $\mu ^\prime $ equal to $ \log \left (\text {c}^{-1}\frac {Y}{1-Y}\right )/2$ . Here, $\mu ^\prime $ is a sample from the generalized beta-prime distribution. We accept the proposed value with probability

$$\begin{align*}\theta = \min\left\lbrace 1\text{, } \frac{\prod_{v=1}^n\left(\text{g}_v + \text{q}_v\exp\left(\hat{\mu}\right)\right)}{\prod_{v=1}^n\left(\text{g}_v + \text{q}_v\exp\left(\mu^\prime\right)\right)}\, \frac{\left(1+\exp\left(\hat{\mu}\right)\right)^{\alpha +\beta}}{\left(1+\exp\left(\mu^\prime\right)\right)^{\alpha +\beta}} \frac{\left(1+\text{c}\,\exp\left(\mu^\prime\right)\right)^{\alpha+\beta+n}}{\left(1+\text{c}\,\exp\left(\hat{\mu}\right)\right)^{\alpha+\beta+n}} \right\rbrace, \end{align*}$$

and retain the current value $\hat {\mu }$ otherwise.

Block II: Updating the edge indicators and interactions. We base the proposal for Metropolis on the current state of the edge and interaction (i.e., $\gamma ^\ast $ and $\sigma ^\ast $ ), and we factor our proposal as

$$\begin{align*}f(\sigma^\prime\text{, }\gamma^\prime \mid \sigma^\ast\text{, }\gamma^\ast) = f(\sigma^\prime \mid \sigma^\ast\text{, }\gamma^\prime) \, p(\gamma^\prime \mid \gamma^\ast). \end{align*}$$

We will also use the MoMS formulation for the two proposal distributions. First, we consider the proposal distribution for the edge indicator,

$$\begin{align*}p(\gamma^\prime \mid \gamma^\ast) = (1 - \gamma^\ast)\, 1_{\{1\}}(\gamma^\prime) + \gamma^\ast\, 1_{\{0\}}(\gamma^\prime), \end{align*}$$

which proposes to change the edge, i.e., it sets $\gamma ^\prime = 1 - \gamma ^\ast $ to make the between model move. Next, we define the conditional proposal distribution for the interaction effect,

$$\begin{align*}f(\sigma^\prime \mid \sigma^\ast\text{, }\gamma^\prime) = (1 - \gamma^\prime) \, 1_{\{0\}}\, (\sigma^\prime) + \gamma^\prime \, 1_{\mathbb{R} \backslash \{0\}}\, (\sigma^\prime) \, f_{\text{proposal}}(\sigma^\prime \mid \sigma^\ast). \end{align*}$$

Here, the proposed value $\sigma ^\prime $ is equal to $0$ if $\gamma ^\prime = 0$ , and otherwise it is drawn from a proposal density $f_{\text {proposal}}(\sigma \mid \sigma ^\ast )$ otherwise. We use a normal proposal distribution centered on the current state (i.e., a random walk) with variance $\nu $ . We need to specify this variance parameter. We follow Kolovsky and Vanucci (Reference Kolovsky and Vanucci2020) and Wadsworth et al. (Reference Wadsworth, Argiento, Guindani, Galloway-Pena, Shelburne and Vanucci2017) and use adaptive Metropolis to learn the variance $\nu $ from the past performance of the Markov chain.Footnote ⁴ We are now ready to formulate our Metropolis–Hastings approach.

Let $\gamma ^\ast $ and $\sigma ^\ast $ denote the current state of the edge indicator and the interaction between a variable i and j, and let $\gamma ^\prime $ and $\sigma ^\prime $ denote its proposed state. We accept the proposed states with probability

$$\begin{align*}\pi = \min\left\lbrace 1\text{, } \overbrace{\frac{p^\ast(\mathbf{X} \mid \boldsymbol{\Sigma}^\prime\text{, }\boldsymbol{\mu})}{p^\ast(\mathbf{X} \mid \boldsymbol{\Sigma}^\ast\text{, }\boldsymbol{\mu})}}^{\substack{\text{Pseudolikelihood}\\\text{ratio}}}\, \times\, \overbrace{\frac{f(\sigma^\prime \mid \gamma^\prime)\,p(\gamma^\prime)}{f(\sigma^\ast \mid \gamma^\ast)\,p(\gamma^\ast)}}^{\substack{\text{Prior ratio}}}\, \times\, \overbrace{\frac{f(\sigma^\ast \mid \sigma^\prime\text{, }\gamma^\ast)\, p(\gamma^\ast \mid \gamma^\prime)}{f(\sigma^\prime \mid \sigma^\ast\text{, }\gamma^\prime)\, p(\gamma^\prime \mid \gamma^\ast)}}^{\text{Proposal ratio}}\right\rbrace, \end{align*}$$

where $\boldsymbol {\Sigma }^\ast $ is the current state of the pairwise interaction matrix, and $\boldsymbol {\Sigma }^\prime $ is the proposed state, with elements in row i, column j, and row j, column i, set equal to the proposed value $\sigma ^\prime $ . If $\gamma ^\ast =0$ , we propose $\gamma ^\prime = 1$ and $\sigma ^\prime \sim N(0\text {, }\nu _{ij})$ and accept the proposed values with probability

$$\begin{align*}\pi = \min\left\lbrace 1\text{, } \frac{p^\ast(\mathbf{X} \mid \boldsymbol{\Sigma}^\prime\text{, }\boldsymbol{\mu})}{p^\ast(\mathbf{X} \mid \boldsymbol{\Sigma}^\ast\text{, }\boldsymbol{\mu})}\, \times\, \frac{f_{\text{slab}}(\sigma^\prime)}{f_{\text{proposal}}(\sigma^\prime \mid \sigma^\ast)}\right\rbrace. \end{align*}$$

Conversely, if $\gamma ^\ast =1$ , we propose $\gamma ^\prime = 0$ and $\sigma ^\prime =0$ . We accept $\gamma ^\prime $ and $\sigma ^\prime $ with probability

$$\begin{align*}\pi= \min\left\lbrace 1\text{, } \frac{p^\ast(\mathbf{X} \mid \boldsymbol{\Sigma}^\prime\text{, }\boldsymbol{\mu})}{p^\ast(\mathbf{X} \mid \boldsymbol{\Sigma}^\ast\text{, }\boldsymbol{\mu})}\, \times\, \frac{f_{\text{proposal}}(\sigma^\ast \mid \sigma^\prime)}{f_{\text{slab}}(\sigma^\ast)}\right\rbrace. \end{align*}$$

The pseudolikelihood ratios in the expressions above boil down to

$$ \begin{align*} \frac{p^\ast(\mathbf{X} \mid \boldsymbol{\Sigma}^\prime\text{, }\boldsymbol{\mu})}{p^\ast(\mathbf{X} \mid \boldsymbol{\Sigma}^\ast\text{, }\boldsymbol{\mu})} &= \exp\left(2(\sigma^\prime-\sigma^\ast) \sum_{v=1}^n x_{vi}x_{vj}\right)\\ &\quad\times \frac{ \prod_{v=1}^n\left(1 + \sum_{u=1}^{m}\exp\left(\mu_{iu} + u\,\left\lbrace x_{vj}\sigma^\ast + r_{vi}\right\rbrace \right)\right) }{ \prod_{v=1}^n\left(1 + \sum_{u=1}^{m}\exp\left(\mu_{iu} + u\,\left\lbrace x_{vj}\sigma^\prime + r_{vi}\right\rbrace \right)\right) }\\ &\quad\times \frac{ \prod_{v=1}^n\left(1 + \sum_{u=1}^{m_j}\exp\left(\mu_{ju} + u\,\left\lbrace x_{vi}\sigma^\ast + r_{vj}\right\rbrace \right)\right) }{ \prod_{v=1}^n\left(1 + \sum_{u=1}^{m_j}\exp\left(\mu_{ju} + u\,\left\lbrace x_{vi}\sigma^\prime + r_{vj}\right\rbrace \right)\right) }, \end{align*} $$

where $r_{vi} = \sum _{h \neq j\neq i} \sigma _{ih}\, x_h$ and $r_{vj} = \sum _{h \neq j\neq i} \sigma _{jh}\, x_h$ .

Two implementations of MoMS Metropolis have been reported in the literature. The Gibbs scheme in Wadsworth et al. (Reference Wadsworth, Argiento, Guindani, Galloway-Pena, Shelburne and Vanucci2017) applies MoMS Metropolis to each indicator-parameter pair in turn, and the add-delete scheme in Kolovsky and Vanucci (Reference Kolovsky and Vanucci2020) applies it to a randomly selected indicator-parameter pair. The add-delete scheme adds an additional update of the “active” interactions (i.e., we update the $\sigma _{ij}^\ast $ for which $\gamma _{ij}^\ast = 1$ ) to speed up convergence (Gottardo & Raftery, Reference Gottardo and Raftery2008; Vanucci, Reference Vanucci, Tadesse and Vanucci2022). Our approach is a hybrid of the two MoMS Metropolis approaches; we first consider a Gibbs scheme in which we update each indicator-parameter pair in turn, and then do an additional update of the active interaction parameters. Here we use a random walk Metropolis, i.e. we propose $\sigma ^\prime _{ij}\sim N(\sigma _{ij}^\ast \text {, }\nu _{ij})$ . This scheme is implemented in the bgms package and is the one used in the numerical and empirical illustrations in the next two sections.

5.1 A comparison of evidence estimates

The pseudolikelihood-based estimates in Appendix B appear to have smaller variance than the stochastic estimates based on the exact likelihood, which could make the posterior distributions based on a pseudolikelihood less sensitive to the prior distribution than posterior distributions based on the exact likelihood. However, in our Bayesian edge selection procedure, the statistical evidence for the inclusion or exclusion of individual edges in the network is modeled by the spike and slab prior distributions on the pairwise interaction parameters. The reduced sensitivity to the prior distribution is likely to affect the estimated evidence for the inclusion or exclusion of edges in the network. Before analyzing the quality of our Bayesian edge selection procedure in the next two sections, we compare here how PL-MoMS-PL and DMH-MoMS estimate the statistical evidence.

In Bayesian edge selection, the evidence for the inclusion or exclusion of individual links in the network is measured by the inclusion probability. The posterior inclusion probability is the model-averaged probability of including an edge in the network, given the data

$$\begin{align*}P(\gamma_{ij} = 1 \mid \mathbf{X}) = \sum_{\boldsymbol{\gamma}: \gamma_{ij} = 1} \,P(\boldsymbol{\gamma} \mid \mathbf{X}). \end{align*}$$

It indicates the probability that the edge is included in the network that generated the observed data and is essential for Bayesian analysis of graphical models. The edge inclusion probability is used to define the median probability structure—a single, optimal structure for predicting new observations (Barbieri & Berger, Reference Barbieri and Berger2004)—and the inclusion Bayes factor. The inclusion Bayes factor is the change from the prior inclusion odds to the posterior inclusion odds,

$$\begin{align*}\text{BF}_{10}^{(ij)} = \frac{P(\gamma_{ij} = 1 \mid \mathbf{X})}{P(\gamma_{ij} = 0 \mid \mathbf{X})} \Bigg/ \frac{P(\gamma_{ij} = 1)}{P(\gamma_{ij} = 0)}. \end{align*}$$

It indicates how much more likely it is that the observed data came from a network that included the edge i–j than from one that did not. The exclusion Bayes factor $\text {BF}_{01}^{(ij)} = {1}/{\text {BF}_{10}^{(ij)}}$ gives the evidence for the conditional independence between nodes i and j in the data. Sekulovski et al. (Reference Sekulovski, Keetelaar, Huth, Wagenmakers, van Bork, van den Bergh and Marsmanin press) show that the inclusion Bayes factor, as a test for conditional dependence or independence of pairs of variables in the network, is robust to assumptions about the structure of the rest of the network. It also appears to be optimal for detecting evidence of conditional independence compared to Bayes factors that compare the same network structure with and without the link, or to credible interval-based tests.

We investigate how PL-MoMS and DMH-MoMS estimate the inclusion Bayes factors using simulated data, following the setup for the bias analysis in Appendix B, where we generated $500$ datasets, each with 300 observations ( $n = 300$ ) on 24 ordinal variables ( $p= 24$ ) with five response categories ( $m = 4$ ). The datasets were generated from the ordinal MRF in Eq. (2). For each variable, four threshold parameters were sampled from a uniform distribution ranging from $-2.0$ to $-0.5$ , sorted in decreasing order. The pairwise interaction parameters were sampled from a normal distribution with a mean of zero and a standard deviation of $1/25$ . For each of these $500$ datasets, we estimate the joint posterior distribution of the model parameters and the edge inclusion indicators with PL-MoMS and DMH-MoMS, as described in Appendix D. It is generally recommended to run the inner Gibbs sampler of the DMH algorithm for $n \times p = 7, 200$ iterations, for each model parameter and edge indicator and pairwise interaction pair (i.e., $m \times p + {p \choose {2}} = 372$ times) in each iteration of the Gibbs sampler. Clearly, this would make its application intractable, as confirmed by the runtime analysis in Appendix E. For this analysis, we therefore use a fixed number of 10 iterations for the inner Gibbs sampler. We ran PL-MoMS and DMH-MoMS for $20,000$ iterations for each of the $500$ datasets. Analysis of a single dataset using PL-MoMS took about 5 min, and using DMH-MoMS with ten iterations of the inner Gibbs sampler took about 11 hr on a single core of a MacBook Pro with an M1 chip. We used the Snellius supercomputer to distribute the analysis of the 500 datasets in batches across 124 cores. On average, it took about 16 min to analyze a single dataset on a single core on Snellius with PL-MoMS and about 19 hr with DMH-MoMS.

We computed the expected a posteriori (EAP) estimates or posterior means of the edge indicator variables (i.e., the inclusion probability) and the pairwise interaction parameters and averaged these estimates across the datasets. Figure 1 shows scatterplots of the averaged EAP estimates of the interaction parameters against the estimated inclusion probabilities, with estimates based on PL-MoMS in the left panel and those based on DMH-MoMS in the right panel. The gray dashed lines indicate the evidence thresholds corresponding to $\text {BF}_{10} = 10$ at the top, $\text {BF}_{10} = 1$ in the middle, and $\text {BF}_{10} = {1}/{10}$ at the bottom. Note that the two scatterplots show a similar relationship between the interaction effects and the inclusion probabilities, but that PL-MoMS tends to show more evidence for edge inclusion and less evidence for edge exclusion than DMH-MoMS, and we therefore expect PL-MoMS to show higher sensitivity and lower specificity than DMH-MoMS. Note that the estimates for the pairwise interactions shrink to zero more for DMH-MoMS than for PL-MoMS, which is a consequence of the posterior distribution based on the pseudolikelihood being less influenced by the prior distribution. Because the EAP estimates based on the pseudolikelihood are further from zero, the pseudolikelihood approach results in increased sensitivity or decreased specificity.

Figure 1 Scatterplots comparing EAP (expected a posteriori) estimates for the pairwise interaction parameter with the estimated inclusion probabilities, each averaged over 500 datasets. The left panel shows estimates based on the pseudolikelihood approach, while the right panel focuses on estimates based on the DMH algorithm. Each scatterplot includes dashed lines representing specific evidence thresholds under a uniform prior; the top line indicates an inclusion Bayes factor value of 10 (i.e., evidence for inclusion), the middle line indicates an inclusion Bayes factor value of 1 (i.e., absence of evidence), and the bottom line indicates an inclusion Bayes factor value of 0.1 (i.e., evidence for exclusion).

5.2 Bayesian edge selection with binary data

Now that we have a better understanding of how PL-MoMS and DMH-MoMS accumulate evidence, we want to evaluate their impact on Bayesian edge selection. We also want to compare our proposed PL-MoMS and DMH-MoMS procedures with two alternative approaches for the Ising model, the EBIC-Lasso approach of van Borkulo et al. (Reference van Borkulo, Borsboom, Epskamp, Blanken, Boschloo, Schoevers and Waldorp2014) and the Bayesian edge selection method of Park et al. (Reference Park, Jin and Schweinberger2022), which uses a continuous spike and slab prior and the DMH algorithm to sample from the posterior distribution. We refer to the latter as DMH-CSaS (i.e., DMH combined with a Continuous Spike and Slab). Based on the results of Sekulovski et al. (Reference Sekulovski, Keetelaar, Huth, Wagenmakers, van Bork, van den Bergh and Marsmanin press) that we discussed earlier, we expect that our discrete spike and slab prior approach will be similar to EBIC-Lasso in that the method will have a high specificity or true negative rate and a low sensitivity or true positive rate. However, given what we learned about PL-MoMS and DMH-MoMS in the previous section, we expect this effect to be partially mitigated by the pseudolikelihood approach, and thus expect PL-MoMS to result in lower specificity and higher sensitivity than DMH-MoMS. Finally, given the above expectations and the results in Park et al. (Reference Park, Jin and Schweinberger2022), we expect PL-MoMS and DMH-MoMS to have higher specificity but possibly lower sensitivity than DMH-CSaS.

We wanted to compare our methods with the DMH-CSaS method proposed by Park et al. (Reference Park, Jin and Schweinberger2022), but because we were unable to use their software implementation, we replicated their simulation setup for a direct comparison. This involved generating 500 datasets from the Ising model (i.e., the ordinal MRF with $m =1$ ), each with 24 variables and 300 observations. Each dataset was simulated with distinct parameter values: category threshold parameters were uniformly sampled between $-2$ and $-0.5$ , and for the 276 edges, 69 were assumed present. The interaction parameters for these edges were sampled from a uniform distribution, with 41 of them sampled uniformly between $.5$ and $2$ , and the remainder between $-1$ and $-0.5$ . Datasets containing variables with zero variance (67/500) were excluded from our analyses.

We estimated the joint posterior distribution of the model parameters and the edge inclusion indicators with PL-MoMS and DMH-MoMS. As in the previous simulations, the inner Gibbs sampler for DMH was run for ten iterations, and we ran our Gibbs samplers for $20,000$ iterations on each of the $500$ datasets. As before, the slab distribution is a Cauchy distribution with a scale of 2.5. We used the IsingFit R package (van Borkulo et al., Reference van Borkulo, Epskamp and Robitzsch2016) with package defaults and the AND rule to estimate the parameters of the Ising model using EBIC-Lasso. Our results with IsingFit are almost identical to those reported by Park et al. (Reference Park, Jin and Schweinberger2022), confirming that our results should be comparable.

We estimated the median probability model using the two MoMS Gibbs samplers, which includes an edge in the network if its corresponding inclusion probability exceeds .5, and excludes it from the network otherwise. For EBIC-Lasso, we selected the edges that had a nonzero estimate and excluded the others. For each method, we compare how well they were able to recover the generating network structure by computing its specificity or true negative rate —TN / (TN + FP)—, sensitivity or true positive rate —TP / (TP + FN)—, and the Rand index (Rand, Reference Rand1971) —(TN + TP) / (TN + FP + TP + FN). The results are in Table 1; we also copied the original results from Park et al. (Reference Park, Jin and Schweinberger2022).

Table 1 Performance metrics—specificity, sensitivity, and Rand index—of different edge selection methods applied to 500 simulated Ising model data sets, following the setup of Park et al. (Reference Park, Jin and Schweinberger2022)

Note: It includes the proposed PL-MoMS and DMH-MoMS approaches, the ${\mathrm {DMH}}_{10,n}$ -CSaS approach of Park et al. (Reference Park, Jin and Schweinberger2022) ( $^{\ast \ast }$ results copied from their Table 1), and the EBIC-Lasso method of van Borkulo et al. (Reference van Borkulo, Borsboom, Epskamp, Blanken, Boschloo, Schoevers and Waldorp2014).

Table 1 shows important differences in performance between the Bayesian edge selection methods and EBIC-Lasso. EBIC-Lasso performed best at identifying which edges were missing from the generating structures, i.e., has a high specificity, while DMH-CSaS performed worst at this task. DMH-MoMS was close to EBIC-Lasso on this metric, while PL-MoMS performed significantly worse. Conversely, DMH-CSaS performed best at identifying which edges were present in the generating structures, i.e., has a high sensitivity, while EBIC-Lasso performed worst on this metric. Interestingly, the DMH-MoMS method was not close to EBIC-Lasso on this metric, as it did much better, but it did not perform as well as PL-MoMS on this metric. These effects were as expected.

The methods thus trade off their specificity and sensitivity differently. Methods that perform well in terms of specificity tend to perform poorly in terms of sensitivity, and vice versa. We used the Rand index to measure how these trade-offs occur and to characterize the overall performance of the methods in recovering the generating network structure. The DMH-MoMS Gibbs sampler performed best on this overall metric with 85% correct identifications, while DMH-CSaS performed worse with 77% correct identifications. Interestingly, although EBIC-Lasso and PL-MoMS trade off specificity and sensitivity differently, they perform quite similarly on this metric with 80% correct identifications.

5.3 Bayesian edge selection with ordinal data

We have seen that PL-MoMS and DMH-MoMS are very capable of recovering the underlying network structure when applied to binary data. In the next analysis, we will see if this good performance generalizes to the analysis of ordinal data with $m>1$ . Ideally, we would like to keep the simulation setup for $m=1$ and $m>1$ as similar as possible to render results comparable. However, there are two complications with this. First, the parameters of the Ising model and the ordinal MRF are not standardized, and the generating parameter values we used for the Ising model lead to degenerate versions of the ordinal MRF. Like the Ising model, when the parameters of the ordinal MRF exceed certain critical thresholds, the model exhibits a phase transition and from that point on produces data with zero variance. Unfortunately, it is generally unknown what these critical thresholds are. We address this by rescaling the parameter values. This brings us to the second complication. Our Bayesian edge selection procedure imposes a fixed spike and slab prior distribution on the pairwise interaction parameters; in our previous analysis, we used a Cauchy with a scale of 2.5 as the slab distribution, but this prior does not account for this rescaling of the interaction parameters. Because of this, and because we expect smaller interaction effects as the number of categories increases, we expect the two MoMS procedures to show increased specificity and decreased sensitivity, all else being equal.

We investigate the Bayesian edge selection performance of PL-MoMS and DMH-MoMS on $500$ datasets, each with 300 observations ( $n = 300$ ) on 24 ordinal variables ( $p = 24$ ) with five response categories ( $m = 4$ ) simulated from the ordinal MRF in Eq. (2) as follows. For each variable, four threshold parameters were sampled from a uniform distribution ranging from $-2.0$ to $-0.5$ , sorted in decreasing order. As in the binary case, we assume that 69 out of 276 edges were present and the rest were absent. However, where the products $x_i\times x_j$ were zero or one in the binary case, they now range from zero to $m^2 = 16$ . To account for this change in value, we sampled the interaction parameters in the same way as before but multiplied them by a factor of $m^{-2} = {1}/{16}$ . Thus, the interaction parameters for the current edges were sampled from a uniform distribution, with 41 of them uniformly sampled between ${1}/{32}$ and ${1}/{8}$ , and the rest between $-{1}/{16}$ and $-{1}/{32}$ . There were no datasets containing variables with zero variance.

Similar to the analysis in the previous section, we estimated the joint posterior distribution of the model parameters and the edge inclusion indicators with PL-MoMS and DMH-MoMS. As before, the inner Gibbs sampler for DMH was run for ten iterations, and we ran our Gibbs samplers for $20,000$ iterations on each of the $500$ datasets. The slab distribution is again a Cauchy distribution with a scale of 2.5. We expect that using the same scale for the Cauchy slab distribution while shrinking the generating interaction parameters would lead to increased specificity and decreased sensitivity, a scaling effect that is essentially a form of the Jeffreys–Lindley paradox in Bayesian hypothesis testing (Jeffreys, Reference Jeffreys1961; Lindley, Reference Lindley1957). To investigate this, we estimate the joint posterior distribution using PL-MoMS with a rescaled Cauchy $(0\text {, }{5}/{32})$ slab distribution. We did not do this for the DMH-MoMS approach because i) we think we can judge the effect of rescaling on this procedure from the effect of rescaling on the PL-MoMS approach, and ii) running the DMH-MoMS approach twice is quite time consuming even on a supercomputer (it took us about 2.5 days to analyze the 500 datasets on Snellius) and expensive.

Table 2 shows that, as expected, when we do not correct the scale of the slab distribution for the increased pairwise interaction effects, or for the fact that the parameters must now shrink in size, the specificity of the variable selection methods increases while their sensitivity decreases compared to the results in the binary case shown in Table 1. It becomes more difficult to detect the smaller effects because the prior is much more diffuse compared to the size of the effect in the binary case. Of the absent edges, PL-MoMS correctly identified about 85% in the previous analysis and now 97%. Similarly, the DMH-MoMS approach, which had 95% correct before, now has a near perfect score on this metric. However, as mentioned above, they performed less well in detecting the edges that now have very small edge weights. Where PL-MoMs correctly identified 66% of the present edges, this has now shrunk to 42%, and for DMH-MoMs it has shrunk from 54% to 20%. While the relative performance of PL-MoMs compared to DMH-MoMs is the same as before, in that DMH-MoMs has higher specificity and lower sensitivity, the fact that DMH-MoMs has such low sensitivity now makes its overall performance of 80% less than PL-MoMs’ 83%.

Table 2 Performance metrics—specificity, sensitivity, and Rand index—of the PL-MoMS and the DMH-MoMS methods applied to 500 simulated ordinal MRF datasets

We also performed an analysis with PL-MoMS correcting for the rescaling of the size of the pairwise interactions. This correction worked perfectly in that, despite the additional parameters that need to be estimated for the ordinal model, it shows almost identical performance on all three metrics as before. The specificity was 85% in both scenarios, the sensitivity decreased slightly from 66% to 64%, and the overall performance of 80% was also the same in both scenarios. This result suggests that correcting for the number of categories in the scale of the slab distribution can effectively mitigate differences in the size of the interaction effects.

6.1 Structure estimation in the presence of uncertainty

The two classical methods each provide a single network estimate but do not quantify how certain we are about this structure. In contrast, our Bayesian variable selection approach in principle provides the uncertainty (i.e., posterior probability) associated with each possible structure. In this setting, the model or structure with the highest posterior probability is then often selected. Under certain prior specifications, choosing the model with the highest posterior probability is, in principle, the estimate provided by classical methods (see Marsman, Huth, et al. (Reference Marsman, Huth, Waldorp and Ntzoufras2022), for a brief discussion). However, selecting the model with the highest posterior probability can be problematic for several reasons. First, we must estimate the posterior structure probabilities, and since we have no analytic expression for their values and there are too many possible structures to enumerate, we are often uncertain about which structure has the highest posterior probability, especially if we are uncertain about the underlying structure. Second, even if we were certain about their exact values, it is often the case that the most likely structure is not that much more plausible than the second most likely structure. To address this latter issue, Barbieri and Berger (Reference Barbieri and Berger2004) suggest using the median probability model, as it has an optimal predictive value even in the face of model uncertainty. The median probability model selects only those variables (edges) with a posterior inclusion probability greater than $.5$ .

To measure our uncertainty about the underlying structure for the problem at hand, we tracked how many different structures the MCMC procedure visited. Assuming the unit information prior for the interaction effects, we ran the MCMC procedure for one million iterations, which visited $994,700$ different structures. The most likely structure accounted for less than $0.01$ percent (one hundredth of one percent) of the posterior probability. These results show that we are very uncertain about the structure of the network; many structures seem plausible for the data at hand, and no structure stands out. We therefore consider the median probability model shown in Figure 2(c). Note that the median probability model is more densely connected than the thresholded partial correlation network, but unlike the EBICglasso network, it retains the negative associations.

6.2 Edge and structural evidence

Importantly, the fact that we are massively uncertain about which overall network structure is correct does not mean that we cannot be highly certain about specific edges. For example, while a number of edges may be present in some structures and not in others, leading to high uncertainty about the overall structure, other edges may be present (or absent) in almost all structures, and we can therefore be highly certain about their presence (or absence). A major advantage of Bayesian model averaging is that it allows us to express this uncertainty locally in terms of the Bayes factors for edge inclusion. Can we interpret the missing edge as evidence of conditional independence, or should we be more cautious with this interpretation?

The Bayesian methodology allows us to answer this question in a straightforward way. Figure 3 shows a scatterplot of the model-averaged posterior means of the interaction parameters (x-axis) against the logarithm of the corresponding inclusion Bayes factor (y-axis). This illustrates how the value of the Bayes factor increases as the corresponding posterior distribution moves away from zero. In the figure, we interpret an inclusion Bayes factor less than ${1}/{10}$ as evidence of edge exclusion, i.e., conditional independence. These inclusion Bayes factors are colored red. Similarly, we interpret inclusion Bayes factors greater than $10$ as evidence of edge inclusion, i.e., conditional dependence. These inclusion Bayes factors are colored in blue. The blue triangles or arrows at the top indicate the inclusion Bayes factors that we estimate to be greater than $148$ . This means that we can be confident that these edges are present in the data generating model. On the other hand, for the many grey edges indicating absence of evidence, we likely would not want to draw strong conclusions about their presence or absence in the data-generating model. Another way to present these results is with the three network plots in Figure 4, which show for which relations there is evidence of absence (left panel), absence of evidence (middle panel), and evidence of presence (right panel).

Figure 3 The mean of the model-averaged posterior distribution of the interaction parameters (x-axis) plotted against the log of the corresponding inclusion Bayes factor (y-axis). Positive values of the log Bayes factors indicate evidence for inclusion, and negative values indicate evidence for exclusion. We use Bayes factor values between ${1}/{10}$ and $10$ to indicate weak (or no) evidence. For this analysis, we assumed a unit information prior on the interaction parameters. Triangles indicate log Bayes factors greater than five.

Figure 4 Three networks showing the level of evidence for each edge. In the network on the left, the edges reflect inclusion Bayes factors less than ${1}/{10}$ ; in the middle network, the edges reflect Bayes factors between ${1}/{10}$ and $10$ ; the edges in the network on the right reflect Bayes factors greater than $10.$ For this analysis, we specified a unit information prior on the interaction parameters.

The inclusion Bayes factors convey the evidence for individual edges. But they can also help us identify parts of the network about which we are confident and parts about which we are uncertain. For example, the right panel of Figure 4 shows that we have conclusive evidence for each of the edges between the symptoms “flashbacks,” “dreams,” and “intrusions.” Let $\mathcal {H}_1$ denote the hypothesis that the three variables form a clique (i.e., there is an edge between each variable), and let $\mathcal {H}_0$ denote its complement (i.e., at least one of the edges is missing). The Bayes factor in which we pit the two hypotheses against each other is equal to

$$\begin{align*}\text{BF}_{10} =\underbrace{\frac{ P(\mathcal{H}_1 \mid \text{data}) }{ P(\mathcal{H}_0 \mid \text{data}) }}_{\text{Posterior Odds}} \Big/ \underbrace{\frac{ P(\mathcal{H}_1}{ P(\mathcal{H}_0) }}_{\text{Prior Odds}}. \end{align*}$$

There are eight different configurations of edges between the three variables, one corresponding to $\mathcal {H}_1$ and seven corresponding to $\mathcal {H}_0$ . Since all structures are equally likely a priori, the prior probability is ${1}/{7}$ in favor of $\mathcal {H}_1$ . Our estimate for the posterior probability that the data come from a network in which the three variables form a clique is close to one, i.e., $P(\mathcal {H}_1 \mid \text {data}) = .999987$ , and the estimated posterior odds in favor of $\mathcal {H}_1$ are $76,922.08$ . The Bayes factor is thus estimated to be $\text {BF}_{10} = 76,922.08 \times 7 = 538,454.6$ , indicating massive evidence in favor of $\mathcal {H}_1$ . In contrast, we have inconclusive evidence between the symptoms “avoidth” (avoidance of thinking about a stressful experience), “amnesia,” and “flashbacks.” We saw all $2^3 = 8$ possible configurations of edges between these three variables, and no structure stood out.

6.3 Prior robustness

The qualitative conclusions drawn from the Bayesian analysis of the ordinal MRF may be sensitive to the choice of the prior distribution. To assess the robustness of the results, we perform a robustness analysis comparing the results of the unit information prior to the Cauchy prior with different scale values. We report the details of this analysis in the Appendix F. Some of the conclusions of our analysis would change depending on the choice of a more or less diffuse prior, especially the BFs that are close to our cutoffs of $10$ or ${1}/{10}$ . At the same time, however, we found that many results were robust to changes in our prior specifications, especially the BFs that are far from our chosen cutoffs. Of the 52 BFs that showed compelling evidence for conditional independence, i.e. edge exclusion, using the unit information prior, 74 showed compelling evidence after adopting the more diffuse Cauchy $(0\text {, }2.5)$ prior density. Thus, while for 22 edges we now conclude that we have evidence of absence rather than absence of evidence, for 52 edges we did not change our conclusions for either prior specification. Conversely, of the 35 BFs that showed compelling evidence for conditional dependence, i.e., edge inclusion, using the unit information prior, 32 showed compelling evidence after adopting the more diffuse Cauchy $(0\text {, }2.5)$ prior density. Thus, while we now conclude that we have absence of evidence rather than evidence of presence for three edges, we have not changed our conclusions for 32 edges under either prior specification. In sum, the unit information prior appears to be the more conservative choice in this analysis.