Bayesian multiple Gaussian graphical models for multilevel variables from unknown classes

Abstract

Gaussian graphical models are a powerful tool for investigating the conditional dependency structure between random variables by estimating sparse precision matrices and can infer networks among variables from multiple classes. Many studies assume that classes of observations are given and use methods to learn the network structures within one level (e.g. pathways or genes). In most cases, however, heterogeneous data may be obtained at different levels. Therefore, in this paper, we consider the learning of multiple connected graphs with multilevel variables from unknown classes. We estimate the classes of the observations from the mixture distributions by evaluating the Bayes factor and learn about the network structures by fitting a neighborhood-selection algorithm. This approach can be used to identify the class memberships and reveal the network structures for lower level and higher level variables simultaneously. Unlike most existing methods, which solve this problem using frequentest approaches, we assess an alternative and novel hierarchical Bayesian approach for incorporating prior knowledge. We demonstrate the unique advantages of our methods through several simulations. A breast cancer application shows that our model’s results can provide insight into biological studies.

Keywords

Bayesian method Gaussian graphical model gene co-expression network multilevel network precision matrix estimation

1. Introduction

Gaussian graphical models have been applied widely in science to describe the conditional dependencies between variables. Learning the graph topology via Gaussian graphical models (GGMs) is equivalent to estimating nonzero entries in the inverse covariance matrix (or precision matrix) $Σ^{- 1}$ . However, estimating a precision matrix is statistically challenging, especially for high-dimensional data. One can compute the maximum likelihood estimator (MLE) for a non-decomposable Gaussian graphical model, also known as covariance selection.⁷ In a high-dimensional case in which the number of observations is less than the number of variables, the sample covariance matrix is singular and thus non-invertible. Another problem will also arise because the MLE can hardly yield exact zeros in the precision matrix. In other words, the MLE does not result in sparse graphs, thus making it difficult to interpret for network analysis.

Numerous methods of estimating the covariance matrix in the high-dimensional case have been developed to learn sparse graphs through penalty functions that force zeros in the corresponding precision matrix. Meinshausen and Bühlmann¹⁶ proposed doing so by performing lasso regression of all nodes to a target node. A penalized log-likelihood approach can be optimized through a variety of methods.^1,19,23 The most popular one is known as graphical lasso (glasso), developed by Friedman et al.⁹ An alternative is to form a Bayesian graphical lasso, which uses Laplace priors on off-diagonal entries of the precision matrix.²² In recent years, researchers have realized that it is more efficient to learn multiple graphs together because the graphs may share certain characteristics. Guo et al.¹¹ learned multiple graphs jointly by decomposing precision matrices into a multiplication of common factors across groups and unique factors for each group. Danaher et al.⁶ generalized this jointly learning multiple graph approach and used two penalty functions: the fused lasso and the group lasso.

However, the standard formulation for learning GGMs assumes that classes of observations are given and that the network structures will be learned within one level (e.g. pathways or genes). However, this assumption might be too optimistic. In some cases, the data at different classes may be heterogeneous. Often, these different classes can have ambiguous categorizations. For example, present-day gene-expression data can involve hundreds of pathway and gene variables while lacking the case–control status. Hence, we propose a technique for inferring multiple connected networks from observations belonging to unknown classes following mixture distributions. For genomic data, it is reasonable that genes can be clustered or grouped into pathways for particular functions. On the other hand, pathways are not isolated either but instead work together to accomplish certain tasks. One may also find it necessary to study their connectivity or relationship patterns. To tackle the problem of multilevel networks, we propose using multiple GGMs for multilevel variables from unknown classes to not only estimate the class membership but also learn the networks at a different level.

In this paper, we develop multiple GGMs for multilevel variables from unknown classes under the Bayesian hierarchical framework in order to infer multiple connected networks from observations belonging to unknown classes following mixture distributions. Our goal is to explore conditional dependency structures among the variables of two levels simultaneously (i.e. the set and element levels) when classes are drawn from the mixture distributions. We investigate a two-step solution to this problem using Bayesian methods. The basic idea of the proposed model is to assess class membership first and then learn the multilevel networks.

Our approach has several novel features: (1) it can reveal networks for multilevel variables simultaneously; (2) it evaluates class membership of the observations from GGMs; and (3) it is fully model based and thus has nice probabilistic interpretations. Clustering data for GGMs is challenging, and only a few works so far have focused on the relevant problems.

The rest of the paper is organized as follows. In Section 2, we first describe Bayesian multiple Gaussian graphical models (denoted as BMGGMs) for multilevel variables from unknown classes. Section 3 proposes a probability model for inferring the high-level variable network from a lower level variable network. Sections 4 and 5 describe the Bayesian generative model and posterior inference, respectively. In Section 6, we illustrate the performance of the proposed model through simulated data. Section 7 demonstrates the model’s application to a case study. Section 8 contains concluding remarks.

2. Multiple Gaussian graphical model for unknown classes

2.1. Problem setup and notation

We start by describing the problem setup and the notation that will be used throughout this paper.

Suppose we are given a data set that has C possible heterogeneous classes with a total of P variables, which belong to K pre-specified groups. Each group has p_k variables such that $\sum_{k = 1}^{K} p_{k} = P$ . We refer to the groups as “sets” (or higher level variables) and the variables within the higher-level variables as elements (or lower level variables). Because our data have a hierarchical structure, we shall refer to these variables as “multilevel variables.”

