A Stochastic Block Ising Model for Multi-Layer Networks with Inter-Layer Dependence

Abstract

Community detection has attracted tremendous interests in network analysis, which aims at finding group of nodes with similar characteristics. Various detection methods have been developed to detect homogeneous communities in multi-layer networks, where inter-layer dependence is a widely acknowledged but severely under-investigated issue. In this paper, we propose a novel stochastic block Ising model (SBIM) to incorporate the inter-layer dependence to help with community detection in multi-layer networks. The community structure is modeled by the stochastic block model (SBM) and the inter-layer dependence is incorporated via the popular Ising model. Furthermore, we develop an efficient variational EM algorithm to tackle the resultant optimization task and establish the asymptotic consistency of the proposed method. Extensive simulated examples and a real example on gene co-expression multi-layer network data are also provided to demonstrate the advantage of the proposed method.

1. Introduction

Community is one of the most widely observed structures in network data, where objects in the same community tend to share some common characteristics and similar connectivity. In a typical gene interaction network, genes sharing similar functions or participating in a common cellular process may form gene communities. For instances, BRCA1, BRCA2, and ATM are widely known DNA damaged response (DDR) genes, as they are deeply implicated in the function of repairing damaged DNA (Lang et al., 2019); epidermal growth factor receptor (EGFR), rearranged during transfection (RET), and mesenchymal epithelial transition factor (MET) have similar functions in non-small cell lung cancer (NSCLC), as EGFR Exon T790M, RET rearrangement and MET exon 14 skipping are well-known NSCLC driver mutations. In the past two decades, detecting community structures have attracted tremendous interests in the literature of network analysis (Lei & Rinaldo, 2015; Newman, 2006; Zhang et al., 2022).

Although most existing community detection methods focus on single-layer networks (Holland et al., 1983; Newman & Girvan, 2004; Newman, 2006; Rohe et al., 2011), multi-layer networks have been popularly encountered in biomedical data. For example, the gene co-expression networks of rhesus monkeys at different developmental stages (Lei et al., 2020) and the brain fMRI imaging networks from schizophrenia and normal subjects (Yuan & Qu, 2021). Yet, methods for multi-layer networks are still in their infancy, and most of the existing methods are extended from those on single-layer networks by assuming the layers are independent from each other. Examples include spectral clustering on the mean of adjacency matrices (Kumar et al., 2011; Paul & Chen, 2020; Tang et al., 2009), maximum likelihood estimation method (Han et al., 2015), least-square estimation methods (LSE; Lei et al., 2020), and the restricted multi-layer SBM model (Paul & Chen, 2016). Although success has been widely reported for discovering communities in multi-layer networks, these methods largely ignore the dependence structure among different layers, which may lead to a sub-optimal performance in community detection.

Inter-layer dependence is a widely acknowledged but severely under-investigated issue in multi-layer network data. For example, some time-varying biological network data (Lei et al., 2020; Zhang et al., 2020) exhibit strong inter-layer dependence among different layers, especially the adjacent layers. Spatial correlation may also exist in multi-layer network constructed from different brain subregions (Lei et al., 2020). Modeling inter-layer dependence in multi-layer network is analogous to modeling covariance structure in dependent data (Bester et al., 2011; Yuan & Qu, 2021), which shall lead to a more appropriate network model and benefit the analysis of multi-layer biological network data.

One of the key challenges for the multi-layer network model is that incorporating inter-layer dependence may substantially increase model complexity as well as computational expenses (Panagiotelis et al., 2012; Sklar, 1959). Only a few attempts have been made in the literature to characterize the inter-layer dependence. For instance, linked matrix factorization (LMF; Tang et al., 2009) approximates each layer by matrix factorization with a layer-specific factor and a common factor across all layers. Co-regularized spectral clustering (COR; Kumar et al., 2011) uses a regularization term to force the eigenvectors of each layer to align with some common matrix. However, these methods are mostly algorithm-based, and characterizing the inter-layer dependence structure with just a common matrix is rather ad hoc and lacks theoretical foundation.

