Inference of Edge Correlations in Multilayer Networks

TL;DR

This paper introduces a correlated multilayer stochastic block model that captures edge dependencies across layers, improving community detection and edge prediction in complex multilayer networks.

Contribution

It relaxes the independence assumption in multilayer SBMs by incorporating edge correlations, providing maximum-likelihood estimation and a new layer correlation measure.

Findings

Edge correlation modeling enhances prediction accuracy.

Method outperforms independent models on synthetic and real data.

Provides a new framework for analyzing multilayer network structures.

Abstract

Many recent developments in network analysis have focused on multilayer networks, which one can use to encode time-dependent interactions, multiple types of interactions, and other complications that arise in complex systems. Like their monolayer counterparts, multilayer networks in applications often have mesoscale features, such as community structure. A prominent type of method for inferring such structures is the employment of multilayer stochastic block models (SBMs). A common (but {potentially} inadequate) assumption of these models is the sampling of edges in different layers independently, conditioned on the community labels of the nodes. In this paper, we relax this assumption of independence by incorporating edge correlations into an SBM-like model. We derive maximum-likelihood estimates of the key parameters of our model, and we propose a measure of layer correlation that reflects the similarity between connectivity patterns in different layers. Finally, we explain how to use correlated models for edge "prediction" (i.e., inference) in multilayer networks. By taking into account edge correlations, prediction accuracy improves both in synthetic networks and in a temporal network of shoppers who are connected to previously-purchased grocery products.

Tables

Table 1. Table 1: Edge-prediction probabilities for various correlated multilayer and monolayer network models.

Model	$\mathbb{P}\left(A_{ij}^2=1|A_{ij}^1=1\right)$	$\mathbb{P}\left(A_{ij}^2=1|A_{ij}^1=0\right)$
Corr. ER	$q/p_1$	$\left(p_2-q\right)/\left(1-p_1\right)$
Corr. SBM	$q_{rs}/p_{rs}^1$	$\left(p_{rs}^2-q_{rs}\right)/\left(1-p_{rs}^1\right)$
Corr. CM	$\sqrt{{\theta }_i^1{\theta }_j^1{\theta }_i^2{\theta }_j^2}q/\left({\theta }_i^1{\theta }_j^1{p}_1\right)$	$\left({\theta }_i^2{\theta }_j^2{p}_2-\sqrt{{\theta }_i^1{\theta }_j^1{\theta }_i^2{\theta }_j^2}q\right)/\left(1-{\theta }_i^1{\theta }_j^1{p}_1\right)$
Corr. DCSBM	$\sqrt{{\theta }_i^1{\theta }_j^1{\theta }_i^2{\theta }_j^2}{q}_{rs}/\left({\theta }_i^1{\theta }_j^1{p}_{rs}^1\right)$	$\left({\theta }_i^2{\theta }_j^2{p}_{rs}^2-\sqrt{{\theta }_i^1{\theta }_j^1{\theta }_i^2{\theta }_j^2}{q}_{rs}\right)/\left(1-{\theta }_i^1{\theta }_j^1{p}_{rs}^1\right)$
SBM	$p_{rs}^2$
DCSBM	${\theta }_i^2{\theta }_j^2{p}_{rs}^2$

Table 2. Table 2: Pairwise layer correlations in several multiplex networks.

Domain	Network	Number	Mean	Largest correlation (corresponding layers)
		of layers	correlation
Social	CS Aarhus [71]	5	$0.27$	$0.45$ (“work” and “lunch” layers)
	Lazega law firm [72]	3	$0.39$	$0.48$ (“advice” and “co-work” layers)
	YouTube [73]	$5$	$0.12$	$0.20$ (“shared subscriptions” and “shared subscribers”)
Biological	C. elegans connectome [74]	3	$0.47$	$0.85$ (“MonoSyn” and “PolySyn” layers)
	P. falciparum genes [75]	9	$0.08$	$0.25$ (“HVR7” and “HVR9” layers)
	Homo sapiens proteins [76]	7	$0.04$	$0.29$ (“direct interaction” and “physical association”)
Other	FAO international trade [34]	364	$0.13$	$0.74$ (“Pastry” and “Sugar confectionery”)
	EU air transportation [22]	37	$0.03$	$0.39$ (“Scandinavian Airlines” and “Norwegian Air Shuttle”)
	ArXiv collaborations [77]	13	$0.07$	$0.73$ (“physics.data-an” and “cs.SI”)

Table 3. Table 3: Predictive performance of different models on the shopping data set, as measured by the AUC.

Model	AUC	AUC
	(ShoppingMod)	(ShoppingSBM)
Monolayer SBM	0.549	0.633
Correlated ER	0.724	0.743
Correlated SBM	0.742	0.793
Monolayer DCSBM	0.725	0.797
Correlated CM	0.817	0.870
Correlated DCSBM	0.818	0.875

Equations

$$\mathbb{P}(A_{ij}^1=1)$$

$$\mathbb{P}(A_{ij}^2=1)$$

$$q:=\mathbb{P}(A_{ij}^1=1,A_{ij}^2=1)$$

$$\mathbb{P}(A_{ij}^1=1,A_{ij}^2=1)$$

$$\mathbb{P}(A_{ij}^1=1,A_{ij}^2=0)$$

$$\mathbb{P}(A_{ij}^1=0,A_{ij}^2=1)$$

$$\mathbb{P}(A_{ij}^1=0,A_{ij}^2=0)$$

$$\mathbb{P}(A_{ij}^2=1|A_{ij}^1=1)$$

$$\mathbb{P}(A_{ij}^2=1|A_{ij}^1=0)$$

$$\mathbb{P}(A_{ij}^2=0|A_{ij}^1=1)$$

$$\mathbb{P}(A_{ij}^2=0|A_{ij}^1=0)$$

$$=\frac{1-p_1-p_2+q}{1-p_1}.$$

$$\rho$$

$$=\frac{q-p_1{p}_2}{\sqrt{p_1(1-p_1)p_2(1-p_2)}},$$

$$\mathbb{P}({\mathbf{A}}^1,{\mathbf{A}}^2|p_1,p_2,q)=\prod _{(i,j)\in \mathcal{E}}{q}^{A_{ij}^1{A}_{ij}^2}(p_1-q{)}^{A_{ij}^1(1-A_{ij}^2)}(p_2-q{)}^{(1-A_{ij}^1)A_{ij}^2}(1-p_1-p_2+q{)}^{(1-A_{ij}^1)(1-A_{ij}^2)}.$$

$$e_{11}$$

$$e_{10}$$

$$e_{01}$$

$$e_{00}$$

$$\mathcal{L}$$

$$+e_{00}\log (1-p_1-p_2+q).$$

$${\hat{p}}_1$$

$${\hat{p}}_2$$

$$\hat{q}$$

$$\hat{\rho }=\frac{e_{00}{e}_{11}-e_{10}{e}_{01}}{\sqrt{(e_{11}+e_{10})(e_{11}+e_{01})(e_{10}+e_{00})(e_{01}+e_{00})}}.$$

$$\mathbb{P}(A_{ij}^1=1)=p_{g_i{g}_j}^1,$$

$$\mathbb{P}(A_{ij}^2=1)=p_{g_i{g}_j}^2,$$

$$\mathbb{P}(A_{ij}^1=1,A_{ij}^2=1)=q_{g_i{g}_j}.$$

$$\mathbb{P}({\mathbf{A}}^1,{\mathbf{A}}^2|\mathbf{g},{\mathbf{p}}^1,{\mathbf{p}}^2,\mathbf{q})=\prod _{(i,j)\in \mathcal{E}}$$

$$\displaystyle\times(1-p_{g_{i}g_{j}}^{1}-p_{g_{i}g_{j}}^{2}+q_{g_{i}g_{j}})^{(1-A_{ij}^{2})(1-A_{ij}^{2})}\Bigr{]}\,.$$

$$e_{rs}^{ab}:=\left|\{(i,j)\in \mathcal{E}:A_{ij}^1=a,A_{ij}^2=b,g_i=r,g_j=s\}\right|$$

$$\mathcal{L}$$

$$+e_{rs}^{10}\log (p_{rs}^1-q_{rs})+e_{rs}^{01}\log (p_{rs}^2-q_{rs})$$

$$\displaystyle\qquad\qquad\quad+e_{rs}^{00}\log(1-p^{1}_{rs}-p^{2}_{rs}+q_{rs})\Bigr{]}\,.$$

$${\hat{p}}_{rs}^1$$

$${\hat{p}}_{rs}^2$$

$${\hat{q}}_{rs}$$

$${\hat{\rho }}_{rs}=\frac{e_{rs}^{00}{e}_{rs}^{11}-e_{rs}^{10}{e}_{rs}^{01}}{\sqrt{(e_{rs}^{11}+e_{rs}^{10})(e_{rs}^{11}+e_{rs}^{01})(e_{rs}^{10}+e_{rs}^{00})(e_{rs}^{01}+e_{rs}^{00})}}.$$

$$\mathrm{corr}(A_{IJ}^1,A_{IJ}^2)=\frac{\mathbb{E}[A_{IJ}^1{A}_{IJ}^2]-E[A_{IJ}^1]\mathbb{E}[A_{IJ}^2]}{\sigma [A_{IJ}^1]\sigma [A_{IJ}^2]},$$

$$\mathbb{E}[A_{IJ}^l]=\mathbb{P}(A_{IJ}^l=1)$$