We index classes by c, where $c = 1, \dots, C$ ; index the matrix entries for lower level variables by ${j, j'}$ , where $1 < j \neq j' < P$ ; and index the matrix entries for higher level variables by ${k, k'}$ , where $1 < k \neq k' < K$ .

2.2. Background and general formulation

For a given class c, we have an $n_{c} \times P$ data matrix, where n_c is the sample size of class c and P is the number of variables. In other words, the cth class contains n_c observations ${x_{c, 1}, \dots, x_{c, n_{c}}}$ , where $x_{c, i} = {x_{c, i, 1}, x_{c, i, 2}, \dots, x_{c, i, p}} \in ℝ^{p}, i = 1, \dots n_{c}$ . We assume that within each class c, $x_{c, 1}, \dots, x_{c, n_{c}} \in ℝ^{p}$ are independent and identically drawn from a P-variate multivariate Gaussian distribution within each class

x_{c, i} \sim N (μ_{c}, Σ_{c}), i = 1, \dots, n_{c}

where

μ_{c} = {μ_{c, 1}, \dots, μ_{c, p}} \in ℝ^{p}

is the mean vector and

Σ_{c}

is a symmetric and positive definite P × P matrix.

Furthermore, we assume that all of the data points are generated from a mixture of C-component Gaussian distributions with unknown class membership for each individual. Following standard practice, we introduce a C-dimensional indicator z_i to represent the latent state. Then, the prior distribution of z can be specified in terms of the mixing coefficients π_c by

p (z_{c, i} = c) = π_{c}

where

\sum_{c = 1}^{C} π_{c} = 1

As each base distribution in the mixture is a multivariate Gaussian distribution with mean $μ_{c}$ and covariance matrix $Σ_{c}$ , the joint distribution has the form

p (x_{c, i} | μ_{c}, Σ_{c}) = \sum_{c = 1}^{C} π_{c} N (x_{c, i} | μ_{c}, Σ_{c}) i = 1, \dots, n_{c}

2.3. Learning a sparse Gaussian graphical model

Our goal is to infer graphs with a sparse structure so that we can describe the relationships among elements. Equivalently, we can estimate the non-zeros in the inverse matrix $Σ_{c}^{- 1}$ , which correspond to the edges in the graphs

Σ_{c}^{- 1} = (\begin{matrix} ω_{c, 1, 1} & \dots & ω_{c, 1, P} \\ ⋮ & ⋱ & ⋮ \\ ω_{c, P, 1} & \dots & ω_{c, P, P} \end{matrix})

(1)

The basic idea of learning a sparse GGM is to represent the inverse covariance matrix with a modified Cholesky decomposition and then apply the neighborhood selection with priors that can lead to sparsity.

Specifically, in order to learn a sparse covariance matrix with statistically interpretable parameterization, we adopt a modified Cholesky decomposition²⁰ of the precision matrix $Σ_{c}^{- 1} = L_{c}^{⊤} D_{c} L_{c}$ , where $L_{c}$ is an upper triangular matrix and $D_{c}$ is a diagonal matrix. With this decomposition, the multivariate normal model can be transformed into a regression problem.

For a given element j, $x_{c, i, j}$ can be written as follows

x_{c, i, j} = μ_{c, j} + \sum_{j' < j} β_{c, j' j} (x_{c, i, j'} - μ_{c, i, j'}) + ϵ_{c, j}, j = 2, \dots, P,

where

ϵ_{c, j}

is the random error for a given class c and element j.

This expression can be expanded further to

\begin{array}{l} x_{c, i, 2} = μ_{c, 2} + β_{c, 12} (x_{c, i, 1} - μ_{c, i, 1}) + ϵ_{c, 2}, \\ x_{c, i, 3} = μ_{c, 3} + \sum_{j' = 1}^{2} β_{c, j' 3} (x_{c, i, j'} - μ_{c, i, j'}) + ϵ_{c, 3}, \\ ⋮ \\ x_{c, i, p} = μ_{c, p} + \sum_{j' = 1}^{p - 1} β_{c, j' p} (x_{c, i, j'} - μ_{c, i, j'}) + ϵ_{c, p} \end{array}

Through these formulations, the upper triangular L has the following entries of regression coefficients

L_{c}^{⊤} = (\begin{matrix} 1 \\ - β_{c, 12} & 1 \\ - β_{c, 13} & - β_{c, 23} & 1 \\ ⋮ & ⋮ & ⋱ \\ - β_{c, 1 P} & - β_{c, 2 P} & \dots & - β_{c, P P - 1} & 1 \end{matrix})

(2)

For simplicity, each dataset is centered such that $μ_{c} = 0$ . From equation (2.3), we have $L_{c} X_{c} = ϵ_{c}$ . Define $D_{c} = var (ϵ_{c})$ ; then, $D_{c}$ is a diagonal matrix of variance. The advantage of this parameterization is that it enforces the symmetry and positive-definiteness of the precision matrix.

The precision matrix obtained from Meinshausen and Bühlmann¹⁶ is not necessarily positive definite. Other approaches^1,23 require that a precision matrix be positive definite. Therefore, to avoid these obstacles, we want the estimated precision matrix to be unconstrained and statistically meaningful. The variables of a modified Cholesky decomposition obviously depend on the order in which the variables appear in the data matrix. However, this problem becomes less important when we aim to find meaningful patterns across the classes, while $Σ_{c}^{- 1}$ in each class is guaranteed to be valid.

Note that we explore the ordering problem with a toy example. We generated data based on the following element network and corresponding precision matrix using the approach defined in Section 6.1.1. The diagram is as follows:

element 1 \to element 2 and element 2 \to element 3.

(See Figure 1 of the Supplementary Materials.)

Σ^{- 1} = (\begin{matrix} 1.00 & 0.50 & 0.00 \\ 0.50 & 1.00 & 0.50 \\ 0.00 & 0.50 & 1.00 \end{matrix})

Our BMGGMs give an estimation

{\hat{Σ}}^{- 1} = (\begin{matrix} 1.00 & 0.50 & - 0.01 \\ 0.50 & 1.25 & 0.54 \\ - 0.01 & 0.54 & 1.30 \end{matrix})

Now, we switch column 2 to column 3. Then, the true diagram becomes

element 1 \to element 3 and element 2 \to element 3.

Our BMGGMs give an estimation.

{\hat{Σ}}^{- 1} = (\begin{matrix} 1.00 & 0.00 & 0.53 \\ 0.00 & 1.00 & 0.52 \\ 0.53 & 0.52 & 1.56 \end{matrix})

The modified Cholesky decomposition gives a different estimation of the precision matrix but can detect the correct adjacency matrix, as shown in Figure 2 of the Supplementary Materials.

3. Multilevel network probability model

In this section, we propose a multilevel network probability model. We illustrate how to infer the network structure at the set level. In our example, the sets are high-level variables, and the elements are the lower level variables within the sets. We estimate set-interaction probabilities using element-interaction probabilities, by implementing the idea of Kim et al.,¹³ who originally proposed set and set interaction estimation through element and element probabilities.

We build a multilevel network model using the following assumption: (A1) two sets will interact if at least one element pair from the two sets interacts.

First, we define some notations. Define $g_{c, j, j'}^{(k, k')} = 1$ if element j interacts with element $j'$ in set pair k and $k'$ and $g_{c, j, j'}^{(k, k')} = 0$ otherwise. Then, $Pr (g_{c, j j'} = 1)$ is the probability that element j interacts with element $j'$ . Let $P_{c, k k'}$ represent the interaction event between set pair k and $k'$ , with $P_{c, k k'} = 1$ if they interact and $P_{c, k k'} = 0$ otherwise.

By using (A1), we can calculate the interaction probabilities between two sets using element interaction probabilities

Pr (P_{c, k k'} = 1) = 1 - \prod_{{g_{c, j, j'}^{(k, k')} \in P_{c, k k'}}} {1 - P r (g_{c, j, j'})^{ρ}},

(3)

where

{g_{c, j, j'}^{(k, k')} \in P_{c, k k'}}

is all pairs of elements from the sets k and

k'

and ρ is an adjusting constant value using the number of element pairs, which will be explained shortly.

We consider the following four conditions:

If $Pr (g_{c, j, j'}^{(k, k')} = 1) = 1$ for at least one lower level variable pair, then $Pr (P_{c, k k'} = 1) = 1$ .

If $Pr (g_{c, j, j'}^{(k, k')} = 1) = 0$ for all lower level pairs, then $Pr (P_{c, k k'} = 1) = 0$ .

If $Pr (g_{j, j'}^{(k, k')} = 1) = 1 / 2$ for all higher level pairs, then $Pr (P_{k k'} = 1) = 1 / 2$ .

If $0 < Pr (g_{c, j, j'}^{(k, k')} = 1) < 1$ , then $0 \leq Pr (P_{c, k k'} = 1) \leq 1$ .

Because the number of elements in a set varies, we adjust the set probability. The third condition is important for controlling the case in which the set probability (3) becomes larger as the number of element pairs becomes large. By using both equation (3) and this third condition, we can obtain the following equation

1 / 2 = 1 - \prod_{{g_{c, j, j'}^{k, k'} \in P_{c, k, k'}}} 1 - {(1 / 2)}^{ρ}

Then, ρ can be obtained by

ρ = \frac{\log {1 - {(1 / 2)}^{\frac{1}{M_{c, k, k'}}}}}{\log (1 / 2)}

where

M_{c, k, k'}

represents the total number of element pairs between sets k and

k'

This adjustment is motivated by the observation that the value of equation (3) increases as the number of elements increases. We also compared our result with ρ = 1. The result did not have significant changes. This similar result also was reported by Kim et al.¹³

4. Bayesian hierarchical framework

We seek to learn the network (both the set and element) structure for each class in order to obtain the inverse covariance matrix $Σ_{1}^{- 1}, \dots, Σ_{C}^{- 1}$ . Because we do not know which classes belong to which clusters, we estimate the class membership for each observation by fitting a finite mixture model.

In this section, we describe how to estimate the parameters under the Bayesian hierarchical framework and the Markov chain Monte Carlo (MCMC) algorithm to fit this model.

4.1. Estimating class membership using Gaussian mixture models

We treat class-membership estimation as a clustering problem. With the distribution of $p (x_{c, i} | μ_{c}, Σ_{c})$ , we can fit the mixture model to compute $p (z_{c, i} | x_{c, i})$ , which indicates the posterior probability that an entry belongs to the cth component given that $x_{c, i}$ observed in it. This value can be evaluated easily by applying Bayes’ theorem

p (z_{c, i} | x_{c, i}) = \frac{π_{c} N (x_{c, i} | μ_{c}, Σ_{c})}{\sum_{c' = 1}^{C} N (x_{c', i} | μ_{c'}, Σ_{c'})}

The cth mixing proportion π_c can be viewed as the prior probability that an entity belongs to the cth component of the mixture $c = 1, \dots, C$ . The posterior probability $p (x_{c, i} | μ_{c}, Σ_{c})$ is the corresponding conditional probability, known as the responsibility, once we have observed it.

Note that a label-switching problem arises if we want samples from the MCMC algorithm. That is, if we calculate $p (z_{c, i} | x_{c, i})$ for each iteration, we are free to relabel all of the classes without it affecting the data segments. This is normal because this is basically a clustering problem in which the parameter estimation does not rely on how the components are labeled. Ignoring this problem will result in severe consequences because the algorithm will never converge. One way to avoid this is to impose an identifiability constraint on the parameter space, such as

π_{1} \leq π_{2} \leq \dots \leq π_{C},

which is applied to our study, as described in Section 6.

4.2. Prior specification

In this subsection, we provide the specification of the priors in our models.

Let $θ_{c, k, k'}$ represent the prior probability of interaction between two sets, k and $k'$ . Let $g_{c, j, j'}$ be a latent indicator variable representing whether elements j and $j'$ are connected. These elements can be within the same set or across different higher levels. Specifically, if the pair of elements j and $j'$ comes from the same set k, then $θ_{c, k, k'}$ becomes $g_{c, j, j'}^{k}$ ; if the pair of elements j and $j'$ comes from different sets, k and $k'$ , then $θ_{c, k, k'}$ becomes $g_{c, j, j'}^{k, k'}$ . The matrix with entries $g_{j, j'}$ is called an adjacency matrix.

For each pair of elements $1 \leq j \leq j' \leq p$ , the prior distribution of the regression coefficient ( $β_{c, j, j'}$ ) for each edge depends on $g_{c, j, j'}$ and is considered to have the following distribution

p (β_{c, j, j'} | g_{c, j, j'}) = {\begin{array}{l} N (0, v_{1}^{2}), & if g_{j, j'} = 1 \\ N (0, v_{0}^{2}), & if g_{j, j'} = 0 \end{array}

(4)

where v₁ is far from 0 but v₀ is close to 0; then,

v_{1} \geq v_{0} > 0

If $g_{c, i, j} = 0$ , then $β_{c, j, j'}$ will have a prior with small variance v₀. We set $v_{1} = λ v_{0}$ . Parameter λ_c influences the prior’s probability of edge inclusion. In the context of Bayesian variable selection,¹⁰ λ_c should be chosen carefully so that it is large enough to support the edges in graphs but not large enough to support non-edges. We adopted the following suggestion in this paper: we set λ on the grid $10, 100, 1000, 1000, 1, 000, 000$ and chose λ = 1000 for illustration throughout this work. The prior called Bernoulli-Gaussian model has been used widely in works on Bayesian variable selection.¹⁴ This is very similar to the standard “spike and slab” prior but does not require summing over the “irrelevant” parameters. Thus, by using the Bernoulli-Gaussian model on the entries of the adjacency matrix, we can estimate the network structure and thus obtain the sparse inverse covariance matrix.

To facilitate the design of the MCMC scheme, we would like to require conjugate priors for the rest of the parameters. Thus, the prior distributions of $v_{0}^{2}, v_{1}^{2}, g_{c, j, j'}$ , and $θ_{c, k, k'}$ are specified as inverse-gamma, Bernoulli, and beta distributions, respectively, defined as

\begin{array}{l} p (v_{0}^{2}) = \frac{f^{e}}{Γ (e)} (v_{0}^{2})^{- e - 1} \exp (- \frac{f}{v_{0}^{2}}), \\ p (g_{c, j, j'} | θ_{c, k, k'}) = (θ_{c, k, k'})^{g_{c, j, j'}} (1 - θ_{c, k, k'})^{(1 - g_{c, j, j'})} \\ p (θ_{c, k, k'}) = \frac{Γ (a + b)}{Γ (a) Γ (b)} (θ_{c, k, k'})^{a - 1} {(1 - θ_{c, k, k'})}^{b - 1} \end{array}

With these priors, one can easily derive a simple Gibbs sampling scheme, which will be discussed in full in Section 6.

The choice of hyperparameters plays an important role in the performance of our BMGGMs. First, the variance $v_{0}^{2}$ depends on hyperparameters e and f. It increases as e decreases and f increases. Conversely, it decreases as e increases and f decreases. Then, a large value of f and small value of e will yield a sparse low-variable network. Preliminary cross-validation experiments showed that setting e = 10, 000 and f = 100 led to good performance with regard to adjacency matrix estimation for a given class. These settings will be applied in the simulations for all of the case studies. Therefore, we set the beta hyperparameters on $θ_{c, k, k'}$ to $a = 1, b = 2$ , regardless of whether $k = k'$ or $k \neq k'$ . A similar setting was used by Marlin and Murphy,¹⁵ but in our work, we want both the lower and higher variable networks to be sparse.

5. Posterior sampling

In this section, we describe a sampling scheme for BMGGMs. To estimate BMGGMs from unknown classes, we conduct two steps: we first compute class membership probabilities and then estimate the network structure for the elements and infer the network structure for the sets using those elements.

This procedure is model based, with only the data observed and with other parameters unobserved. The unobserved parameters inferred from Section 2 have statistical meanings. The procedure of the first step can be illustrated as a directed graph, as shown in Figure 1. Note that the data x can be affected by both $μ$ and $Σ$ . The focus of this paper is an approach based on the covariance matrix $Σ$ .

Regarding the second step, the procedure for estimating element networks is described in Section 2.3. The procedure of learning set networks is demonstrated in the block array below. According to the notation from Section 3, $P_{c, k k'}$ is the interaction event between set pair k and $k'$ for a given class c.

\begin{matrix} 1 & \dots & k & \dots & K \end{matrix}

\begin{matrix} (\begin{matrix} 1 & \dots & P_{c, 1 k} & \dots & P_{c, 1 K} \\ \dots & ⋱ & \dots & \dots & \dots \\ P_{c, k 1} & \dots & 1 & \dots & P_{c, k K} \\ \dots & \dots & \dots & ⋱ & \dots \\ P_{c, K 1} & \dots & P_{c, K k} & \dots & 1 \end{matrix}) \begin{matrix} 1 \\ ⋮ \\ k \\ ⋮ \\ K \end{matrix} \end{matrix}

5.1. Full conditional distribution

Let the following denote the parameters of interest in a vectorized fashion: $β_{c, j} = {[β_{c, j, 1}, β_{c, j, 2,}, \dots, β_{c, j, j - 1}]}^{⊤}, g_{c, j} = {[g_{c, j, 1}, g_{c, j, 2,}, \dots, g_{c, j, j - 1}]}^{⊤}, x_{c, j} = {[x_{c, 1, j}, x_{c, 2, j}, \dots, x_{c, N, j}]}^{⊤}$ . Let $Ψ_{c}$ be the set of parameters of interest.

We let $p (x_{c, j} | β_{c, j})$ represent the likelihood obtained from a multivariate normal distribution. In addition, $p (β_{c, j} | g_{c, j})$ is the prior distribution of $β_{c, j}, p (g_{c, j} | θ_{c, k, k'})$ is the prior distribution of $g_{c, j}$ , and $p (θ_{c, k, k'})$ is the prior distribution. Then, the joint posterior has the following form

\begin{array}{l} p (Ψ_{c}) \propto \prod_{j = 1}^{P - 1} p {(x_{c, j} | β_{c, j}) \cdot p (β_{c, j} | g_{c, j})} \cdot \prod_{k = 1}^{K - 1} {p (g_{c, j} | θ_{c, k, k'}) \times p (θ_{c, k, k'})} . \end{array}

Because all of the prior distributions are conjugate distributions, the full conditional distributions have the closed forms summarized as follows:

The full conditional distribution of $β_{j}$

\begin{array}{l} p (β_{c, j} | -) \propto p (x_{c, j} | β_{c, j}) \cdot p (β_{c, j} | g_{c, j}); \\ \propto \exp {- \frac{1}{2 σ_{c}^{2}} {[x_{c, j} - x_{c, - j} β_{c, j}]}^{⊤} [x_{c, j} - x_{c, - j} β_{c, j}]} \times \exp {- \frac{1}{2} β_{c, j}^{⊤} Σ_{c, g}^{- 1} β_{c, j}}; \\ \propto \exp {- \frac{1}{2} [β_{c, j}^{⊤} (σ^{- 2} x_{c, - j} x_{c, - j} + Σ_{c, g}^{- 1}) β_{c, j} - 2 β_{c, j}^{⊤} (σ_{c}^{2} x_{c, - j} x_{c, j})]} . \end{array}

One can show that $[β_{c, j} | -] \sim N (μ_{c}^{N}, V_{c}^{N}),$ where

\begin{array}{l} μ_{c}^{N} = {(σ_{c}^{- 2} x_{c, - j} x_{c, - j} + Σ_{g}^{- 1})}^{- 1} (σ_{c}^{2} x_{c, - j} x_{c, j}), \\ V_{c}^{N} = {(σ_{c}^{- 2} x_{c, - j} x_{c, - j} + Σ_{c, g}^{- 1})}^{- 1}, and \\ Σ_{c, g} = diag {v_{1}^{2} g_{c, j} + v_{0}^{2} (1 - g_{c, j})} . \end{array}

The full conditional distribution of $g_{j}$

\begin{array}{l} p (g_{c, j} | -) \propto p (β_{c, j} | g_{c, j}) \cdot p (g_{c, j} | θ_{c, k, k'}) \\ \propto | Σ_{c, g} |^{- \frac{1}{2}} \exp {- \frac{1}{2} β_{c, j}^{⊤} Σ_{c, g}^{- 1} β_{c, j}} \times {(θ_{c, k, k'})}^{g_{c, j}} {(1 - θ_{c, k, k'})}^{1 - g_{c, j}} \end{array}

and sample $g_{c, j}$ with probability $Pr (g_{c, j} = 1)$ .

The full conditional distribution of $θ_{c, k, k'}$

\begin{array}{l} p (θ_{c, k, k'} | -) \propto p (g_{c, j} | θ_{c, k, k'}) \cdot p (θ_{c, k, k'} | a, b) \\ \propto {(θ_{c, k, k'})}^{\sum_{j' > j} g_{c, j}} {(1 - θ_{c, k, k'})}^{\sum_{j' > j} (1 - g_{c, j})} \times {(θ_{c, k, k'})}^{a - 1} {(1 - θ_{c, k, k'})}^{b - 1} \\ \propto {(θ_{c, k, k'})}^{\sum_{j' > j} g_{c, j} + a - 1} {(1 - θ_{c, k, k'})}^{\sum_{j' > j} (1 - g_{c, j}) + b - 1} \\ \propto {(θ_{c, k, k'})}^{a^{N} - 1} {(1 - θ_{c, k, k'})}^{b^{N} - 1}, \end{array}

where $a^{N} = \sum_{j' > j} g_{c, j} + a$ and $b^{N} = \sum_{j' > j} (1 - g_{c, j}) + b$ .

We now present MCMC algorithm 1 for solving BMGGMs with unknown classes in greater detail.

Algorithm 1. MCMC algorithm for Bayesian multiple GGMs with unknown classes

Initialize $Ψ_{c} = {π_{c}, μ_{c}, β_{c}, g_{c}, Θ_{c}}$ for each class c.

Select the prior.

for each iteration do:

Evaluate the data class membership by $p (z_{c, i} | x_{c, i}) = \frac{π_{c} N (x_{i} | μ_{c}, Σ_{c})}{\sum_{c' = 1}^{C} N (x_{i} | μ_{c'}, Σ_{c'})}$ .

for each class $c = 1, \dots, C$ do

for each elements $j = 2, \dots, P$ do

Update $β_{c, j j'}$ and $g_{c, j j'}$

end for

for each set $k = 1, \dots, K$ do

Update $θ_{c, k, k'}$

end for

Compute the estimated precision matrix for each class using $Σ_{c}^{- 1} = L_{c}^{⊤} D_{c} L_{c}$

end for

Return: posterior samples of $g_{c}, z_{c}$ for each class c.

Our R package implementing BMGGMs with unknown classes is available on our GitHub repository.

5.2. Posterior inference

In this subsection, we represent two ways of estimating the posterior of the edge inclusion (i.e. $g_{c, j j'} = 1$ ) for each edge. Both approaches select the edges only marginally.

For the first practical approach, we obtained the posteriors of the edge inclusion after the burn-in. Then, we chose the edges whose posteriors of inclusion were larger than 0.5, which was the threshold used by Barbieri et al.³ and Peterson et al.¹⁸ We refer to this approach as BMGGM1.

The second approach is based on testing. We treat the edges selection as a hypothesis testing problem for a Bernoulli distribution in which ${\hat{g}}_{c, j j'} \sim N (g_{c, j j'}, S E ({\hat{g}}_{c, j j'}))$ , where $S E ({\hat{g}}_{c, j j'}) = \sqrt{\frac{{\hat{g}}_{c, j j'} (1 - {\hat{g}}_{c, j j'})}{N}}$ and N is the number of MCMC simulations. For large N, we assume

Z = \frac{g_{c, j j'} - {\hat{g}}_{c, j j'}}{S E ({\hat{g}}_{c, j j'})} \sim N (0, 1) .

Then, the approximate $1 - α$ confidence intervals for $g_{c, j j'}$ are

[{\hat{g}}_{c, j j'} - Z_{α / 2}^{*} \times S E ({\hat{g}}_{c, j j'}), {\hat{g}}_{c, j j'} + Z_{α / 2}^{*} \times S E ({\hat{g}}_{c, j j'})],

where

Z_{α / 2}

is a critical value and α is the significance level. Finally, we select the edges if the corresponding confidence interval does not contain 0. We refer to this approach as BMGGM2.

Furthermore, we may compare differential networks among different classes. Here, we define differential element networks if the value of $| g_{c, j j'} - g_{c', j j'} | = 1$ . This means the element pair $(j, j')$ exists in the network of either class c or class $c'$ but not both. Similarly, we test the following null hypothesis against the given alternative

H_{0} : g_{c, j j'} = g_{c', j j'} versus H_{a} : g_{c, j j'} \neq g_{c', j j'}

Thus, the test statistic is

Z = \frac{{\hat{g}}_{c, j j'} - {\hat{g}}_{c', j j'}}{\sqrt{{\hat{g}}_{j j'} (1 - {\hat{g}}_{j j'}) (\frac{1}{N} + \frac{1}{N})}}

where

{\hat{g}}_{j j'} = \frac{N {\hat{g}}_{c, j j'} + N {\hat{g}}_{c', j j'}}{N + N} = \frac{1}{2} ({\hat{g}}_{c, j j'} + {\hat{g}}_{c', j j'})

Then, the approximate $1 - α$ confidence intervals for $g_{c, j j'}$ are

[{\hat{g}}_{c, j j'} - {\hat{g}}_{c', j j'} - Z_{α / 2}^{*} \times \sqrt{{\hat{g}}_{j j'} (1 - {\hat{g}}_{j j'}) (\frac{1}{N} + \frac{1}{N})}, {\hat{g}}_{c, j j'} - {\hat{g}}_{c', j j'} + Z_{α / 2}^{*} \times \sqrt{{\hat{g}}_{j j'} (1 - {\hat{g}}_{j j'}) (\frac{1}{N} + \frac{1}{N})}]

Likewise, we claim that the edge is differential if the corresponding confidence interval does not contain 0. A similar idea can be extended to more than two differential networks by using multiple comparisons.

In the later experiments, we found that the first approach (BMGGM1), with a fixed threshold of 0.5, resulted in better performance than the second approach (BMGGM2) did. We will focus on the first approach throughout this paper unless otherwise specified.

6. Simulation

In this section, we present two simulation studies to evaluate the performance of our BMGGMs. In the first experiment, we estimated all of the parameters of interest and infer GGMs, assuming that the class labels of the observations are given. We also compared our BMGGM approach to the existing methods. The second experiment was similar to the first one, except that we assumed that the class labels were unknown. We evaluated the model’s ability to make predictions by learning differential graphs across classes.

6.1. Simulation study for known classes

6.1.1. Simulation setting

In this subsection, we run a simulation to evaluate the performance of our BMGGMs when the class is given. We start by drawing n = 100 samples independently and identically from multivariate Gaussian distribution $N (0, Σ_{c}^{- 1})$ . Because the class is given, we simplify the notation by letting $Σ = Σ_{c}$ .

The inverse covariance matrix $Σ$ depends on the corresponding network structure. To the form adjacency matrix, we created a matrix with ones on its diagonal and with zeros on entries not corresponding to network edges. To simulate the adjacency matrix for both elements and sets, we took the following steps:

Step S.1: We first created K = 6 set networks based on one of the network types (an AR(2) network, chain network, or scale-free network). Each set network had p_k = 10 elements. Thus, we could have corresponding adjacency matrices on the diagonal block of $Σ$ .

Step S.2: For the off-diagonal block of $Σ$ , we set the entries to 0 except for $g_{4, 17} = 1$ , where the 4th and 17th elements belong to sets 1 and 2, respectively. In this way, we connected sets 4 and 17 based on assumption A1 from Section 3. That is, $g_{4, 17}^{(1, 2)} = 1$ because element 4 interacts with element 17 in set pairs 1 and 2 for a given class.

For the elements, we considered three types of simulated networks – AR(2), chain, and scale-free networks-as shown in Figure 2. We explain how to generate these networks in detail:

AR(2) network: This is also called a second-order autoregression model. For a given set, the corresponding precision matrix $Ω$ is set with entries $ω_{j, j} = 1,$ for $j = 1, \dots, P$ ; $ω_{j, j + 1} = ω_{j + 1, j} = 0.5$ , for $j = 1, \dots, P - 1$ ; and .. for $j = 1, \dots, P - 2.$

Chain network: This network structure corresponds to a band matrix, which can be created by adding nonzero entries to a diagonal matrix. It resembles a tridiagonal precision matrix.⁸

Scale-free network: We generate set networks using the Barabasi-Albert algorithm,² each with a power-law degree distribution. That is, for a given set, the degree distribution is

p (k) \propto k^{- α},

where α is some prefixed parameter. Power-law degree distributions can mimic the network structure of biological data⁴ and are usually more difficult to learn than the other types of network structures.¹⁷

To create a symmetric and positive-definite covariance matrix $Σ$ , we first replaced the adjacency matrix with non-zero values from a uniform distribution on $[- 1, - 0.5] \cup [0.5, 1]$ . Then, we proceed with the Steps C.1–C.4 used by Danaher et al.,⁶ with which we divide each off-diagonal element in its row and average the matrix with its transpose:

Step C.1: Convert the network structure to the corresponding adjacency matrix (i.e. a (0, 1) matrix with zeros on its diagonal).

Step C.2: Replace the first entries with other nonzero entries $\sim uniform (0.5, 1)$ .

Step C.3: Make the matrix positive definite by dividing each off-diagonal entry by the sum of the absolute values of the off-diagonal entries in its row; then, average the matrix with its transpose.

Step C.4: Each observation is drawn independently and identically from the multivariate Gaussian distribution $N (0, Σ^{- 1})$ .

6.2. Evaluation metrics

Following the standard practice for the Gaussian graphical model, we assess the performance of the Gaussian graphical model by focusing on the estimation of the adjacency matrix or network structure for each set ks, (where $k = 1, \dots, 6$ ). The key measure of the model should reveal predictive power in obtaining an accurate sparse matrix, with the ability to detect not only the edges but also non-edges.

Several metrics are worth assessing: the true positive rate (TPR), true negative rate (TNR), and accuracy (ACC). To calculate these metrics, we first define the false positive (FP), true positive (TP), false negative (FN), and true negative (TN) for the edge status of a pair of nodes $g_{j, j'}$

\begin{array}{l} {TP}^{k} = \sum_{1 < j < j' < P} I {g_{j, j'}^{k} = 1, {\hat{g}}_{j, j'}^{k} = 1}; \\ {FP}^{k} = \sum_{1 < j < j' < P} I {g_{j, j'}^{k} = 0, {\hat{g}}_{j, j'}^{k} = 1}; \\ {TN}^{k} = \sum_{1 < j < j' < P} I {g_{j, j'}^{k} = 0, {\hat{g}}_{j, j'}^{k} = 0}; \\ {FN}^{k} = \sum_{1 < j < j' < P} I {g_{j, j'}^{k} = 1, {\hat{g}}_{j, j'}^{k} = 0} . \end{array}

Accordingly, we can calculate the TPR, TNR, and ACC for each higher level variable and average them

\begin{array}{l} TPR = \frac{1}{K} \sum_{k = 1}^{K} (\frac{{TP}^{k}}{{TP}^{k} + {FN}^{k}}); \\ TNR = \frac{1}{K} \sum_{k = 1}^{K} (\frac{{TN}^{k}}{{FP}^{k} + {TN}^{k}}); \\ ACC = \frac{1}{K} \sum_{k = 1}^{K} (\frac{{TP}^{k} + {TN}^{k}}{{TP}^{k} + {FP}^{k} + {FN}^{k} + {TN}^{k}}) . \end{array}

Because the inverse covariance matrix is rather sparse, we have far more non-edges than edges for a given graph. We also introduce precision and recall, which are widely used in information retrieval and binary classification. We can define them as follows

\begin{array}{l} Precision = \frac{1}{K} \sum_{k = 1}^{K} (\frac{{TP}^{k}}{{TP}^{k} + {FP}^{k}}); \\ Recall = \frac{1}{K} \sum_{k = 1}^{K} (\frac{{TP}^{k}}{{TP}^{k} + {FN}^{k}}) . \end{array}

Precision measures how many selected items are relevant, and recall measures how many relevant items are selected.

6.3. Simulation results

We implemented the Bayesian MCMC procedure described in Section 5 to obtain samples from the posterior probabilities. To make inferences, we ran our Gibbs sampler for 10,000 iterations, with another 10,000 as burn-in. We compared our method with the other three methods using the same dataset. First, we applied the graphical lasso by Friedman et al.⁹ (referred as Glasso). Then, we apply the group graphical lasso (referred to GGL) and fused graphical lasso (referred to FGL) proposed by Danaher et al.⁶ All of the tuning parameters were selected to give the lowest average cross-validation error.

We began by assessing whether the MCMC procedure could converge to the stationary distribution. For the scale-free network case, Figure 3 of the Supplementary Materials demonstrates the trace plots of the number of edges for each set network, which indicates good mixing rates of Markov chains. To estimate the adjacency matrix, we obtained the posterior probability of edge inclusion by calculating the average for the MCMC samples of $g_{j, j'}$ . The heat maps for the posterior probability of edge inclusion are displayed in Figures 4 to 6 of the Supplementary Materials, with their patterns showing good recovery of the true network structure. Figure 7 of the Supplementary Materials shows the receiver operator curves (ROC) for each set network along with the TPR and FPR, as described in Section 6. The MCMC trace plots, heat maps, and ROC curves for the chain network case are displayed in Figures 8 to 12 of the Supplementary Materials, and those for AR(2) network case are also displayed in Figures 12 to 17 of the Supplementary Materials.

We compared our method with the other three methods – Glasso, GGL, and FGL on the same data set. The results of the estimated network structure are summarized in Table 1. We assessed the accuracy of the estimated network structure via precision, recall, and ACC together, instead of ACC alone, because the graphs we obtained are rather sparse. We averaged these results over six set networks to make comparisons with other methods.

Table 1.

Comparison of five methods in terms of four measures with standard error (SE) over 100 simulated runs under scale-free, AR(1), and AR(2) network case.

Network	Method	TNR	Precision	Recall	ACC
	GLasso	0.45(0.07)	0.31(0.03)	1.00(0.00)	0.56(0.06)
Scale	FGL	0.98(0.02)	0.82(0.12)	0.52(0.00)	0.88(0.05)
-Free	GGL	0.92(0.03)	0.77(0.08)	1.00(0.00)	0.93(0.02)
	BMGGM1	0.98 (0.01)	0.92(0.07)	0.78(0.10)	0.94(0.02)
	BMGGM2	0.79 (0.01)	0.57(0.05)	0.98(0.01)	0.82(0.05)
	FGL	1.00(0.10)	NaN	0.00(0.03)	0.89(0.06)
AR(1)	GGL	1.00(0.04)	1.00(0.04)	0.15(0.04)	0.90(0.09)
	GLasso	0.62(0.03)	0.20(0.10)	1.00(0.01)	0.67(0.08)
	BMGGM1	0.90(0.03)	0.50(0.03)	0.50(0.10)	0.85(0.09)
	BMGGM2	0.89(0.01)	0.48(0.10)	0.52(0.07)	0.84(0.08)
AR(2)	FGL	1.00(0.03)	1.00(0.03)	0.88(0.02)	0.95(0.01)
	GGL	1.00(0.10)	1.00(0.09)	0.92(0.06)	0.97(0.08)
	GLasso	0.54(0.11)	0.58(0.07)	1.00(0.10)	0.72(0.07)
	BMGGM1	0.98(0.08)	0.97(0.10)	1.00(0.07)	0.98(0.10)
	BMGGM2	0.88(0.03)	0.83(0.04)	0.99(0.09)	0.92(0.09)

GLasso: graphical lasso⁹; GGL: group graphical lasso⁶; FGL: fused graphical lasso⁶; BMGGM1: our proposed method with thresholding; BMGGM2: our proposed method with testing.

The results suggest that our BMGGMs had the best accuracy overall. GGL performed slightly better than FGL, which was implied by Danaher et al.⁶ Glasso has the lowest accuracy because it learns each graph individually without considering the connections between the set networks. With these simulation settings, our BMGGMs will overestimate the edges, as compared to joint glasso, which led to relatively higher TPR and recall, and relatively lower TNR and precision.

Lastly, our BMGGMs could identify the connections between the sets, as we can clearly see in Figure 3, where there is an edge across sets 1 and 2.

6.4. Simulation study for unknown classes

6.4.1. Simulation setting

In this subsection, we run a simulation to evaluate the performance of our algorithm when the class is not given. For each class, we started by drawing samples independently and identically from the multivariate Gaussian distribution $N (0, Σ_{c}^{- 1})$ , where c = 1, 2. The sample size here is 1000, with 200 cases and 800 controls.

For the covariance matrices, we considered K = 10 sets and p_k = 10 elements each. $Σ_{1}$ was created in the same way as described in Section 6.1. To form $Σ_{2}$ , we removed the edges randomly in $Σ_{1}$ (i.e. setting the non-zero entries in the corresponding adjacency matrix to 0) and added new edges randomly by replacing zero entries with non-zero values. We set the proportion of differential edges for these two graphs at 20%.

6.4.2. Evaluation metrics

The evaluation criteria were the same as those used in Section 6.2. Likewise, we redefined the FP, TP, FN, TN, and ACC for the class of observations

\begin{array}{l} TP = \sum_{i = 1}^{n} I {c = 1, \hat{c} = 1}, FP = \sum_{i = 1}^{n} I {c = 2, \hat{c} = 1}; \end{array}

\begin{array}{l} TN = \sum_{i = 1}^{n} I {c = 2, \hat{c} = 2}, FN = \sum_{i = 1}^{n} I {c = 1, \hat{c} = 2}; \\ ACC = \frac{1}{n} \sum_{i = 1}^{n} (\frac{TP + TN}{TP + FP + FN + TN}) . \end{array}

6.4.3. Performance when the unknown class size is 2

Because this was a clustering problem, a label-switching problem may arise in each iteration. We solved this by adding the constraint $n_{1} < n_{2}$ ; that is, $π_{1} < π_{2}$ . Figure 4 displays how the algorithm converges as posterior predictions of the numbers of cases and controls are become more spread out at both ends. With this, the algorithm was able to achieve 81.2% accuracy over 100 simulated runs.

6.4.4. Results for unknown classes with different sizes

We also conducted an additional simulation when C was 3 with a scale-free network. We calculated ACC and the standard deviation (SD) of ACC obtained from 100 simulations. The rest of the settings were same as those used in Section 6.4.1. The experiments showed that the performance of the BMGGM in terms of ACC decreased as C increased. Table 2 represents the performance of BMGGM when C varied from 1 to 3 with the scale-free network. As expected, the value of ACC decreases as C increases from 1 to 3. In addition, we observed that the BMGGM improved, as the sample size among classes became more unbalanced. A BMGGM does not break down even if the sample is small, such as 50. When C = 2 and the sample sizes of the two classes are (n₁, $n_{2})$ =(50, 950), the ACC values is 78.4%. However, when C = 3 and the sample size of three classes is set as (n₁, n₂, n₃)=(800, 150, 50), the ACC values is 64.1%.

Table 2.

The performance of BMGGM1 in terms of accuracy (ACC) and its the standard deviation (SD) of ACC obtained from 100 simulations when C varies from 1 to 3.

C	Network	Sample size of each class	ACC (SD)
1	Scale-Free	n = 1000	$1.00 (0)$
2	Scale-Free	$n_{1} = 200, n_{2} = 800$	$0.80 (0.01)$
		$n_{1} = 100, n_{2} = 900$	$0.79 (0.01)$
		$n_{1} = 50, n_{2} = 950$	$0.78 (0.02)$
3	Scale-Free	$n_{1} = 800, n_{2} = 150, n_{3} = 50$	$0.64 (0.08)$
		$n_{1} = 700, n_{2} = 200, n_{3} = 100$	$0.60 (0.07)$

7. Application

In this section, we apply our BMGGMs to reveal the dependence structure for the breast cancer gene expression data, in which genes can be viewed as elements and sets can be viewed as sets.

According to the literature, American white women are slightly more likely to develop breast cancer than African American women are.⁵ Therefore, our goal was to compare differential networks among racial groups to better understand the genetic differences. These differential networks included not only the pathway network but also the gene network with each pathway. Specifically, we considered a gene–gene pair to be differential if the true value of $| g_{1, j, j'} - g_{2, j, j'} | = 1$ in our example, which reflects that edge $(j, j')$ is included among either white women or non-white women but not both. This means the edge between a gene–gene pair is different in white women than in non-white women.

The human breast cancer data set was collected from the University of Texas M.D. Anderson Cancer Center²¹ and contains 22,283 gene-expression measurements across 176 white patients and 102 non-white patients. Furthermore, the genes were mapped into 1320 pathways using the canonical pathways (CP) from the Molecular Signatures Database (MsigDB), with the number of genes within each pathway ranging from 4 to 778. We applied our method to the top three most significant pathways (P711, P717, and P956) expressed between white and non-white women breast cancer patients. In this data set, we had n = 278 samples, c = 2 classes, and K = 3 sets. The numbers of genes in these three pathways (P711: REACTOME ORC1 REMOVAL FROM CHROMATIN, P717: REACTOME SIGNALING BY ERBB2, and P956: REACTOME DOWNSTREAM SIGNAL TRANSDUCTION) were $p_{P 956} = 59, p_{P 711} = 21$ , and $p_{P 717} = 17$ , respectively. The dimension of the data matrix was 97.

We normalized all of the genes with the mean of 0 and standard deviation of 1 for each class. Because goals were to generate meaningful patterns and create a hypothesis, we set the hyperparameters to yield a sparse pathway network and a gene network. We chose the beta hyperparameters for $θ_{c, k, k'}$ of $a = 1, b = 2$ , regardless of whether $k \neq k'$ or $k \neq k'$ . We set the gamma hyperparameters e = 10, 000 and f = 1000 to further enforce the sparsity.

In Figures 5 and 6, we present the gene networks for white women and non-white women, in which the colors represent the pathways (because the genes could be grouped into pathways). The results show that the element structures for these two groups were very different. Furthermore, the graph of white women had more edges (144 edges) than that of non-white women (114 edges), and only 51 edges overlapped. In comparing the vertex degrees (i.e. the number of nearest neighbors for a vertex) of the two networks, we found that the white women network had more genes with large degrees (see Figure 18 of the Supplementary Materials).

At the pathway level, we identified pathway P956 as being conditionally dependent on pathway P717 among white women, while P956 and P711 were conditionally independent. The adjacency matrix for these two groups are summarized below

\begin{array}{l} Θ_{white women} = (\begin{matrix} 1.000 & 1.000 & 1.000 \\ 1.000 & 1.000 & 1.000 \\ 1.000 & 1.000 & 1.000 \end{matrix}) \\ Θ_{non white women} = (\begin{matrix} 1.000 & 1.000 & 1.000 \\ 1.000 & 1.000 & 0.006 \\ 1.000 & 0.006 & 1.000 \end{matrix}) \end{array}

Analysis of the gene and pathway networks reveals resulted consistent with the breast cancer genomics,¹² suggesting that American white women are slightly more likely to develop breast cancer than African American women are. Overall, the novel application of our BMGGMs to understand breast cancer networks yielded results consistent with the known literature and has identified potential biomarkers and pathways for future research.

On the other hand, we are also interested in applying BMGGMs when the class is unknown. We treated race (i.e. white and non-white women) as a class and evaluated its performance in terms of prediction rate.

We applied the same settings as the previous section. That is, we chose the P956, P711, and P717 pathways. We normalized all of the genes with means of 0 and 1 standard deviation of 1. The hyperparameters were also the same to yield sparse pathway and gene networks. We set the beta hyperparameters on $θ_{c, k, k'}$ to $a = 1, b = 2$ , regardless of $k \neq k'$ or $k \neq k'$ and used the gamma hyperparameters e = 10,000 and f = 1000, to further enforce the sparsity.

In addition, we forced a constraint ( $n_{1} < n_{2}$ , where n₁ is the number of cases and n₂ is the number of controls) to prevent the label-switching problem. The proposed method achieved a prediction rate of 64.7% because all of the samples had breast cancer. We discovered a very unbalanced clustering result. The prediction accuracy and clustering may have been affected by the non-zero off-diagonal estimation.

8. Discussion

We have proposed the BMGGMs as a framework for learning the network structure of multiple groups that could be possibly connected to each other. Given a data set with a hierarchical structure, these models can uncover the network structure for sets and elements simultaneously. Furthermore, our approach can be used to evaluate the Bayes factor to decide the class membership for each class. We designed and employed a simple Gibbs sampling scheme to solve the BMGGMs.

Our method is comparable to the existing methods in terms of the accuracy of adjacency matrix estimation. Our method also provides good statistical interpretations through posterior probabilities for each parameter because our generative model is more explainable and interpretable than deterministic models. We defined the model by explaining the generative mechanism in a top-down fashion.

In this paper, we made several novel contributions: we (1) developed BMGGMs that can reveal set and element networks simultaneously; (2) established a Bayesian factor approach to evaluate class membership for observations from GGMs; and (3) estimated the models by using a Bayesian rather than a frequentist approach. This Bayesian formulation is fully model based and thus has nice probabilistic interpretations. Clustering the data for the GGMs is very difficult, and only a few works have focused on the relevant problems so far.

There are several potentials for future work. For instance, when it comes to clustering the observations, one could discuss the convergence properties of the Bayesian method. Another possibility is to extend our model to allow discrete data.

We focused on network estimation. However, it is difficult to estimate the non-zero off-diagonal entries of precision matrix accurately. Because these entries play an important role in making statistical inferences on test hypotheses regarding networks as well as predictions, it will be worthwhile to develop statistical procedures to estimate them accurately in order to we can conduct better statistical inferences regarding networks.

Notably we investigated the performance of our BMGGM as the cluster size increased. As expected, the performance decreased as the cluster size increased. However, we also observed that the performance improved as the sample size increased. Also, the more unbalanced the sample size among classes, the better the BMGGMs performance. If some cluster has extremely small size such as 1, the BMGGM might not able to detect it. Hence, it is important to study the relationship between the optimal sample size and the class size. This will require intensive theoretical derivation and also simulation studies for future research.

Supplemental Material

sj-pdf-1-smm-10.1177_09622802211022405 - Supplemental material for Bayesian multiple Gaussian graphical models for multilevel variables from unknown classes

Supplemental material, sj-pdf-1-smm-10.1177_09622802211022405 for Bayesian multiple Gaussian graphical models for multilevel variables from unknown classes by Jiali Lin and Inyoung Kim in Statistical Methods in Medical Research

Footnotes

Acknowledgements

We are grateful to the associate editor and reviewer for their valuable suggestion and constructive input.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Inyoung Kim

Supplemental Materials

Supplementary materials are available in a separate file as part of the article. An R package implementing BMGGMs with unknown classes is available on the author’s GitHub repository.

References

Banerjee

O, El

Ghaoui

d’Aspremont

Model selection through sparse maximum likelihood estimation for multivariate gaussian or binary data. J Mach Learn Res 9(Mar): 485–516, 2008.

Barabási

A-L

Albert

Emergence of scaling in random networks. Science 1999; 286: 509–512.

Barbieri

Berger

, et al. Optimal predictive model selection. Ann Stat 2004; 32: 870–897.

Chen

Sharp

BM.

Content-rich biological network constructed by mining Pubmed abstracts.

BMC Bioinform 2004; 5: 147.

Chlebowski

Chen

Anderson

, et al. Ethnicity and breast cancer: factors influencing differences in incidence and outcome. J Natl Cancer Inst 2005; 97: 439–448.

Danaher

Wang

Witten

DM.

The joint graphical lasso for inverse covariance estimation across multiple classes.

J R Stat Soc: Series B (Statistical Methodology) 2014; 76: 373–397.

Dempster

AP.

Covariance selection. Biometrics 1972; 28: 157–175.

Fan

Variable selection via nonconcave penalized likelihood and its oracle properties. J Am Stat Assoc 2001; 96: 1348–1360.

Friedman

Hastie

Tibshirani

Sparse inverse covariance estimation with the graphical lasso.

Biostatistics 2008; 9: 432–441.

10.

George

McCulloch

RE.

Variable selection via Gibbs sampling. J Am Stat Assoc 1993; 88: 881–889.

11.

Guo

Levina

Michailidis

, et al. Joint estimation of multiple graphical models. Biometrika 2011; 98: 1–15.

12.

Keenan

Moy

Mroz

, et al. Comparison of the genomic landscape between primary breast cancer in African American versus white women and the association of racial differences with tumor recurrence. J Clin Oncol 2015; 33: 3621.

13.

Kim

Liu

Zhao

Bayesian methods for predicting interacting protein pairs using domain information.

Biometrics 2007; 63: 824–833.

14.

Kuo

Mallick

Variable selection for regression models. Sankhyā 1998; 60: 65–81.

15.

Marlin

Murphy

KP.

Sparse gaussian graphical models with unknown block structure. In: Proceedings of the 26th annual international conference on machine learning. New York: ACM, 2009, pp.705–712.

16.

Meinshausen

Bühlmann

High-dimensional graphs and variable selection with the lasso.

Ann Stat 2006; 34: 1436–1462.

17.

Peng

Wang

Zhou

, et al. Partial correlation estimation by joint sparse regression models. J Am Stat Assoc 2009; 104: 735–746.

18.

Peterson

Stingo

Vannucci

Bayesian inference of multiple gaussian graphical models. J Am Stat Assoc 2015; 110: 159–174.

19.

Rothman

Bickel

Levina

, et al. Sparse permutation invariant covariance estimation. Electron J Stat 2008; 2: 494–515.

20.

Rue

Held

Gaussian Markov random fields: theory and applications. Boca Raton, FL: CRC Press, 2005.

21.

Shi

Campbell

Jones

, et al. The microarray quality control (MAQC)-ii study of common practices for the development and validation of microarray-based predictive models. Nature Biotechnol 2010; 28: 827.

22.

Wang

, et al. Bayesian graphical lasso models and efficient posterior computation. Bayesian Anal 2012; 7: 867–886.

23.

Yuan

Lin

Model selection and estimation in the gaussian graphical model. Biometrika 2007; 94: 19–35.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

5.38 MB