In this paper, we propose a novel SBIM model to incorporate the inter-layer dependence for community detection in multi-layer networks. It models the homogeneous community structure with the SBM model, and employs the popular Ising model (Cheng et al., 2014; Ravikumar et al., 2010) to characterize the inter-layer dependence. The resultant pseudo-likelihood function is then optimized via a variational expectation maximization (EM) algorithm that optimizes a separable lower bound of the objective function. The asymptotic estimation consistency of the SBIM model is established, and its superior numerical performance is also compared against some popular multi-layer community detection methods on extensive simulated examples and a real application on a gene co-expression multi-layer network data.

The rest of the paper is organized as follows. Section 2 introduces the proposed SBIM model. Section 3 develops an efficient computing algorithm based on a variational EM scheme to tackle the resultant optimization task. Section 4 establishes the asymptotic estimation consistency of the SBIM model. In Section 5, we examine the performance of the proposed method with various simulated networks and compare it against some popular competitors in literature. A real application on a gene co-expression multi-layer network data is provided to further demonstrate the strength of the proposed SBIM model. A brief summary is provided in Section 6. Supporting information is also provided for computing details and all technical proofs.

2. Stochastic Block Ising Model

Consider a multi-layer network consisting of n common nodes, whose adjacency matrices are denoted as ⁠, where denotes the mth layer, and if there is an edge between nodes i and j in the mth layer and otherwise. Suppose that these n nodes possess a homogeneous community structure with K communities across all layers, and the community membership matrix is denoted as with ⁠, where if node i belongs to the kth community and 0 otherwise. Further, each column of contains at least one 1, indicating there is no empty community in the multi-layer network.

To model the joint likelihood of ⁠, we turn to model and separately. First, each row of is independently modeled with a multinomial distribution

where with for each and ⁠. Therefore, each row of will only have exactly one 1, corresponding to the non-overlapped community structure. Thus, the log-likelihood of can be written as

Next, we propose a SBIM model to characterize the joint distribution of given ⁠. Specifically, 's are mutually independent given ⁠, and the log-likelihood of given becomes

where is an order-4 tensor with being a symmetric matrix, ⁠, and is the partition function ensuring to be a proper probability density function.

The proposed SBIM model has a number of compelling features. On the one hand, depends on nodes i and j only through their community memberships, echoing the popular SBM model (Holland et al., 1983). On the other hand, can efficiently capture the dependence among 's for each node pair ⁠, due to the properties of the Ising model (Cheng et al., 2014). Given ⁠, provides a concise quantification of the conditional dependence between and given ⁠, where ⁠. More specifically, and are independent given all other 's if and only if ⁠, and the strength of their conditional dependence changes monotonically with the value of ⁠. Furthermore, if only when ⁠, the proposed SBIM model reduces down to the case with no inter-layer dependence (Han et al., 2015).

Let ⁠, then the joint log-likelihood of is

It can be difficult to directly optimize due to the intractable partition function in ⁠. A natural treatment is to approximate with a pseudo-likelihood function (Besag, 1974; Ravikumar et al., 2010)

where the second equality derives from the binary logistic regression model of over ⁠. Then, the joint pseudo log-likelihood of becomes

where ⁠. The pseudo log-likelihood of is then ⁠.

3. Variational EM

Maximizing yields the maximum pseudo-likelihood estimate ⁠. However, note that involves the summation over possible values of ⁠, and thus the traditional EM algorithm (Dempster et al., 1977) cannot be directly applied to optimize ⁠. As a remedy, we exploit a variational EM approach (Mariadassou et al., 2010), which can greatly facilitate computation by considering a lower bound of and restricting the distribution of to a factorized form.

To proceed, we consider a lower bound of that

where denotes the Kullback–Leibler (KL) divergence, ⁠, and can be any distribution of to approximate ⁠. Particularly, we set to be a completely factorized distribution that

where denotes the multinomial distribution with parameters satisfying for each i. Then, and for any ⁠. Denote as with for simplicity, then it follows from the definition of the KL divergence that

where is the entropy of ⁠, and

To optimize ⁠, we proceed to update and alternatively until convergence. Specifically, given at the tth step, we solve the following two sub-optimization tasks separately:

In (Equation 5), taking partial derivative with respect to yields that for each

Further, note that

With step size γ₀, is then updated as

The update of in (Equation 6) can be solved by the Lagrange multiplier. Let

where are the Lagrangian multipliers. Then, is the solution of the following fixed-point problem:

The computational details about the derivative of and are provided in the Supporting information.

It is important to remark that the variational EM algorithm is guaranteed to converge to a local maximum (Mariadassou et al., 2010). Denote the variational estimate as with ⁠. Then, for each ⁠, can be regarded as the estimated membership of node i. To prevent it from a degenerated optimum, it is helpful to warm start the algorithm with a proper initializer. A number of initialization strategies have been developed for the variational EM algorithm, such as the LSE (Lei et al., 2020) and the layer-wise maximum likelihood estimate (MLE; Han et al., 2015). In our numerical experiments, we use the estimates from MLE to initialize in the alternative updating algorithm, which has satisfactory performance. Furthermore, to initialize ⁠, we simply set all of its elements to be 0 as little prior knowledge about is available.

4. Theory

In this section, we establish some theoretical results to quantify the asymptotic behavior of the proposed method in terms of estimation of the true parameter and the posterior distribution of ⁠.

Let with indicating whether node i belongs to true community k. Let ⁠, and we assume the limits of exist for each ⁠, denoted as ⁠. Then, the multi-layer network is generated from the joint density function ⁠. The estimation error of and is measured by

where σ can be any permutation on [K], and for any with ⁠. Then, is a variant of KL divergence. Lemma 1 shows that the pseudo log-likelihood attains its maximum at ⁠, suggesting that it is an appropriate surrogate loss to facilitate estimation of and implying that ⁠. This leads to an interesting analogy to the popular Fisher consistency in the machine learning literature (Lin, 2004).

Lemma 1. For any ⁠, if ⁠, then we have

for any and ⁠, and the equality holds if and only if ⁠.

The following technical assumptions are made. We will write if and if ⁠.

Assumption 1. For any and ⁠, for some and there exists a positive sequence such that ⁠.

Assumption 2. For any ⁠, there exists an such that for any ⁠.

Assumption 1 assumes a lower bound for ⁠. It also admits sparse network structure as well as flexible dependence structure among different layers by allowing a diverging upper bound on ⁠. Assumption 2 guarantees that the community structure is distinguishable in that there exists some ⁠, whose rows are sufficiently different one from another.

Theorem 1. Suppose Assumptions 1 and 2 hold. If and for some ⁠, then for any ⁠, it holds true that

Theorem 1 establishes the asymptotic estimation consistency of ⁠. Note that is allowed to diverge at a fast rate, indicating strong dependence among different layers. In fact, Theorem 1 still holds true for unbalanced communities when tends to 0 not too fast. The technical proofs of this theorem and all other theorems are given in the Supporting information.

Next, the asymptotic consistency of the variational estimate can also be established.

Theorem 2. Suppose all assumptions in Theorem 1 hold. For any ⁠, it holds true that

Theorem 2 assures that the estimates from the variational EM algorithm is also asymptotically consistent. This theoretical result is particularly interesting, as the global optimizer of the joint pseudo-likelihood (3) is difficult to obtain, if not impossible. The variational estimate, on the other hand, provides a computationally feasible and theoretically justifiable estimation of the proposed SBIM model.

We now turn to establish some consistency results of the proposed method in terms of community detection under additional assumptions.

Assumption 3. For any ⁠, if ⁠, then

for any ⁠. Further, there exist ⁠, such that is not permutation invariant.

Note that it follows from Lemma 1 that is always positive if and ⁠. Thus, Assumption 3 further assures that it cannot vanish too fast.

Theorem 3. Suppose all assumptions in Theorem 1 and Assumption 3 hold. Let ⁠, then we have ⁠. Furthermore, let subject to and ⁠, then it also holds true that as ⁠.