$$=\sum _{(r,s)\in \mathcal{B}}\mathbb{P}(A_{IJ}^l=1|g_I=r,g_J=s)\mathbb{P}(g_I=r,g_J=s)$$

$$=\sum _{(r,s)\in \mathcal{B}}{\hat{p}}_{rs}^1\frac{e_{rs}}{|\mathcal{E}|}=\sum _{(r,s)\in \mathcal{B}}\frac{m_{rs}^l}{e_{rs}}\frac{e_{rs}}{|\mathcal{E}|}=\frac{m_l}{|\mathcal{E}|},$$

Full text

I Introduction

A network is an abstract representation of a system in which entities called “nodes” interact with each other via connections called “edges” [1]. Most traditionally, in a type of network called a “graph”, each edge encodes an interaction between a pair of nodes. Networks arise in many domains and are useful for numerous practical problems, such as detecting bot accounts on Twitter [2], finding vulnerabilities in electrical grids [3], and identifying potentially harmful interactions between drugs [4]. A common feature of many networks is mesoscale (i.e., intermediate-scale) structures. Detecting such structures amounts to a type of coarse-graining, providing representations of a network that are more compact than listing all of the nodes and edges. Types of mesoscale structures include community structure [5], core–periphery structure [6, 7], role similarity [8], and others. An increasingly popular approach for modeling and detecting such structures is by using stochastic block models (SBMs) [9], a generative model that can produce networks with community structure or other mesoscale structures.

For many applications of network analysis, it is important to move beyond ordinary graphs (i.e., “monolayer networks”) to examine more complicated network structures, such as collections of interrelated networks. One can study such structures through the flexible lens of multilayer networks [10, 11, 12, 13]. Similar to monolayer networks, a multilayer network consists of a collection of “state nodes” that are connected pairwise by edges. A state node is a manifestation of a “physical node” (which we will also sometimes call simply a “node”), which represents some entity, in a specific layer. Different layers may correspond to interactions in different time periods (yielding a temporal network), different types of relations (yielding a multiplex network), or other possibilities. As in the setting of monolayer networks, modeling and inferring mesoscale structures in multilayer networks is a prominent research area.

A key assumption of almost all existing models of multilayer networks with mesoscale structure is that edges are generated independently, conditioned on a multilayer partition [14, 15, 16, 17, 18, 19, 20, 21]. This independence condition applies both within each layer (which is inconsistent with the fact that real networks often include $3$-cliques and other small-scale structures) and across layers (which is inconsistent with the fact that the same nodes are often adjacent to each other in multiple layers). In this paper, we focus on relaxing the edge-independence assumption that applies to edges between the same two physical nodes in different layers. We still consider each pair of nodes independently.

In Fig. 1, we show an example of a two-layer network with both strong positive and strong negative edge correlations. Incorporating such correlations into a network model is beneficial for many applications. For example, a multiplex network of air routes, where each layer corresponds to one airline, is likely to include some popular routes that appear in multiple layers (and unpopular routes may appear only in one layer)[22]. In a temporal social network, we expect people to have repeated interactions with other people [23]; this is a stronger statement than just saying that they tend to interact more within the same community over time. Such edge persistence is also common in many bipartite user–item networks: shoppers tend to buy the same grocery products over time [24], customers of a music-streaming platform listen repeatedly to their favorite songs [25], and Wikipedia users edit specific pages several times [26].

Multilayer network models that incorporate edge correlations have many important applications. One is the inference task of edge prediction (also called link prediction), where one seeks to assign probabilities of occurrence to unobserved edges. SBMs have often been used for edge prediction for both monolayer networks [27, 4, 9] and multilayer networks [15, 20, 28]. The models that we propose should yield better performance than past efforts, as they take advantage of interlayer edge correlations in data. We discuss this in more detail in Sec. III. Another application is to graph matching [29], where one seeks to infer a latent correspondence between nodes in two different networks when one does not know the identities of the nodes. For example, one may wish to match common users between anonymized Twitter and Facebook networks. In a series of papers [30, 31, 32], Lyzinski and collaborators established conditions under which graph matching is successful. They tested their methods on correlated Erdős–Rényi (ER) networks and correlated SBMs, which we also investigate in this paper (albeit for a different purpose). We take these models further by incorporating degree correction [33], which generates networks with heterogeneous degree distributions and is important for inference using SBM. This extension may allow one to study the graph-matching problem on more realistic network models. A third application of our work is efficient computation of correlations between pairs of layers of a multilayer network. One can use such correlation estimates to quantify the similarity between different layers and potentially to compress multilayer network data by discarding (or merging) layers that are strongly positively correlated with an existing layer [34]. Previous papers have focused primarily on node-centric notions of layer similarity [35, 36, 37], whereas our correlated models yield edge-centric measures of similarity. Benefits of our approach over related studies [38, 34, 39] include the fact that correlation values cover an intuitive range (between $-1$ and $1$) and that they work equally well for quantifying layer similarity and dissimilarity.

Our edge-correlated network models are also useful for community detection. Given a multilayer network, one can design an inference algorithm that determines both the parameters that describe the edge probabilities (and correlations) and a multilayer community structure that underlies these probabilities. Solving this inference problem enables the detection of “correlated communities”. Because of the additional complexity in the model, this is bound to be more difficult than standard multilayer community detection, so we leave this inference problem for future work. Instead, for the rest of the paper, we assume that we know the block structure of a network; we infer the remaining parameters, including the correlations that are the core element of our model 111In fact, previous algorithms for SBM inference [27, 82, 15, 16] have suggested that one can formulate a full inference algorithm using an iterative approach that alternates between estimating a multilayer network’s mesoscale structure and estimating its edge probabilities and correlations. Therefore, the latter task — which is the focus of the present paper — constitutes one half of an algorithm for correlated-community detection.. We also assume that the block assignments $\bm{g}$ are the same for all layers; the case in which communities can vary arbitrarily across layers is significantly more difficult [41], and we leave its consideration for future work. With these restrictions, one can determine $\bm{g}$ using any method of choice.

Some existing models of multilayer networks incorporate interlayer dependencies by prescribing joint degree distributions [42, 43, 36], by incorporating edge overlaps [44, 45], or by modeling the appearance of new edges through preferential-attachment mechanisms [46]. The models that we describe in this paper are similar to those that were introduced in [30, 31, 32] for graph-matching purposes. Another noteworthy paper is one by Barucca et al. [47] that described a generalized version of the temporal SBM of Ghasemian et al. [17]. This generalization includes an “edge-persistence” parameter $\xi$, which gives the probability that an edge from one layer also occurs in the next temporal layer. For several reasons, we take a different approach. First, the model of [47] is specific to temporal networks, whereas we are also interested in other types of multilayer networks. Second, their model does not easily incorporate degree correction. Third, we want to include correlations explicitly in the model, rather than implicitly using the edge-persistence parameter $\xi$.

Our paper proceeds as follows. In Sec. II, we describe our models of multilayer networks with edge correlations. We start with a simple example of correlated Erdős–Rényi (ER) graphs in Sec. II.1 to make our exposition for more complicated models easier to follow. In Sec. II.2, we integrate mesoscale structures by incorporating correlations in an SBM-like model. We then introduce degree correction in Sec. II.3. For all of these models, we derive maximum-likelihood (ML) estimates both of the marginal edge-existence probabilities in each layer and of the interlayer correlations. ML estimation is common for SBMs and DCSBMs in both monolayer networks [33] and multilayer networks [16]. Although ML estimation is less powerful than performing Bayesian inference [9], the former is consistent for both SBMs and DCSBMs [48] and it recovers many common techniques for detecting mesoscale structures in networks [49]. In Sec. III, we describe how to use our models for edge prediction, and we provide some results for synthetic networks. We then proceed with applications in Sec. IV. In Sec. IV.1, we use our models to estimate pairwise layer correlations in several empirical networks. In Sec. IV.2, we use our correlated models for edge prediction in a temporal network of grocery purchases. We summarize our results in Sec. V and discuss a few ideas for future work.

II Correlated Models

In our derivations, we consider just two network layers at a time. Although this may seem limiting, we can apply our framework to generate correlated networks with more than two layers in a sequential manner (see the discussion in Sec. II.1), and we can determine pairwise layer correlations for a network with arbitrarily many layers (see the applications in Sec. IV). In Sec. V, we briefly discuss the challenges that arise when modeling three or more layers simultaneously, rather than in a sequential pairwise fashion.

Consider a network with two layers and identical sets of nodes in each layer; this is known as a node-aligned multilayer network [10]. Let $\bm{A}^{1}$ and $\bm{A}^{2}$ denote the adjacency matrices of our two network layers. As in many generative models of networks, we assume that edges in these two layers are generated by some random process, so the entries $A_{ij}^{1}$ and $A_{ij}^{2}$ are random variables. Imposing some statistical correlation between these two sets of random variables introduces interlayer correlations in the resulting multilayer network structure.

Our goal is to propose a model of correlated networks in which each layer is, marginally, a degree-corrected stochastic block model (DCSBM) [33]. However, it is instructive to first consider the simpler cases in which each layer is marginally an Erdős–Rényi random graph (see Sec. II.1) or an SBM without degree correction (see Sec. II.2). Correlated ER models and correlated SBMs have been studied previously, most notably in work by Lyzinski et al. [30, 31, 32] on the graph-matching problem. However, our use of these models for estimating layer correlations is novel, as are the correlated DCSBMs that we propose in Sec. II.3.