Theorem 3 shows that the pseudo posterior distribution is consistent in the sense that its maximizer exactly recovers the true community structure as defined by ⁠. Even though is computationally intractable, Theorem 3 further assures that it can be well approximated by when is given. In the challenging case with unknown ⁠, it would be necessary to quantify the behavior of as in Assumption 3 to obtain similar consistency results as in Theorem 3 under additional assumptions.

5. Numerical Experiments

In this section, we examine the numeric performance of the proposed SBIM model in both synthetic and real networks, and compare it against five popular competitors in the literature, including spectral clustering on mean adjacency matrix (SCM; Paul & Chen, 2020), LSE of multi-layer block models (Lei et al., 2020), LMF Paul & Chen, 2020), co-regularized spectral clustering (COR; Paul & Chen, 2020), and layer-wise MLE of SBM (Han et al., 2015). The implementations for LSE, LMF, and COR are available at the authors' personal website, whereas SCM and MLE are implemented by ourselves.

The performance of each method is evaluated by the community detection error (Wang, 2010), which is defined as

where is the indicator function, and denote the true and estimated community labels, respectively. Clearly, measures the disagreement between and ⁠, and is invariant to the permutation of community labels.

5.1 Synthetic Examples

We consider three synthetic networks concerning various aspects of the competing methods.

Example 1. In this example, we vary the network size, the layer number, and the dependence strength. We set ⁠, ⁠, and ⁠. Then, for ⁠, where is the identity matrix and is the symmetric hollow matrix with for ⁠, and for ⁠. Let and ⁠, where can measure the strength of the dependence among different layers. Following the treatment in Ravikumar et al. (2010) and Cheng et al. (2014), Gibbs sampling is exploited to generate ⁠. Before is generated, we add a disturbance on to slightly blur the community structure, where is a symmetric matrix with for ⁠. The averaged community detection errors of the competing methods over 50 independent replications, as well as their standard errors are summarized in Table 1.

It is evident from Table 1 that SBIM outperforms the other five methods in Example 1. The performances of SCM, LMF, and COR appear to be less satisfactory, whereas MLE and LSE are more competitive but still substantially worse than SBIM. Almost all methods perform better as the network size grows. This is due to the fact that more information is available as the number of nodes increases. It is worth noting that compared to the cases of ⁠, the performances of SBIM, MLE, and LSE on cases of are worse. One possible reason is that the dependence relationship is difficult to be fully captured by LSE when the number of layers increases, leading to less accurate initializers for SBIM and MLE. It is also interesting to note that the performances of all methods become worse when μ gets larger or the dependence becomes stronger. Furthermore, the difference between SBIM and MLE when is smaller than the difference when μ is large, which demonstrates the importance of incorporating the dependence structure.

Example 2. In this example, we test the performances of the competing methods on networks with unbalanced communities. The data-generating scheme is similar to Example 1, except that will vary and other parameters are fixed. We set ⁠, ⁠, and for ⁠, where is the diagonal matrix with independently for ⁠, and is the symmetric hollow matrix with independently for and for ⁠. Let with controlling the unbalance degree. The averaged community detection errors of the competing methods over 50 independent replications, as well as their standard errors are summarized in Table 2.

It is clear from Table 2 that SBIM outperforms the other five competitors in all cases. It is evident that the community detection errors of almost all methods decreases as α₁ increases, or the communities become more balanced. It is also interesting to note that SCM, LMF, and COR are less sensitive to α₁ compared with the other three methods.

Example 3. In this example, we test the performances of the competing methods on sparse networks. The data-generating scheme is same as Example 2 except that we set and ⁠) for ⁠, where is the diagonal matrix with for ⁠. Let control the network sparsity. Then, the network densities, defined as the average of elements in over 50 independent replications, are 0.11, 0.06, 0.03, 0.01. It is clear that the networks are sparser as μ increases. The averaged community detection errors of the competing methods over 50 independent replications, as well as their standard errors are summarized in Table 3.

Table 3 shows that SBIM outperforms the competitors in all cases. It is clear that the community detection errors of all methods increase as μ increases or the networks become sparser. Yet, SBIM still produces superb numerical performance compared with all five competitors, where MLE appears to perform better than the other four.