In monolayer SBMs, it is common to use either Bernoulli or Poisson random variables to generate edges between nodes. The former is generally more accurate, because it does not yield multiedges; however, the latter is more common, as it often simplifies calculations considerably [48, 50]. Nevertheless, we have found that Bernoulli models are simpler when incorporating correlations [24]. They also have several other advantages, including the fact that they work for both sparse and dense networks and that they can handle the entire correlation range between $-1$ and $1$. Therefore, we consider only Bernoulli models in this paper.

II.1 Correlated Erdős–Rényi Layers

II.1.1 Forward model

Assume that the intralayer networks that correspond to $\bm{A}^{1}$ and $\bm{A}^{2}$ are ER graphs from the $\mathcal{G}(n,p)$ ensemble [1] with edge probabilities $p_{1}$ and $p_{2}$. For each pair of nodes $(i,j)$, we therefore have

[TABLE]

To couple edges that connect the same pair of nodes in different layers, let

[TABLE]

denote the joint probability for an edge to occur in both layers. Unless $q=p_{1}p_{2}$, this construction implies that the random variables $A_{ij}^{1}$ and $A_{ij}^{2}$ are not independent.

The parameters $p_{1}$, $p_{2}$, and $q$ (which lie in the interval $[0,1]$) fully specify a forward model of networks with correlated ER layers. To generate a network from this model, one considers each node pair $(i,j)$ and, independently of all other node pairs, assigns values to $A_{ij}^{1}$ and $A_{ij}^{2}$ according to the following probabilities:

[TABLE]

These expressions follow from the laws of probability and from the definitions of $p_{1}$, $p_{2}$, and $q$. For these probabilities to be well-defined, it is both necessary and sufficient that $0\leq q\leq\min(p_{1},p_{2})$ and $p_{1}+p_{2}\leq 1+q$. In Fig. 2, we illustrate the feasible region for $p_{1}$ and $p_{2}$, given a value of $q$.

It is also possible to generate a correlated ER network in a sequential manner. First, one generates the adjacency matrix $\bm{A}^{1}$ by placing edges with probability $p_{1}$. One then determines the probabilities of edges in the second layer by conditioning on the first layer:

[TABLE]

With this approach, it is possible to generate networks with arbitrarily many layers by first sampling edges in the first layer, and then sampling edges in each subsequent layer by conditioning on the previous one. This kind of process is especially well-suited to temporal networks, in which layers have a natural ordering. For multiplex networks, it is more appropriate to extend Eqns. (4) to handle more than two layers.

It is also possible to parametrize correlated ER graphs in terms of the marginal Bernoulli probabilities, $p_{1}$ and $p_{2}$, and the Pearson correlation

[TABLE]

where $\mathbb{E}[\cdot]$ and $\sigma[\cdot]$, respectively, denote the mean and standard deviation of a random variable. One benefit of using $\rho$, rather than $q$, as the third model parameter is that its value is easier to interpret. A value of $\rho$ that is close to [math] indicates a weak correlation between layers, whereas values that are close to the extremes of $+1$ and $-1$ indicate a strong positive correlation and a strong negative correlation, respectively.

We can gain further intuition by considering the cases $\rho=0$, $\rho=1$, and $\rho=-1$. First, $\rho=0$ if and only if $\mathbb{P}(A_{ij}^{1}=1,A_{ij}^{2}=1)=\mathbb{P}(A_{ij}^{1}=1)\mathbb{P}(A_{ij}^{2}=1)$. That is, the correlation is [math] if and only if edges are generated independently in the two layers with marginal probabilities of $p_{1}$ and $p_{2}$. For $\rho=1$, one can show (see [24]) that $p_{1}=p_{2}=q$, which corresponds to the two layers having identical network structure. Lastly, for $\rho=-1$, we have $q=0$ and $p_{1}=1-p_{2}$ (see [24]), and two nodes are adjacent in one layer if and only if they are not adjacent in the other layer.

II.1.2 Maximum-likelihood parameter estimates

We now derive ML estimates of the parameters $p_{1}$, $p_{2}$, and $q$. Let $\mathcal{E}$ denote the set of node pairs that can form edges. For undirected networks without self-edges, there are $|\mathcal{E}|=N(N-1)/2$ such node pairs to consider, where $N$ is the number of physical nodes. By contrast, $|\mathcal{E}|=N(N-1)$ when generating directed networks without self-edges. With this general notation, all our derivations in Sec. II are valid for both directed and undirected networks, with or without self-edges. (They are also valid for bipartite networks [24].) We consider each pair of nodes $(i,j)\in\mathcal{E}$ independently when generating edges, so the likelihood of observing adjacency matrices $\bm{A}^{1}$ and $\bm{A}^{2}$ is

[TABLE]

It is helpful to introduce the following notation:

[TABLE]

These quantities correspond, respectively, to the number of node pairs that are adjacent in both layers, are adjacent in the first layer but not in the second, are adjacent in the second layer but not in the first, and are not adjacent in either layer. Using this notation and taking the logarithm of (7), we arrive at the following expression for the log-likelihood:

[TABLE]

When fitting our model to network data, the quantities $e_{11}$, $e_{10}$, $e_{01}$, $e_{00}$ are all known; and we seek to determine the values of $p_{1}$, $p_{2}$, and $q$ that are best explained by the data. To do so, we maximize the log-likelihood (II.1.2) by setting its partial derivatives to [math] 222The Hessian is negative definite at the critical point that one obtains by setting these derivatives to [math], so this point is a local maximum (as opposed to a local minimum or saddle point) of the log-likelihood function.. We obtain

[TABLE]

In all three expressions, the denominator is equal to the number of potential edges (i.e., the cardinality of $\mathcal{E}$). Additionally, let $m_{1}=e_{11}+e_{10}$ and $m_{2}=e_{11}+e_{01}$ denote the number of observed edges in the first and the second layers, respectively. It follows that the ML estimate $\widehat{p}_{1}$ is equal to the number of observed edges in layer 1 divided by the number of potential edges, and an analogous relation holds for $\widehat{p}_{2}$. The estimate $\widehat{q}$ is equal to the number of node pairs that are adjacent in both layers divided by the total number of node pairs. These results match our intuition.

We obtain an estimate of the Pearson correlation $\rho$ between the two layers by substituting the ML estimates $\widehat{p}_{1}$, $\widehat{p}_{2}$, and $\widehat{q}$ into Eqn. (II.1.1) to obtain

[TABLE]

One can show that maximizing the log-likelihood (II.1.2) with respect to $p_{1}$, $p_{2}$, and $\rho$ (rather than with respect to $p_{1}$, $p_{2}$, and $q$) gives the same expression for $\widehat{\rho}$, confirming that this is indeed an ML estimate of the correlation. Note that $\widehat{\rho}$ is not defined when either layer is an empty or a complete graph, as the corresponding Bernoulli random variable has a standard deviation of [math].

In App. A, we calculate the variances that are associated with the ML estimates $\widehat{p}_{1}$, $\widehat{p}_{2}$, $\widehat{q}$. We then show using a synthetic example that these scale as $1/N^{2}$, where we recall that $N$ is the number of physical nodes. These results quantify the uncertainty around the ML estimates for correlated ER models as a function of network size.

II.2 Correlated SBMs

One of the ways in which real-world networks differ from ER random graphs is that the former have mesoscale structures, such as communities [5]. We use SBMs to incorporate such structures into our correlated models.

Let $\bm{g}$ be a vector of block assignments, which we take to be identical for both network layers, and let $K$ denote the number of blocks. As we explained in Sec. I, we assume throughout the present paper that we are given $\bm{g}$, and we aim to estimate the remaining model parameters. Following terminology from [27], let $\mathcal{B}=\{1,\ldots,K\}\times\{1,\ldots,K\}$ be the set of “edge bundles” $(r,s)$, each of which is described by its own set of parameters $p_{rs}^{1}$, $p_{rs}^{2}$, and $q_{rs}$. The $K\times K$ matrices $\bm{p}^{1}$, $\bm{p}^{2}$, and $\bm{q}$ play an analogous role to $p_{1}$, $p_{2}$, and $q$ in the ER layers.

Let $g_{i}$ denote the block assignment of node $i$. A correlated two-layer SBM is described by the following set of equalities:

[TABLE]

Lyzinski et al. proposed this forward model in [31] to study the graph-matching problem. By contrast, we focus on the inverse problem of estimating the parameters $\bm{p}^{1}$, $\bm{p}^{2}$, and $\bm{q}$, given some network data.

II.2.1 Maximum-likelihood parameter estimates

As in Sec. II.1, suppose that we consider each node pair $(i,j)$ independently. The likelihood of observing adjacency matrices $\bm{A}^{1}$ and $\bm{A}^{2}$ is then

[TABLE]

In this product, each factor depends on $i$ and $j$ only via their block memberships $g_{i}$ and $g_{j}$, so we can combine several terms. First, define

[TABLE]

for $(a,b)\in\{(1,1),(1,0),(0,1),(0,0)\}$, in analogy with $e_{11}$, $e_{10}$, $e_{01}$, and $e_{00}$ from Sec. II.1. We can then write the log-likelihood as

[TABLE]

The advantage of writing the log-likelihood as in (II.2.1) is that it clearly separates the contribution from different edge bundles. Using the results for ER layers from Sec. II.1, we immediately obtain (without further calculations) the following ML parameter estimates:

[TABLE]

where $m_{rs}^{1}$ and $m_{rs}^{2}$ denote the number of edges between blocks $r$ and $s$ in layers $1$ and $2$, respectively, and $e_{rs}$ is the number of possible edges between nodes in block $r$ and block $s$. When there is a single edge bundle (i.e., when we do not assume any block structure in a network), the ML estimates (15)–(17) recover those that we obtained for correlated ER networks in Sec. II.1. Each edge bundle also has a corresponding Pearson correlation, whose ML estimate is

[TABLE]

In applications to temporal consumer–product networks, we find that different edge bundles have vastly different correlation values [24]. We anticipate that other empirical multilayer networks have similar properties.

II.2.2 Effective correlation

Although having different correlation values for different edge bundles can be useful, it is also helpful to have a single correlation measurement for a given multilayer network. For example, one may wish to use such a network diagnostic for one of the purposes that we outlined in Sec. I. One way to define an “effective correlation” is to first sample two node indices, $I$ and $J$, uniformly at random and then compute the Pearson correlation of the random variables $A_{IJ}^{1}$ and $A_{IJ}^{2}$. That is,

[TABLE]

where we use capital letters for the node indices $I$ and $J$ to emphasize that they are random variables.

We can calculate each term on the right-hand side of (19) by conditioning on the block assignments of the randomly chosen nodes $I$ and $J$. First, for $l\in\{1,2\}$, we have

[TABLE]

where $m_{l}$ denotes the number of edges in layer $l$. The expression (II.2.2) is the same as the probability $p_{l}$ of generating an edge in layer $l$ for the ER case. Because $A_{IJ}^{l}$ is a Bernoulli random variable (in other words, it can only take the values $1$ or [math]), its standard deviation is

[TABLE]

Lastly,

[TABLE]

The estimated value of the effective correlation is thus

[TABLE]

which recovers the value in (12) for ER layers (i.e., without any block structure in the model). We stress that there is no reason a priori to expect this outcome. In fact, the analogous result does not hold for Poisson models [24]. In the present case, the fact that there is such a correspondence between models is convenient for practical reasons, as it implies that one can perform the simpler calculations from Sec. II.1 to obtain correlation estimates between network layers, even for networks with nontrivial mesoscale structure.

II.3 Correlated Degree-Corrected SBMs

The models that we have discussed thus far generate networks in which nodes in the same block have the same expected degree. SBMs that make this kind of assumption tend to perform poorly when they are used to infer mesoscale structure in real networks, many of which have highly heterogeneous degree distributions. This observation led to the development of degree-corrected SBMs (DCSBMs) [33]. We expect that such adjustments can also make a difference when modeling edge correlations, so we now extend the model from Sec. II.2 to incorporate degree correction.

We continue to work with two-layer networks, which we again specify in terms of two intralayer adjacency matrices, $\bm{A}^{1}$ and $\bm{A}^{2}$, with a common block structure that we specify with a vector $\bm{g}$. For each node pair $(i,j)\in\mathcal{E}$, we place edges in the two layers according to the probabilities

[TABLE]

We will soon justify the expression in (24). The quantities $\theta_{i}^{l}$ and $\theta_{j}^{l}$, with $l\in\{1,2\}$, are the degrees of nodes $i$ and $j$, normalized by the mean degrees. We calculate these quantities directly from an input degree sequence, so they are not model parameters. For undirected and unipartite networks, $\theta_{i}^{l}=d_{i}^{l}\,/\,\langle d^{\,l}\rangle$, where $i\in\mathcal{N}$ and $\langle d^{\,l}\rangle$ is the mean degree in layer $l$. This normalization recovers the model in Sec. II.2 when $\theta_{i}^{l}=1$ (i.e., when all nodes have the same degree). The model parameters $p_{rs}^{1}$, $p_{rs}^{2}$, and $q_{rs}$ are now edge “propensities” that, together with the degrees, control the probabilities of edges in the layers.

The probabilities in Eqns. (22)–(23) ensure that, marginally, $\bm{A}^{1}$ and $\bm{A}^{2}$ are generated according to monolayer DCSBMs [33]. It is not obvious how to model the joint probability $\mathbb{P}(A_{ij}^{1}=1,A_{ij}^{2}=1)$. In particular, it is not clear how it should depend on the observed degrees of nodes $i$ and $j$ in layers $1$ and $2$ 333One desiderata of a correlated degree-corrected model is that it satisfies the inequality $\mathbb{P}(A^{1}_{ij}=1,A^{2}_{ij}=1)\leq\mathbb{P}(A^{1}_{ij}=1)$ for all node pairs $(i,j)$. This requirement rules out some possibilities for $\mathbb{P}(A^{1}_{ij}=1,A^{2}_{ij}=1)$, for instance one where this quantity has no dependence on node degrees.. Part of the complication is that there are four such quantities for each node pair $(i,j)$. The choice from (24) works particularly well when $\rho=1$ and the normalized degree sequences $\bm{\theta}^{1}$ and $\bm{\theta}^{2}$ are the same, as it reduces to a single degree-corrected SBM that generates two identical network layers. Another sensible option is to set $\mathbb{P}(A_{ij}^{1}=1,A_{ij}^{2}=1)=\theta_{i}^{1}\theta_{j}^{1}\theta_{i}^{2}\theta_{j}^{2}q_{g_{i}g_{j}}$. This choice has the nice property that edges in a particular edge bundle $(r,s)\in\mathcal{B}$ are independent if and only if $q_{rs}=p_{rs}^{1}p_{rs}^{2}$, which matches the independence condition from Sec. II.2 for the setting without degree correction. However, this second model underperforms the one from (22)–(24) for edge prediction (see Sec. III and [24]). Consequently, for the rest of this paper, we use the model from Eqns. (22)–(24) as our correlated DCSBM.

II.3.1 Maximum-likelihood parameter estimates

When writing the log-likelihood for correlated DCSBMs, we can ignore any additive terms that only involve known quantities, such as the normalized degrees $\theta_{i}^{l}$. We can thus write

[TABLE]

As in Sec. II.2, we seek to maximize $\mathcal{L}$ with respect to the parameters $p_{rs}^{1}$, $p_{rs}^{2}$, and $q_{rs}$ by setting the corresponding derivatives to [math]. However, degree-corrected models have the crucial complication that node pairs $(i,j)$ in the same edge bundle $(r,s)$ are no longer stochastically equivalent (i.e., the corresponding entries of the adjacency matrix are no longer sampled from independent, identically distributed random variables), so their contributions to the log-likelihood are no longer the same in general. Consequently, the ML equations for correlated DCSBMs involve $\mathcal{O}(N^{2}/K^{2})$ terms, making them more difficult to solve efficiently.

In certain cases, we are able to make some approximations that make these ML equations easier to solve. Recall that $\theta_{i}^{l}=1$ if the degree of node $i$ is equal to the mean degree in layer $l$. For $(i,j)\in\mathcal{E}$, we write

[TABLE]

If the degree distribution is narrow, such that all node degrees are close to the mean degree, then $\varepsilon_{ij}^{1}$ and $\varepsilon_{ij}^{2}$ are small parameters (which can be either positive or negative). In this case, a first-order Taylor expansion yields

[TABLE]

We also calculate

[TABLE]

and

[TABLE]

Using the approximations (26)–(28), we expand the first derivatives of $\mathcal{L}$ to first order in $\varepsilon_{ij}^{1}$ and $\varepsilon_{ij}^{2}$. (See [24] for details.) This calculation yields the following system of equations:

[TABLE]

In these equations, we defined $e_{rs}^{11}$, $e_{rs}^{10}$, $e_{rs}^{01}$, and $e_{rs}^{00}$ as in Sec. II.2. Additionally, we set

[TABLE]

We can efficiently calculate all of these quantities from the matrices $\bm{A}^{1}$ and $\bm{A}^{2}$.

The system of equations (29) reduces to the analogous equations for correlated SBMs if we ignore all of the terms that depend on $\varepsilon_{ij}^{l}$ (i.e., the terms that correspond to perturbations of the degrees from their mean values). The zeroth-order solution that we obtain from ignoring these terms provides a good initialization of a numerical algorithm to solve (29) for the parameters $p_{rs}^{1}$, $p_{rs}^{2}$, and $q_{rs}$. In practice, when using correlated DCSBMs for edge prediction (see Sec. III), we find that using a first-order approximation to determine $p_{rs}^{1}$, $p_{rs}^{2}$, and $q_{rs}$ gives results that are almost identical to those from the zeroth-order approximation [see Eqns. (15)–(17)]. Consequently, we suggest using these zeroth-order approximation for edge-prediction applications, given that they are straightforward to calculate and have negligible impact on the quality of the results. For large networks, we also obtain a noticeable improvement in calculation speed when using these approximations.

In App. B, we compare the parameters that we estimate using the approximate system of equations (29) with those from the log-likelihood (II.3.1). As expected, the quality of the approximation depends on the shape of the degree distribution, with larger discrepancies between the two approaches for broader degree distributions.

II.3.2 Correlation values

For the SBMs without degree correction from Sec. II.2, node pairs $(i,j)$ from a given edge bundle $(r,s)$ have the same Pearson correlation $\rho_{rs}$. This no longer holds for degree-corrected models. Instead, each node pair $(i,j)$ has its own correlation value