5.2 Real Application

We now apply the proposed SBIM model to analyze a gene co-expression data, which describes different developmental stages of medial prefrontal cortex of rhesus monkeys and helps understand developmental brain disorders. The data were first analyzed in Bakken et al. (2016), and its community structure was also studied in Liu et al. (2018), Lei et al. (2020) and Lei and Lin (2022). We choose the six stages for the prenatal period, which are labeled as E40 to E120 indicating the number of embryonic days of age. Further, we focus on the neural projection guidance-enriched genes, which has been shown to be related to autism spectrum disorder by Liu et al. (2018). This results in a six-layer network with 171 nodes. The edges are constructed by thresholding the sample correlations among node pairs, which is a common practice in gene network pre-processing to help keep most likely biologically-significant relationships (Bellec et al., 2017; Buchanan et al., 2020; Larremore et al., 2013; Perkins & Langston, 2009; Stuart et al., 2003). Specifically, let be the mth correlation matrix of the continuous gene co-expression networks with ⁠. Then, the binary adjacent matrix is constructed by setting for some ⁠. We set ξ as the 80% percentile of all 's, and thus the sparsity of the resulting binary networks is 0.2.

Note that each layer in this six-layer network indicates different developmental stages. Thus, the former layers will have effects on the latter layers. Then, it is natural that different layers are correlated, especially the adjacent layers. To validate such an observation, we calculate the correlation coefficient for any two different layers by converting them as 171 × 171-dimensional vectors. As a comparison, we also generate six networks independently with 171 nodes and same edge density, and calculate its correlation coefficient for any two different layers. The results are plotted with heat maps in Figure 1.

It is clear from Figure 1 that the layer-wise correlation coefficient for the gene co-expression multi-layer network is substantially larger than that for the randomly generated independent multi-layer network, which indicates that different layers of the gene co-expression multi-layer network are correlated. More interestingly, the layer-wise correlation coefficient between adjacent layers of the gene co-expression multi-layer network is also larger than others from the left panel of Figure 1, which supports our arguments.

Following the suggestion of Lei et al. (2020), we apply all competing methods to this gene co-expression multi-layer network data with ⁠. Note that the true community structure of this gene network is not available. Thus, we use the averaged internal density (Xu et al., 2022; Yang & Leskovec, 2012) to evaluate the performances of all methods. Specifically, we define

where is the estimated community assignment vector, is the size of kth estimated community, contains all nodes pair in the kth estimated community, and is the edge density and introduced as baseline. The estimated averaged internal densities of all methods are reported in Table 4. It is clear that the proposed SBIM model is the best performer with the most improvement being 8.5%, suggesting that the proposed framework is effective in incorporating the inter-layer dependence to improve the estimated community structure.

6. Discussion

In this paper, we propose a novel SBIM model to detect communities in multi-layer networks with inter-layer dependence. Whereas the community structure across different layers is modeled via the SBM model, the inter-layer dependence among layers is incorporated via the Ising model. We also develop an efficient variational EM algorithm to tackle the resultant optimization task, and establish the asymptotic estimation consistency of the developed model. Its superior numerical performance is demonstrated in various synthetic networks and a real-life gene co-expression multi-layer network data. As further extensions, it would be desirable to establish theoretical results on the exact recovery of with the variational EM estimate in more realistic scenarios without assuming knowledge of true parameters.

Data Availability Statement

The data that support the findings in this paper are available in the website at https://portland-my.sharepoint.com/:u:/g/personal/jingnzhan2-c_my_cityu_edu_hk/EXUx6ii_tOxGq0N6QnL8GtYBsEAXE4Fdnz4DU9w0hub2hQ?e=ybPXOe.

Acknowledgments

The authors are grateful to reviewers, associate editor, and editor for their insightful comments and suggestions which have improved the manuscript significantly. Jingnan Zhang's research is supported in part by “USTC Research Funds of the Double First-Class Initiative” YD2040002020, Junhui Wang's research is supported in part by HK RGC grants GRF-11304520, GRF-11301521, and GRF-11311022, and CUHK Startup Grant 4937091.

References