[TABLE]

As in our earlier expansions of the ML equations, we approximate $\varrho_{ij}$ to first order in $\varepsilon_{ij}^{1}$ and $\varepsilon_{ij}^{2}$. We obtain

[TABLE]

Ignoring terms that depend on $\varepsilon_{ij}^{1}$ and $\varepsilon_{ij}^{2}$ (i.e., terms that correspond to perturbations of the degrees from their mean values), we obtain $\varrho_{ij}\approx\rho_{rs}$. This approximation works especially well when $p_{rs}^{1}$ and $p_{rs}^{2}$ are also small, such that their respective network layers are sparse.

The case $q_{rs}=p_{rs}^{1}p_{rs}^{2}$ requires separate consideration 444Situations when one or more of the denominators in Eqn. (31) are $\mathcal{O}(\varepsilon)$ also require separate derivations. However, we do not treat these exceptional cases in the present paper. to avoid dividing by [math]. First-order approximations in $\varepsilon_{ij}^{1}$ and $\varepsilon_{ij}^{2}$ for this case give

[TABLE]

In particular, the zeroth-order solution gives $\varrho\approx 0$, in agreement with the SBM without degree correction from Sec. II.2.

III Edge Prediction

The aim of edge prediction (also called “link prediction”) in networks is to infer likely missing edges and/or spurious edges [54]. Edge prediction is useful for filling in incomplete data sets, such as protein-interaction networks (in which edges are often established as a result of costly experiments) [55] or terrorist-association networks (which are typically constructed based on partial knowledge) [56]. In the context of bipartite user–item networks, edge-prediction techniques provide candidates for personalized recommendations.

One can perform edge prediction in either a supervised or an unsupervised fashion. We briefly discuss each of these types of approaches.

Supervised methods rely on models that learn how a specified set of features relates to the presence or absence of edges. Existing methods that take advantage of multilayer structure for edge prediction typically do so through the specification of multilayer features. These include aggregations of monolayer features [57], as well as path-based [58] and neighborhood-based [59, 60] features that consider multiple layers. Although many of these features depend indirectly on the similarity between different layers, none of these methods quantify the level of correlation or use it for edge prediction.

With unsupervised methods, one obtains a ranking of node pairs such that edges are more likely among higher-ranked pairs. Common approaches include ones that are based on probabilistic models and ones that are based on similarity indices (like the Jaccard index or the Adamic–Adar index) [54]. An example of the former for multilayer networks is a method that maps each network layer independently to a hyperbolic space and then uses the hyperbolic distance between nodes in one layer to predict edges in another layer [61]. This work used node-centric notions of correlation and thereby complements the edge-centric perspective of our work.

Methods that rely on similarity indices include those that first generate latent states (i.e., so-called “embeddings”) for the nodes and then rank node pairs according to the similarities of these embedding vectors [62, 63]. Tillman et al. [64] used layer-level correlations to combine monolayer similarity indices into a single score. There are also approaches that implicitly take advantage of similarities across layers, such as by extracting common higher-order structures (specifically, subgraphs with three or more nodes) and looking for patterns that differ by exactly one edge [65].

Variants of SBMs are popular choices for unsupervised edge-prediction methods [55, 4, 56, 9], including in multilayer settings [15, 20]. As an example, in a monolayer degree-corrected Bernoulli SBM, the probability that two nodes, $i$ and $j$, are adjacent according to the model is $\mathbb{P}(A_{ij}=1)=\theta_{i}\theta_{j}p_{g_{i}g_{j}}$. The pairs $(i,j)$ for which these probabilities are relatively large but which are not adjacent (i.e., with $A_{ij}=0$) in the actual network of interest produce a list of likely candidates for missing edges. Similarly, pairs $(i,j)$ for which these probabilities are small but which are adjacent (i.e., with $A_{ij}=1$) in the actual network may be spurious edges.

III.1 Edge Prediction Using Correlated Models

There have been several recent attempts to perform edge prediction in multilayer networks [15, 20, 28]. All of these methods use multilayer information to infer mesoscale structures in networks, but then they perform edge prediction independently in each layer, conditioned on the inferred mesoscale structure and any other model parameters. In particular, when using one of these approaches, observing that two nodes are adjacent in one layer has no bearing on their probability to be adjacent in another layer. We aim to use our correlated models to overcome this limitation.

As in our prior discussions, consider a network with two layers with intralayer adjacency matrices $\bm{A}^{1}$ and $\bm{A}^{2}$, and let $\bm{g}$ be the shared block structure of these layers. Our goal is to predict edges in the second layer, conditioned on the adjacency structure of the first layer. For each node pair $(i,j)\in\mathcal{E}$, the key quantities to calculate are the probabilities $\mathbb{P}(A_{ij}^{2}=1|A_{ij}^{1}=1)$ and $\mathbb{P}(A_{ij}^{2}=1|A_{ij}^{1}=0)$ for $i$ and $j$ to be adjacent in the second layer, conditioned on them either being adjacent or non-adjacent in the first layer. For example, using the correlated Bernoulli SBM from Sec. II.2 (which has no degree correction), we have

[TABLE]

This set of probabilities is the same across all node pairs $(i,j)$ from the same edge bundle $(r,s)$. Now suppose that we have a positive correlation in this edge bundle, so $\rho_{rs}>0$. From the definition of the Pearson correlation, it follows that $q_{rs}>p_{rs}^{1}p_{rs}^{2}$. We then find that $\mathbb{P}(A_{ij}^{2}=1|A_{ij}^{1}=1)>p_{rs}^{2}$ and $\mathbb{P}(A_{ij}^{2}=1|A_{ij}^{1}=0)<p_{rs}^{2}$, whereas using a monolayer SBM would entail that $\mathbb{P}(A_{ij}^{2}=1)=p_{rs}^{2}$. Therefore, the effect of incorporating correlations into our edge-prediction model when these correlations are positive is (1) to increase the probability that nodes $i$ and $j$ are adjacent in the second layer when the corresponding edge also exists in the first layer and (2) to decrease this probability when the corresponding edge is absent from the first layer. The effect is reversed for negative correlations.

In Table 1, we summarize the two key probabilities (33)–(34) for four different correlated models, alongside the probabilities $\mathbb{P}(A_{ij}^{2}=1)$ for monolayer SBMs and DCSBMs. We include a correlated “configuration model” (CM) [66], which is a special case of the degree-corrected SBM from Sec. II.3 when there is only one block. (Alternatively, one can think of correlated CMs as extensions of correlated ER models that incorporate degree correction.) We use all of these models for edge prediction in synthetic networks in Sec. III.2 and in consumer–product networks in Sec. IV.2. In particular, the two monolayer models are baselines that we hope to outperform using our correlated models.

III.2 Tests on Synthetic Networks

We use $K$-fold 555Note that this $K$ is different from the one that denotes the number of blocks elsewhere in this paper. cross-validation to assess the performance of the models from Table 1 on the edge-prediction task. In machine learning, this is an effective way to measure predictive performance [54]. After partitioning a given data set into $K$ parts, one fits a model to $K-1$ of these subsets and uses it to make predictions on the remaining (i.e., “holdout”) set. One uses each subset once as a holdout, so one does this process $K$ times in total. For our problem, we perform $5$-fold cross-validation (which is a standard choice in the machine-learning literature) by splitting the data in the second layer of a given network into $5$ subsets. Effectively, this consists of hiding $20\%$ of the entries of the adjacency matrix $\bm{A}^{2}$, such that we do not know whether they are edges or not. We then train a model on $100\%$ of the entries of $\bm{A}^{1}$ and $80\%$ of the entries of $\bm{A}^{2}$, and we use it to make predictions about the $20\%$ holdout data from $\bm{A}^{2}$. We do this $5$ times to cover each choice of holdout data.

A common way to assess the performance of a binary classification model (i.e., a model that assigns one of two possible values to test data) is by using a receiver operating characteristic (ROC) curve. An ROC curve plots the true-positive rate (TPR) of a classifier versus the false-positive rate (FPR) for various choices of a threshold. Many models — including those that are used for edge prediction in networks — make probabilistic predictions, so specifying a threshold is necessary to convert these into binary predictions. Lowering the threshold increases both the TPR and the FPR. A model has predictive power if the former grows faster than the latter, such that the entire ROC curve lies above the diagonal line $\mbox{TPR}=\mbox{FPR}$, which gives the performance of a random classifier. As a single summary measure of a model’s predictive performance, it is common to report the area under an ROC curve (AUC). Larger AUC values are better, with a value of $\mathrm{AUC}=0.5$ indicating equal success as random guessing and $\mathrm{AUC}=1$ corresponding to perfect prediction. Even in the latter case, one still needs to determine a choice of threshold that completely separates true positives from false positives.

The AUC is not the only possible quantity to assess edge-prediction performance, although it is very common [15, 68]. In many networks, the number of edges is much smaller than the number of non-edges (i.e., pairs of nodes that are not adjacent). In situations with such an imbalance, the area under the precision–recall (PR) curve is more sensitive than the AUC to variations in model performance. Nevertheless, we use the AUC because it has an intuitive interpretation (specifically, as the probability that the underlying model ranks a true positive above a true negative) that allows us to establish the result in App. C. Additionally, the conclusions that we draw in Sec. IV.2 about the performance of different models do not change significantly if we use PR curves rather than ROC curves.

We now describe how to generate synthetic network benchmarks that are suitable for testing the models in Table 1. We construct these networks so that they have two tunable parameters: the Pearson correlation $\rho\in[-1,1]$ and a community-mixing parameter $\mu\in[0,1]$ that controls the strength of the planted mesoscale structure. See Bazzi et al. [14] for more details about the definition of $\mu$. One can also explicitly control the degree distribution, such as by including a parameter $\eta_{k}$ for the slope of a truncated power law (e.g., as used in [14] to sample a degree sequence in each layer). For the experiments in this section, we fix $\eta_{k}=-2$ and use a minimum degree of $k_{\mathrm{min}}=10$ and a maximum degree of $k_{\mathrm{max}}=50$. It would be interesting to explore the performance gap between degree-corrected models and models without degree correction as one varies $\eta_{k}$, $k_{\mathrm{min}}$, and $k_{\mathrm{max}}$, although we do not do so in the present paper. Lastly, for our numerical experiments in this section, we use networks with $N=2000$ nodes in each layer and $n_{c}=5$ communities, with community sizes sampled from a flat Dirichlet distribution (i.e., one with $\theta=1$ in the notation of [14]).

We examine two versions, which we call CorrSBM and CorrDCSBM, of a correlated benchmark that is parametrized by the correlation $\rho$ and the community-mixing parameter $\mu$. For both versions of the benchmark, we generate the (undirected and unipartite) adjacency matrix $\bm{A}^{1}$ of the first layer in the same way. Specifically, given $\mu$ and degree-distribution parameters $\eta_{k}$, $k_{\mathrm{min}}$, and $k_{\mathrm{max}}$, we use the code from [69] to generate $\bm{A}^{1}$ and its associated block structure $\bm{g}$. We fit a monolayer model — either an SBM or a DCSBM, depending on the selected version of the benchmark — to $\bm{A}^{1}$ to obtain the marginal edge propensities $\bm{p}^{1}$ for the first layer. We then choose $\bm{p}^{2}$ in one of two ways. For $\rho\in[0,1]$, we set $\bm{p}^{2}=\bm{p}^{1}$, which ensures that we can generate networks with correlations that cover the entire range from [math] to $1$. For $\rho\in[-1,0]$, we set $\bm{p}^{2}=\bm{\mathbbm{1}}-\bm{p}^{1}$, where $\bm{\mathbbm{1}}$ is a matrix with all entries equal to $1$; this ensures that we can generate networks with correlations that cover the entire range from $-1$ to [math]. Given $\bm{p}^{1}$, $\bm{p}^{2}$, and $\rho$, we then determine $\bm{q}$ using either the correlated SBM of Sec. II.2 or the correlated DCSBM of Sec. II.3. For the CorrDCSBM benchmark with $\rho\geq 0$, we set the normalized degrees $\theta_{i}^{2}$ to be equal to the corresponding quantities $\theta_{i}^{1}$ from the first layer. Again, this choice ensures that we can generate networks all the way to $\rho=1$. (We also implemented a version of this benchmark that samples degrees independently in the second layer, and we found qualitatively similar results.) The final step consists of generating $\bm{A}^{2}$ given $\bm{A}^{1}$, the propensities $\bm{p}^{2}$ and $\bm{q}$, and (for the CorrDCSBM benchmark only) the normalized degree sequences $\bm{\theta}^{1}$ and $\bm{\theta}^{2}$ for both layers. To perform this step, we first compute edge probabilities using either of the correlated models from Secs. II.2 and II.3, and we then generate edges independently according to these probabilities.

We now present our results for the two variants of the benchmark. In Fig. 3, we show sample ROC curves for one network that we create using the CorrDCSBM benchmark with $\mu=0.3$ and $\rho=0.5$. We compare the performance of our correlated models with a monolayer DCSBM baseline, which performs edge prediction using only information from the second network layer. Two of the correlated models outperform this baseline, and the correlated CM performs comparably well (i.e., it has a similar AUC).

In Fig. 4, we show results for the CorrSBM benchmark for two choices of the community-mixing parameter $\mu$ and several values (both positive and negative) of the Pearson correlation $\rho$. As expected, the AUC values for monolayer SBMs are independent of $\rho$, whereas the predictive performance of correlated ER models and correlated SBMs improves as we increase $|\rho|$. In particular, when $|\rho|=1$, the two correlated models make perfect predictions. When $\rho=0$, the performance of the correlated ER model is indistinguishable from chance (because $\mathrm{AUC}\approx 0.5$), whereas correlated SBMs have identical performance to monolayer SBMs. The gap between the two correlated models is smaller for $\mu=0.8$ than for $\mu=0.3$, because the underlying block structure is weaker in the former case than in the latter. The AUC of the monolayer baseline is also smaller for $\mu=0.8$ than for $\mu=0.3$.

One striking feature in Fig. 4 is that all curves are approximately straight lines (to within sampling error). It makes sense that the performance of monolayer SBMs does not vary with $\rho$, as these models do not use any information from the other layer, but the linear dependence on $\rho$ of the other two curves is less intuitive. For the correlated ER model, we can establish rigorously (see App. C) that the AUC is approximately equal to $(1+|\rho|)/2$ when $p_{1}\approx p_{2}$ or when $p_{1}\approx 1-p_{2}$. Given that the correlated SBM curves from Fig. 4 also exhibit a linear dependence on $\rho$, we believe that it is possible to establish similar results for correlated models that incorporate mesoscale structure. These results have practical importance, as they allow one to quickly estimate the additional benefits of using correlated models instead of monolayer SBMs for edge prediction.

In Fig. 5, we show results for the CorrDCSBM benchmark for two choices of the community-mixing parameter $\mu$ and nonnegative 666The case $\rho\leq 0$ requires a different way of sampling normalized degrees $\theta_{i}^{2}$ of nodes in the second layer to be able to take $\rho$ all the way to $\rho=-1$. More precisely, when $\rho$ is close to $-1$, the second layer is dense if the first layer is sparse (and vice versa), so setting $\bm{\theta}^{1}=\bm{\theta}^{2}$ is inadequate. values of the Pearson correlation $\rho$. As expected, when $\rho=0$, correlated DCSBMs perform similarly to monolayer DCSBMs. As in Fig. 4, the performance of the monolayer model is roughly independent of $\rho$, whereas the two correlated models perform better as $\rho$ increases. The gap between the two correlated models narrows substantially as one increases $\mu$ from $0.3$ to $0.8$.

IV Applications

We discuss two applications of correlated multilayer-network models to the analysis of empirical networks. In Sec. IV.1, we report pairwise layer correlations for several multiplex networks of different sizes. In Sec. IV.2, we consider a temporal bipartite network of customers and products. Using an approach similar to that from Sec. III, we demonstrate that correlated multilayer models have a better edge-prediction performance than monolayer baselines.

IV.1 Layer Correlations in Empirical Networks

We now calculate pairwise layer correlations using the formula (21). Recall that this expression gives the effective correlation between two layers, assuming that they have identical block structures (although their edge-propensity parameters can be different). Crucially, this calculation does not require that one first determines the underlying block structure. In fact, as we demonstrated in Sec. II.2, the effective correlation of a correlated SBM recovers the correlation of a correlated ER graph, which is straightforward to compute. Accounting for node degrees, as we did in Sec. II.3 for correlated DCSBMs, significantly increases the complexity of such a calculation. Additionally, as we showed in Sec. II.3, correlations using a degree-corrected model are rather similar to those that one obtains without degree correction.

In Table 2, we report the mean pairwise layer correlation for $9$ multiplex networks. (See App. D for descriptions of these data sets.) To provide additional insight into these networks, we also report the two layers with the largest effective correlations.

We make a few observations about some of the results in Table 2. For the C. elegans connectome, the layers that correspond to two types of chemical synapses are highly correlated with each other, and their correlation to the layer of electrical synapses is comparatively lower. For the European Union (EU) air transportation network, the two most correlated layers are those that correspond to Scandinavian Airlines and Norwegian Air Shuttle flights; this is consistent with the findings in [38] (which were based on a different method for quantifying layer similarity). For the network of arXiv collaborations between network scientists, the two most similar categories are “physics.data-an” (which stands for “Data Analysis, Statistics and Probability”) and “cs.SI” (which stands for “Social and Information Networks”). One hypothesis is that these two labels are often used together in papers; such common usage results in a large edge overlap between the corresponding layers and hence in a large correlation value.

To quantify edge correlations at a more granular level, one can first infer block assignments $\bm{g}$ and then calculate correlations $\rho_{rs}$ between all block pairs $(r,s)$. One possible finding from such a calculation may be that correlations between Scandinavian Airlines and Norwegian Air Shuttle routes are significantly larger in certain geographical regions than in others.

IV.2 Edge Prediction in Shopping Networks

The data-science company dunnhumby gave us access to “pseudonymized” transaction data from stores of a major grocery retailer in the United Kingdom. The data were pseudonymized by replacing personally identifiable information with numerical IDs, rendering it impossible to identify individual shoppers. For our analysis, we aggregate transactions over fixed time windows to construct bipartite networks of customers and products. We refer to these structures as “shopping networks”. Because some purchases occur in higher volumes than others, it is useful to incorporate edge weights. Given a customer $i$ and a product $j$, the item-penetration weight is equal to the fraction of all of the items purchased by customer $i$ that are product $j$. The basket-penetration weight is equal to the fraction of all baskets (i.e., distinct shopping trips) of customer $i$ that include product $j$. See the doctoral dissertation [24] for more details about these weighting schemes.

We now apply the edge-prediction methodology from Sec. III to temporal shopping networks, in which edges and edge weights can change from changes in shopping behavior, with a fixed set of customers and a fixed set of products. We construct networks with two layers, which cover the three-month time periods of March–May 2013 and June–August 2013, respectively. Using the same underlying transaction data, we construct two networks for which we determine the vector $\bm{g}$ of block assignments (the same one for both layers 777The assumption that communities are identical across layers is a reasonable one for this data set. When performing community detection without this restriction, we find that more than 90% of nodes stay in the same community across temporal layers. See [24] for details.) in different ways. For the first network (which we call ShoppingMod), we use basket-penetration weights for the edges and apply multilayer modularity maximization [79, 80] to the weighted network to determine community assignments $\bm{g}$ 888We use GenLouvain [80] to perform multilayer modularity maximization. This algorithm requires the specification of (at least) two parameters [41, 79], which we set to $\gamma=1.2$ and $\omega=5.0$. The large value of $\omega$ ensures that community assignments are identical across layers.. For the second network (which we call ShoppingSBM), we initially calculate item-penetration weights, and we then apply a threshold to remove edges whose weight is below the median weight (i.e., approximately $50\%$ of the edges). We fit a degree-corrected SBM to the resulting unweighted network using the belief-propagation algorithm from [82]. We expect better edge-prediction performance for the second network, because we detect its block structure using an SBM (as opposed to using modularity maximization, which is more restrictive).

As with our tests on synthetic networks in Sec. III.2, we use $5$-fold cross-validation to assess edge-prediction performance. We summarize the AUC values of our various correlated multilayer models and the monolayer baselines in Table 3, and we show sample ROC curves in Fig. 6. We make a few observations about these results. First, our approximation $\mathrm{AUC}\approx(1+\rho)/2$ for correlated ER models is very accurate for these two networks, whose correlations are approximately $0.44$ and $0.48$, respectively. Second, our correlated multilayer models outperform the monolayer baselines for both networks. In particular, the very simple correlated ER model — which assigns one of two probabilities to edges, as indicated in Table 1 — performs about as well as the more sophisticated monolayer DCSBM for the ShoppingMod network. Third, as expected, AUC values are systematically larger for ShoppingSBM than for ShoppingMod. Finally, although incorporating mesoscale structure leads to better performance when there is no degree correction, this does not seem to be the case for degree-corrected models, as correlated DCSBMs do not perform significantly better than correlated CMs. This is also apparent in Figs. 6(b,d), where we observe almost identical ROC curves for the two models. This result suggests that, for some networks, taking into account layer correlations and degree heterogeneity alleviates the need to also consider mesoscale structure when performing edge prediction. This observation has practical implications, as a correlated CM is much easier than a correlated DCSBM to fit to data and to use for edge prediction. However, for recommendation systems, there are situations in which fitting a correlated DCSBM is beneficial, even if its edge-prediction performance is similar to that of a correlated CM. For instance, one may wish to identify relevant customers for a chosen product, irrespective of how much they buy (i.e., their degree). SBMs are able to distinguish between customers with equal degrees and identify those with the greatest predisposition to buy a particular product, whereas CMs are not.

We also note a result from [83] that a curve lies completely above (i.e., “dominates”) another ROC curve if and only if the same relationship holds for the associated PR curves. This result implies for almost all of the curves in Fig. 6 that the rankings of our models based on AUCs are almost identical to those that we would obtain if we instead base them on the areas under PR curves. The only curves whose ranking when we use PR curves is unclear from this result are the correlated ER and monolayer SBM curves in Fig. 6(c).

It is informative to consider the source of false positives in Fig. 6. Because the correlation between layers is positive, we are likely to predict an edge where none exists when all of the following conditions hold: (1) nodes $(i,j)$ are adjacent in one layer but not the other; (2) the node pair $(i,j)$ belongs to an edge bundle $(r,s)$ with a large layer correlation $\rho_{rs}$; and (3) nodes $i$ and $j$ have large degrees. Note that condition (c) applies only for correlated models with degree correction.

V Conclusions and Discussion

We introduced models of multilayer networks in which edges that connect the same nodes in different layers are not independent. In comparison to models without edge correlations, our models offer an improved representation of many empirical networks, as interlayer correlations are a common phenomenon: flights between major airports are serviced by multiple airlines, individuals interact repeatedly with the same people, consumers often buy the same products over time, and so on. Among other potential applications, one can use our models to improve edge prediction, to study the graph-matching problem on more realistic benchmark networks, and to calculate layer correlations as insightful summary statistics for networks.

To model layer correlations, we used bivariate Bernoulli random variables to generate edges simultaneously in two network layers. (See [24] for derivations using Poisson random variables.) Correlated Bernoulli stochastic block models were proposed previously [31], although only as forward models for generating networks, rather than for performing inference given empirical data. Another key contribution of our work is a degree-corrected variant of such a model. The maximum-likelihood equations are significantly more difficult to solve in this case, but we were able to make useful simplifications with suitable approximations. Notably, these simplified equations closely approximate those for models without degree correction for networks with almost homogeneous degree distributions.

The models in the present paper that incorporate some mesoscale structure $\bm{g}$ assume that such structure is given. This setup has the benefit that one can use any desired algorithms to produce a network partition, including ones that operate on weighted or annotated networks or that use nonstandard null models in a modularity objective function. This makes our approach for analyzing correlations suitable for a wide variety of applications.

Fitting a correlated SBM to network data yields a correlation value $\rho_{rs}$ for each edge bundle $(r,s)$. We have defined an effective correlation that combines all of these values into a single measure of similarity between two layers. Notably, the value of the effective correlation is independent of a network’s mesoscale structure, making it extremely easy to compute (see Eqn. (12)). We illustrated this method of assessing layer similarity for multiplex networks from social, biological, and other domains.

Another application of our work is to edge prediction in multilayer networks. Our numerical experiments revealed that simple correlated models (e.g., a correlated configuration model or a correlated SBM without degree correction) can outperform monolayer DCSBMs in terms of AUC even for moderate correlation values. We also observed such improved performance for consumer–product networks, which have significant layer correlations ($\rho\approx 0.45$). We expect that a correlated multilayer DCSBM will typically outperform a monolayer DCSBM for most empirical networks, even when there are lower levels of correlation.

There are many interesting ways to build further on our work. For example, it would be useful to be able to model all layers simultaneously, rather than in a pairwise fashion, especially for multiplex networks (in which layers do not have a natural ordering). One challenge is that a multivariate Bernoulli distribution of dimension $L$ has $2^{L}-1$ parameters; this grows quickly with the number $L$ of layers. For a temporal setting, we have proposed generating correlated networks in a sequential way by conditioning each layer on the previous one. In some cases, it will be useful to relax this “memoryless” assumption and condition a layer on all previous layers, rather than only on the most recent layer. For example, the purchases of shoppers in December one year are strongly related not only to their purchases in November, but also to what they bought in December during the previous year.

In the present paper, we have not considered the case of nonidentical mesoscale structure across layers. (See [24] for a possible approach.) Additionally, although one can account for edge weights when fitting a block structure $\bm{g}$ to use with our models, the rest of our derivations apply only to unweighted networks. Modeling correlated weighted networks entails prescribing both edge-existence and edge-weight correlations. Yet another idea is to derive correlated models for networks with overlapping communities. Previous research [15] suggests that this can substantially improve edge-prediction performance. For correlated DCSBMs, for example, it would be useful to learn the dependence of the joint probability $\mathbb{P}(A_{ij}^{1}=1,A_{ij}^{2}=1)$ on node degrees, rather than assume the parametric form in Eqn. (24). This is a challenging problem for which maximum-likelihood estimation is unlikely to be a suitable tool, so it falls outside the scope of the present paper.

Our work is also a starting point for designing algorithms to detect correlated communities in networks by inferring $\bm{g}$ alongside other model parameters. A practical outcome of such an algorithm would be a set of communities that persist across layers if and only if the edges in those communities are sufficiently highly correlated with each other. Such an approach would offer a new interpretation of what it means for a community to span multiple layers [84, 41].

Lastly, although we performed edge prediction in an unsupervised manner, it is also possible to use estimated correlations as features in supervised models [58, 57, 59, 60] to improve their performance.

In summary, our work highlights the importance of relaxing edge-independence assumptions in statistical models of network data. Doing so provides richer insights into the structure of empirical networks, improves edge-prediction performance, and yields more realistic models on which to test community-detection, graph-matching, and other types of algorithms.

Acknowledgements.

ARP was funded by the EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling (EP/L015803/1) in partnership with dunnhumby. We are grateful to dunnhumby for providing access to grocery-shopping data and to Shreena Patel and Rosie Prior for many helpful discussions. We also thank Junbiao Lu for helpful comments.

Appendix A Variance of ML Estimates

In Sec. II, we derived ML parameter estimates for various types of correlated network models. In this appendix, we show how to obtain the variance of these estimates and we illustrate how these variances scale with network size (i.e., the number of nodes). We focus on correlated ER models and the corresponding log-likelihood (II.1.2). The same approach also works for the other types of models that we examined.

Let $\bm{\beta}=[p_{1}~{}p_{2}~{}q]^{\top}$ denote the vector of parameters for a correlated ER model. Under mild conditions [85], the ML estimate $\widehat{\bm{\beta}}$ converges in distribution to a multivariate normal distribution as the number $N$ of nodes tends to infinity. For large but finite $N$, the quantity $\widehat{\bm{\beta}}$ is distributed approximately according to $\mathcal{N}(\bm{\beta}^{*},\mathcal{I}^{-1}(\bm{\beta}^{*}))$, where $\bm{\beta}^{*}$ is the true-parameter vector and $\mathcal{I}^{-1}(\bm{\beta}^{*})$ is the inverse of the Fisher information matrix evaluated at the parameter values in $\bm{\beta}^{*}$. The diagonal entries of $\mathcal{I}^{-1}(\bm{\beta}^{*})$ provide variance estimates for $p_{1}$, $p_{2}$, and $q$.

To illustrate how these variance estimates scale with the number $N$ of nodes, we simulate correlated ER networks with $p_{1}=0.1$, $p_{2}=0.085$, and $q=0.05$ (which correspond to a correlation of $\rho=0.5$) for different network sizes. In Fig. 7, we plot the 95% confidence intervals around the ML estimates for each of the three parameters in our model. We find empirically that the variance scales with $1/N^{2}$, so the standard deviation (and thus the width of the $95\%$ confidence intervals) scales with $1/N$.

Appendix B Maximizing the Full Log-Likelihood Versus our Approximate Log-Likelihood for Correlated DCSBMs

In Sec. II.3, we derived an approximation (29) of the log-likelihood for correlated DCSBMs that enables a more efficient estimation of the parameters than maximizing the exact log-likelihood (II.3.1). We now use simulated data to illustrate how the results that one obtains using this approximation compare with maximizing the original likelihood.

We simulate networks using the CorrDCSBM benchmark (see Sec. III.2) with $N=1000$ nodes, $K=5$ communities, a mixing parameter of $\mu=0.3$, and a correlation of $\rho=0.5$. For $N=2000$ nodes, which we used in other experiments in the paper, we find that maximizing the full log-likelihood is prohibitively slow. Even for $N=1000$, obtaining estimates from the full log-likelihood takes about 30 minutes in total, whereas the calculation runs in about $5$ seconds on a typical laptop when we use the approximate log-likelihood.

The quality of our approximation depends strongly on the shape of the degree distribution: accuracy degrades for distributions with larger variances and heavier tails. We illustrate this behavior in Fig. 8, where we plot the full-likelihood estimates [see Eqn. (II.3.1)] versus the approximate-likelihood estimates [see Eqn. (29)] for two networks. For both networks, we sample degrees from truncated power-law distributions using the code in [69]. In Fig. 8(a), we choose a relatively narrow degree distribution (with $k_{\mathrm{exp}}=0$, a minimum degree of $18$, and a maximum degrees of $22$). In this case, the parameter values that we estimate using the approximate log-likelihood closely match those from the full log-likelihood. For Fig. 8(b), we choose a relatively wide degree distribution (with $k_{\mathrm{exp}}=-2$, a minimum degree of $10$, and a maximum degree of $50$). In this case, the approximate log-likelihood tends to overestimate parameter values, especially for larger values of these parameters. See the top-right corner of Fig. 8(b).

As this example illustrates, there is a trade-off between accuracy and speed. It may be possible to derive better approximations, such as by considering a second-order expansion instead of a first-order expansion with respect to the quantities $\varepsilon_{ij}^{1}$ and $\varepsilon_{ij}^{2}$ or by adding corrections for large-degree nodes to the approximate log-likelihood.

Appendix C Edge-Prediction AUC as a Function of the Pearson Correlation $\rho$

We establish the following result.

Proposition**.**

The AUC for a correlated ER model is an affine function of the Pearson correlation $\rho$. In particular, when $p_{1}\approx p_{2}$ or $p_{1}\approx 1-p_{2}$, we have $\mathrm{AUC}_{\mathrm{ER}}\approx(1+|\rho|)/2$.

Proof.

Suppose that $\rho>0$. (The case $\rho\leq 0$ is similar.) With a correlated ER model, all unobserved interactions $(i,j)$ in the second layer have one of two probabilities: $q/p_{1}$ if $A_{ij}^{1}=1$ and $(p_{2}-q)/(1-p_{1})$ if $A_{ij}^{1}=0$. Because $\rho>0$, we have $q>p_{1}p_{2}$, which implies that $q/p_{1}>(p_{2}-q)/(1-p_{1})$. Selecting a threshold between these two probabilities amounts to predicting that everything that is an edge in the first layer is also an edge in the second layer and that everything that is not an edge in the first layer is also not an edge in the second layer. Let $a$ and $b$ denote the TPR and FPR, respectively, at such an intermediate threshold. In this case, $a$ and $b$ are the coordinates of the point at which the slope of the ROC curve changes. See the illustration in Fig. 9.

By straightforward geometry,

[TABLE]

is the area under this ROC curve. The next step is to estimate $a$ and $b$. With a correlated ER model, the number of true positives is proportional 999It is not exactly equal, because with $K$-fold cross validation, we are performing a prediction on only $1/K$ of the data at a time. to $e_{11}$; the model’s prediction is correct every time that an edge that is present in the first layer is also present in the second layer. To find the TPR, one needs to divide this quantity by the number of edges (there are $e_{11}+e_{01}$ of them) in the second layer. We obtain

[TABLE]

Similarly, every time an edge that is present in the first layer is not present in the second layer counts as an incorrect prediction of the model. Therefore, the number of false positives is proportional to $e_{10}$. Dividing this by the number of non-edges in the second layer yields

[TABLE]

From (35) and (36), it follows that

[TABLE]

which is an affine function of $\rho$. When $p_{1}\approx p_{2}$ or $p_{1}\approx 1-p_{2}$, as is the case in Fig. 3, we obtain $\mathrm{AUC}_{\mathrm{ER}}\approx(1+|\rho|)/2$, as desired. Using a similar argument, one can show that the same result holds (with the same assumptions on $p_{1}$ and $p_{2}$) when $\rho\leq 0$. ∎

Given that the correlated SBM curves from Fig. 4 also appear to depend linearly on $\rho$, we believe that it is possible to establish similar results for correlated models that incorporate mesoscale structure.

Appendix D Data Sets

We provide brief descriptions of the multiplex networks that we analyzed in Sec. IV.1. For weighted networks, we disregard edge weights when calculating layer correlations. We downloaded these networks, aside from the YouTube and P. falciparum data sets, from https://comunelab.fbk.eu/data.php.

D.1 CS Aarhus

This is an undirected and unweighted social network of offline and online relationships between $N=61$ members of the Department of Computer Science at Aarhus University [71]. There are $T=5$ layers: (1) regularly eating lunch together; (2) friendships on Facebook; (3) co-authorship; (4) leisure activities; and (5) working together.

D.2 Lazega Law Firm

This directed, unweighted network encompasses interactions between $N=71$ partners and associates who work at the same law firm [72]. The network has $T=3$ layers that encode co-work, friendship, and advice relationships.

D.3 YouTube

This is an undirected, weighted network of interactions between $N=15088$ YouTube users [73]. There are $T=5$ types of interactions: (1) direct contacts (“friendships”); (2) shared contacts; (3) shared subscriptions; (4) shared subscribers; and (5) shared favorites.

D.4 C. elegans Connectome

This is a directed, unweighted network of synaptic connections between $N=279$ neurons of the nematode C. elegans [74]. There are $T=3$ layers, which correspond to electric, chemical monadic (“MonoSyn”), and chemical polyadic (“PolySyn”) junctions.

D.5 P. falciparum Genes

This is an undirected, unweighted network of $N=307$ recombinant genes from the parasite P. falciparum, which causes malaria [75]. There are $T=9$ layers that correspond to distinct highly variables regions (HVRs), in which these recombinations occur. Two genes are adjacent in a layer if they share a substring whose length is statistically significant.

D.6 Homo sapiens Proteins

This is a directed, unweighted network of interactions between $N=18222$ proteins in Homo sapiens [87]. There are $T=7$ layers, which correspond to the following types of interactions: (1) direct interactions; (2) physical associations; (3) suppressive genetic interactions; (4) association; (5) colocalization; (6) additive genetic interactions; and (7) synthetic genetic interactions. The original data is from BioGRID [76], a public database of protein interactions (for humans as well as other organisms) that is curated from different types of experiments.

D.7 Food and Agriculture Organization (FAO) Trade

This is a weighted, directed network of food imports and exports during the year $2010$ between $N=214$ countries [34]. There are $T=364$ layers, which correspond to different food products.

D.8 European Union Air Transportation

This is an undirected, unweighted network of flights between $N=450$ airports in Europe [22]. There are $T=37$ layers, each of which corresponds to a different airline.

D.9 ArXiv Collaborations

This is an undirected, weighted coauthorship network between $N=14489$ network scientists [77]. There are $T=13$ layers, which correspond to different arXiv subject areas: “physics.soc-ph”, “physics.data-an”, “physics.bio-ph”, “math-ph”, “math.OC”, “cond-mat.dis-nn”, “cond-mat.stat-mech”, “q-bio.MN”, “q-bio”, “q-bio.BM”, “nlin.AO”, “cs.SI”, and “cs.CV”.