Reconstructing Sparse Multiplex Networks with Application to Covert Networks

Abstract

Network structure provides critical information for understanding the dynamic behavior of complex systems. However, the complete structure of real-world networks is often unavailable, thus it is crucially important to develop approaches to infer a more complete structure of networks. In this paper, we integrate the configuration model for generating random networks into an Expectation–Maximization–Aggregation (EMA) framework to reconstruct the complete structure of multiplex networks. We validate the proposed EMA framework against the Expectation–Maximization (EM) framework and random model on several real-world multiplex networks, including both covert and overt ones. It is found that the EMA framework generally achieves the best predictive accuracy compared to the EM framework and the random model. As the number of layers increases, the performance improvement of EMA over EM decreases. The inferred multiplex networks can be leveraged to inform the decision-making on monitoring covert networks as well as allocating limited resources for collecting additional information to improve reconstruction accuracy. For law enforcement agencies, the inferred complete network structure can be used to develop more effective strategies for covert network interdiction.

1. Introduction

Multiplex networks, in which each layer consists of the same set of nodes but different sets of links (Figure 1), are a powerful tool for describing and analyzing the connectivity from interactions of different types among the same set of entities in complex systems [1]. Examples include social networks [2], multi-modal transportation networks [3], multi-tier supply networks [4], etc. Due to the interplay between network structure and dynamics, the knowledge of complete multiplex network structure is essential to a deep understanding of the dynamic behavior of multiplex networks and predicting future interactions between nodes in the networks [5]. Nonetheless, in many practical applications, particularly covert networks, it is extremely difficult to obtain sufficient data for constructing the complete structure of all layers of a multiplex network due to limited observing resources and privacy concerns [6].

Reconstructing the complete topology (i.e., all nodes and links) from incomplete toplogies (partially observed nodes and links) is an exceedingly challenging task primarily because (i) the number of missing nodes and links will be very large; (ii) sparse networks (e.g., illicit supply networks) usually have much lower link density than typical (overt) social networks, thereby it is often inappropriate to use many prediction methods developed for denser networks with a balanced ratio of positive and negative labels; and (iii) observations are limited to task-specific activity, excluding historic relations that are instrumental to network resilience [7,8].

On the application side, this study focuses on illicit networks, such as drug and human trafficking networks. The reason is that illicit networks compromise national security and prosperity as well as public health and safety [9]. Despite considerable advances in understanding and disrupting the operations of covert networks, a deep understanding of the structure and dynamics of covert networks remains lacking. Although the structure is crucial for understanding, monitoring, and controlling covert networks, this topological information is always incomplete [10,11]. There are multiple reasons: (i) Investigations typically focus on current operations with the aim of documenting the involvement of specific individuals in specific criminal activities. (ii) Potential criminal actors may be identified and are then often omitted from the network due to essential legal principles and civil rights protections. (iii) Criminal actors take measures to obfuscate their involvement, hiding, or disguising activities and relations because they are wary of attack by agents of the criminal justice system, as well as competitors or rivals. By modeling illicit networks as a multiplex and reconstructing the structure of multiplex covert networks, this study can be leveraged to help monitor, control, and interdict the illicit activities on covert networks.

The rest of this paper is structured as follows. Section 2 presents the literature review. Details on the EMA framework are presented in Section 3. Details on the real-world multiplex networks used to validate the proposed EMA framework are given in Section 4. The numerical experiments and the resultant results are reported in Section 5, followed by the discussion on the implication of this study on disrupting illicit networks in Section 6 and the concluding remarks in Section 7.

2. Literature Review

Network reconstruction methods can be divided into three categories: similarity-based approach, probabilistic (statistical) approach, and algorithmic (machine learning) approach.

Structural similarity-based methods assume that the nodes tend to connect to other nodes with a higher level of similarity computed by a certain distance function [12]. These approaches give similarity scores for missing or unobserved links between any pair of two nodes, then links with high similarity scores will be predicted to exist. For example, in [13], the eigenvectors of the layer adjacency matrix are used to measure the layer topological similarity and then the topological similarity (element of the layer similarity matrix) of the unconnected nodes is used to predict if a link exists between them. Berlusconi et al. [14] and Calderoni et al. [15] apply multiple similarity metrics, such as common neighbor and resource allocation, to identify missing links in a criminal network. They demonstrated that these similarity metrics can identify possible missing links in criminal networks with noise or incomplete information. However, the performance of similarity-based methods depends on the structure of specific networks [16]; therefore, prior knowledge about the networks under study is usually required to improve the predictive performance [13]. This approach is infeasible for the network reconstruction problem in this study because the nodes and links to be inferred are completely missing without any connection with observed nodes. Therefore, all scores of structural similarity between unobserved nodes and observed nodes are zeros.

Statistical approaches for network structure inference assume a network has a known structure, and then a parametric or nonparametric model is built to fit the structure, such as the stochastic block model or its variants, as well as the Dirichlet network distribution model [17,18,19,20]. Then, the model parameters are estimated using statistical methods, such as maximum likelihood estimation. The model parameters are then used to compute the formation probability of each unobserved link, such as the stochastic block model and the hierarchical structure model. In particular, Kim and Leskovec [21] developed a scalable Kronecker Expectation–Maximization (KronEM) approach to reconstruct a single-layer network with a known number of missing nodes by alternating between inferring the unobserved part of the network (E-step) and estimating the optimal model parameters (M-step) until convergence. The network structure is assumed to be described by the Kronecker graph model, and this approach does not require attributes about nodes or edges. Wu et al. [22] use a Bayesian approach to reconstruct two-layer multiplex networks from the available aggregate topology of all layers (the aggregation of all the layers into a monoplex layer) and partial observations of nodes in each layer, and the inference is conducted by leveraging the Expectation–Maximization algorithm. However, statistical approaches would be computationally challenging for large networks, especially when a large number of samples are required for inference [21].

The previous two approaches are based on computing a score for each unobserved link by defining a similarity or a probability function. However, the inference of network structure can also benefit from other algorithmic (machine learning) approaches, including supervised learning and optimization techniques. In particular, Zhang and Chen [23] develop graph neural network (GNN) models to predict missing links. However, the method is mainly developed for single-layer networks. Gao et al. [24] developed a network embedding method called interactive learning across relations that exploits existing multiple types of relational data to predict inductive links (new or unobserved relation) between existing nodes instead of inter- or intra-layer links in a known layer.

On the application side, a growing number of researchers are resorting to multilayer network models to study illicit supply chains [11,25,26,27,28] since these models allow actors to be placed in different layers representing different characteristics of covert networks. Going beyond mapping, only the co-offending activity, multilayer models help to reveal more about the social processes that connect people, which can be used to infer missing information [29,30]. However, to the best of our knowledge, a few studies have attempted to infer their multiplayer structure (topology), even though the complete structure of multiplex illicit networks is typically missing, as discussed earlier.

In this study, we develop an EMA framework built on the EM algorithm. The EMA framework is advantageous for reconstructing covert networks because:

• The sparsity of links in covert networks, i.e., the negatives (no links between nodes) significantly outnumber the positives, leads to challenges in inference and thus significantly affects the predictive performance of classical learning-based methods [31].
• It can infer the links that connect unobserved nodes of the network. This is hardly possible for methods built on structural similarities, such as Adamic–Adar index, Jaccard index, and Resource Allocation index [32], because the links to be inferred usually connect unobserved nodes that are isolated from the observed part of the network.
• It embeds observed criminal activity within the latent social structure required to sustain operations within a dynamic operational context [33,34,35,36,37]. This informs disruption efforts by improving estimates of network resiliency [38,39].

The EMA framework is an extension to the EM framework proposed in a recent work [22], but the EMA framework is more generic because it can be applied to multiplex networks with more than two layers and unknown aggregate topology, which is more common for real-world incomplete multiplex networks.

3. Methodology

3.1. Problem Description

Mathematically, a multiplex network can be described by a set $\mathit{Z}=\left\{{\mathit{Z}}^{\ell }\right|\ell =1,2,\dots ,m\}$, where m is the number of layers. ${\mathit{Z}}^{\ell }$ consists of the complete set of nodes $\mathcal{N}$ and the complete set of edges ${\mathcal{E}}^{\ell }$ of relation type ℓ. For incomplete multiplex networks, only some parts (partial topologies) of each layer are observed, i.e., a fraction c of nodes and the associated links of the complete network. If we denote the partially observed topologies as $\mathit{X}=\left\{{\mathit{X}}^{\ell }\right|\ell =1,2,\dots ,m\}$, then we have ${\mathit{X}}^{\ell }\subseteq {\mathit{Z}}^{\ell }$. In this study, we use the same notation for topology and the adjacency matrix that encodes the connectivity of the topology. We consider undirected and unweighted multiplex networks; therefore, the observed and complete adjacency matrices for each layer of a multiplex network are symmetric and binary (i.e., ${\mathit{Z}}_{ij}^{\ell },{\mathit{X}}_{ij}^{\ell }\in \{0,1\},\forall i,j\in \mathcal{N}$). Formally, the problem of this study can be given as follows.

Problem. 

Consider an undirected and unweighted multiplex network wherein each layer ℓ have the same set of nodes $\mathcal{N}$ but a different set of links ${\mathcal{E}}^{\ell }$. The multiplex network reconstruction problem is to infer the complete topology of all layers $\mathit{Z}$ given partial topologies of all layers $\mathit{X}$.

Note that the classical link prediction problem [32], in which only links are missing and typically not all links connected a node are missing, is a special case of the network reconstruction problem, in which both nodes and links can be missing, i.e., all links connected to a missing node are missing as well. As such, the network reconstruction problem is naturally more difficult than the link prediction problem.

3.2. Expectation–Maximization–Aggregation Framework

This section provides details on the proposed EMA framework for multiplex network reconstruction, including the classical EM algorithm, the network generation model, the integration of aggregate topology into the EM framework, as well as the corresponding algorithms.

3.2.1. Expectation–Maximization

The EM algorithm [40] is an iterative approach for estimation problems wherein the value of model parameters and latent variables depend on each other. For the network reconstruction problem, where the task is to find the most probable topology of the missing subnetwork, it is required to first estimate the parameters of the network model, which establish a connection between the observed subnetworks $\mathit{X}$ and the complete multiplex network $\mathit{Z}$. This framework can be formulated as follows.

We here consider a model parameterized by $\mathsf{\Theta }$. Given the observed subnetworks $\mathit{X}$ of the complete multiplex network $\mathit{Z}$, the task of classical maximum likelihood estimation is to estimate $\mathsf{\Theta }$. However, since the complete observations $\mathit{Z}$ are not available, we need to marginalize it out. Let $q\left(\mathit{Z}\right)$ be the probability density distribution of $\mathit{Z}$, i.e., ${\sum }_{\mathit{Z}}q\left(\mathit{Z}\right)=1$ and $q\left(\mathit{Z}\right)\ge 0$, then we obtain the likelihood of the observed topology given the model

$$\begin{array}{}\mathrm{(1a)}& lnp\left(\mathit{X}\mid \mathsf{\Theta }\right)& =ln\sum _{\mathit{Z}}p\left(\mathit{X},\mathit{Z}\mid \mathsf{\Theta }\right)\hfill \mathrm{(1b)}& & =ln\sum _{\mathit{Z}}q\left(\mathit{Z}\right)·\frac{p\left(\mathit{X},\mathit{Z}\mid \mathsf{\Theta }\right)}{q\left(\mathit{Z}\right)}\hfill \mathrm{(1c)}& & \ge \sum _{\mathit{Z}}q\left(\mathit{Z}\right)ln\left(\frac{p\left(\mathit{X},\mathit{Z}\mid \mathsf{\Theta }\right)}{q\left(\mathit{Z}\right)}\right).\hfill \end{array}$$

The last inequality follows from Jensen’s equality. Specifically, $ln\mathbb{E}\left[x\right]\ge \mathbb{E}[lnx]$ since $lnx$ is a concave function. The equality is attained when $\frac{p\left(\mathit{X},\mathit{Z}\mid \mathsf{\Theta }\right)}{q\left(\mathit{Z}\right)}$ is a constant, i.e., $q\left(\mathit{Z}\right)\propto p\left(\mathit{X},\mathit{Z}\mid \mathsf{\Theta }\right)$. Typically, $q\left(\mathit{Z}\right)$ is set to the posterior distribution of $\mathit{Z}$, thus we get

$$q\left(\mathit{Z}\right)=\frac{p\left(\mathit{X},\mathit{Z}\mid \mathsf{\Theta }\right)}{{\sum }_{\mathit{Z}}p\left(\mathit{X},\mathit{Z}\mid \mathsf{\Theta }\right)}=p(\mathit{Z}\mid \mathit{X},\mathsf{\Theta }).$$

(2)

The lower bound given by Equation (1c), referred to as the evidence lower bound (ELBO), is used because it is difficult to directly optimize Equation (1a) due to the missing information [41]. The EM algorithm then proceeds by alternating between the following two steps until convergence [40]:

• E-step: updating $q\left(\mathit{Z}\right)$ by setting it to the posterior of $\mathit{Z}$ given the current $\mathsf{\Theta }$, which is used to construct the local ELBO on the log-likelihood.
• M-step: updating $\mathsf{\Theta }$ by maximizing the ELBO, i.e., solving for $\mathsf{\Theta }$ in

  $$\frac{\partial }{\partial \mathsf{\Theta }}\sum _{\mathit{Z}}q\left(\mathit{Z}\right)lnp\left(\mathit{X},\mathit{Z}\mid \mathsf{\Theta }\right)=0.$$

  (3)

3.2.2. Network Generation Model

In this study, the parameters in the model for generating the topology of each layer include the degree sequences of layers $\mathit{d}$ and the adjacency of the aggregate topology from all layers of the multiplex network $\mathit{A}$. The degree sequence of a layer is used to generate the topology of that layer using the configuration model [42,43]. According to this model, for a network with a given degree sequence, the probability of a link between two nodes with degree ${\mathit{d}}_i$ and ${\mathit{d}}_j$ is $\frac{{\mathit{d}}_i{\mathit{d}}_j}{2\left|\mathcal{E}\right|-1}$ where $\left|\mathcal{E}\right|$ is the number of edges. Therefore, after estimating the degree sequences of all layers, we can probabilistically infer the complete topology of each layer. The adjacency matrix of the aggregate topology is included in the model parameters because we add an aggregation step to the EM framework to obtain the aggregated topology of all layers. In this way, information on the interlayer dependence of the multiplex network is leveraged to improve the reconstruction accuracy. Without this step, the topologies of layers are independent of each other and are hence reconstructed separately.

3.2.3. Full Algorithm for the EMA Framework

This section presents details about how to reconstruct multiplex networks using the EMA framework, including the E-step, M-step, and A-step.

In the E-step, the task is to obtain the posterior of $\mathit{Z}$ given by Equation (2). According to the configuration model, the probability of a link between nodes i and j in layer ℓ before it is not observed yet can be given by

$$p_{ij}^{\ell }=\frac{{\mathit{d}}_i^{\ell }{\mathit{d}}_j^{\ell }}{2|{\mathcal{E}}^{\ell }|-1}$$

(4)

Considering that $\frac{{\mathit{d}}_i^{\ell }{\mathit{d}}_j^{\ell }}{2|{\mathcal{E}}^{\ell }|-1}$ can be greater than 1 and certain entries are observed, Equation (4) is changed to

$$p_{ij}^{\ell }=\left\{\begin{array}{cc}{\mathit{X}}_{ij}^{\ell },\hfill & \mathrm{if}i\in {\mathcal{V}}_{\mathrm{obs}}\mathrm{and}j\in {\mathcal{V}}_{\mathrm{obs}}\hfill \\ min(1,\frac{{\mathit{d}}_i^{\ell }{\mathit{d}}_j^{\ell }}{2|{\mathcal{E}}^{\ell }|-1}),\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(5)

Then, $q\left({\mathit{Z}}_{ij}^{\ell }\right)$ is expressed as

$$q\left({\mathit{Z}}_{ij}^{\ell }\right)={\left(p_{ij}^{\ell }\right)}^{{\mathit{Z}}_{ij}^{\ell }}{\left(1-p_{ij}^{\ell }\right)}^{1-{\mathit{Z}}_{ij}^{\ell }},{\mathit{Z}}_{ij}^{\ell }\in \{0,1\},$$

(6)

Next, in the M-step, given the current $q\left({\mathit{Z}}_{ij}^{\ell }\right)$, we solve for ${\mathit{d}}_i^{\ell }$ for any $i\in \mathcal{N},\ell \in \mathcal{L}$ as follows. In the ELBO, the term involving ${\mathit{d}}_i^{\ell }$ is

$$\begin{array}{c}\hfill \sum _{{\mathit{Z}}_{ij}^{\ell }}q\left({\mathit{Z}}_{ij}^{\ell }\right)ln\left(\frac{{\left(\frac{{\mathit{d}}_i^{\ell }{\mathit{d}}_j^{\ell }}{2|{\mathcal{E}}^{\ell }|-1}\right)}^{{\mathit{Z}}_{ij}^{\ell }}{\left(1-\frac{{\mathit{d}}_i^{\ell }{\mathit{d}}_j^{\ell }}{2|{\mathcal{E}}^{\ell }|-1}\right)}^{1-{\mathit{Z}}_{ij}^{\ell }}}{q\left({\mathit{Z}}_{ij}^{\ell }\right)}\right)\end{array}$$

(7)

To simplify the term involving ${\mathit{d}}_i^{\ell }$ in Equation (7), we let $a_j=\frac{{\mathit{d}}_j^{\ell }}{2|{\mathcal{E}}^{\ell }|-1}$. Then differentiate with respect to an arbitrary ${\mathit{d}}_i^{\ell }$ and we obtain

$$\begin{array}{}\mathrm{(8a)}& \sum _{{\mathit{Z}}_{ij}^{\ell }}\frac{\partial }{\partial {\mathit{d}}_i^{\ell }}\left(q\left({\mathit{Z}}_{ij}^{\ell }\right)·ln\left[{\left(a_j{\mathit{d}}_i^{\ell }\right)}^{{\mathit{Z}}_{ij}^{\ell }}{(1-a_j{\mathit{d}}_i^{\ell })}^{1-{\mathit{Z}}_{ij}^{\ell }}\right]\right)\hfill \mathrm{(8b)}& =\sum _{{\mathit{Z}}_{ij}^{\ell }}q\left({\mathit{Z}}_{ij}^{\ell }\right)·\frac{{\mathit{Z}}_{ij}^{\ell }-a_j{\mathit{d}}_i^{\ell }}{{\mathit{d}}_i^{\ell }(1-a_j{\mathit{d}}_i^{\ell })}\hfill \mathrm{(8b)}& =\frac{\mathbb{E}\left[{\mathit{Z}}_{ij}^{\ell }\right]-a_j{\mathit{d}}_i^{\ell }}{{\mathit{d}}_i^{\ell }(1-a_j{\mathit{d}}_i^{\ell })}\hfill \end{array}$$

Let Equation (8c) = 0 and we get

$${\mathit{d}}_j^{\ell }{\mathit{d}}_i^{\ell }=\left(2|{\mathcal{E}}^{\ell }|-1\right)\mathbb{E}\left[{\mathit{Z}}_{ij}^{\ell }\right],\forall i\in \mathcal{N}.$$

(9)

To ensure that the solution to ${\mathit{d}}_i^{\ell }$ does not involve ${\mathit{d}}_j^{\ell }$, we sum both sides of the above equation over j and obtain

$${\mathit{d}}_i^{\ell }=\frac{2|{\mathcal{E}}^{\ell }|-1}{2|{\mathcal{E}}^{\ell }|-{\mathit{d}}_i^{\ell }}\sum _{j\in \mathcal{N}}\mathbb{E}\left[{\mathit{Z}}_{ij}^{\ell }\right],\forall i\in \mathcal{N}.$$

(10)

Solving for ${\mathit{d}}_i^{\ell }$, we have

$${\mathit{d}}_i^{\ell }=\left|{\mathcal{E}}^{\ell }\right|-\sqrt{|{\mathcal{E}}^{\ell }{|}^2-2\left|{\mathcal{E}}^{\ell }\right|\sum _{j\in \mathcal{N}}\mathbb{E}\left[{\mathit{Z}}_{ij}^{\ell }\right]+\sum _{j\in \mathcal{N}}\mathbb{E}\left[{\mathit{Z}}_{ij}^{\ell }\right]}.$$

(11)

Note that Equation (10) has two solutions to ${\mathit{d}}_i^{\ell }$, but the other solution is discarded because it is greater than the number of links in layer ℓ. For ${\mathit{d}}_i^{\ell }$ to have a real-valued solution, it is easy to show that the following condition should be satisfied

$${\left(|{\mathcal{E}}^{\ell }|-\sum _{j\in \mathcal{N}}\mathbb{E}\left[{\mathit{Z}}_{ij}^{\ell }\right]\right)}^2\ge {\left(\sum _{j\in \mathcal{N}}\mathbb{E}\left[{\mathit{Z}}_{ij}^{\ell }\right]\right)}^2-\sum _{j\in \mathcal{N}}\mathbb{E}\left[{\mathit{Z}}_{ij}^{\ell }\right],\forall i\in \mathcal{N},$$

(12)

Since $|{\mathcal{E}}^{\ell }|-{\sum }_{j\in \mathcal{N}}\mathbb{E}\left[{\mathit{Z}}_{ij}^{\ell }\right]\ge 0$, we have

$$\frac{|{\mathcal{E}}^{\ell }|}{{\sum }_{j\in \mathcal{N}}\mathbb{E}\left[{\mathit{Z}}_{ij}^{\ell }\right]}\ge \sqrt{\frac{{\sum }_{j\in \mathcal{N}}\mathbb{E}\left[{\mathit{Z}}_{ij}^{\ell }\right]-1}{{\sum }_{j\in \mathcal{N}}\mathbb{E}\left[{\mathit{Z}}_{ij}^{\ell }\right]}}+1,\forall i\in \mathcal{N},$$

(13)

which indicates that a layer should not have giant hub nodes that are incident with over 50% of the edges of a layer. This condition is typically satisfied for sparse multiplex networks. For layers with a sufficiently large number of edges but no giant hub node, $|{\mathcal{E}}^{\ell }|\gg {\mathit{d}}_i^{\ell }$, thus the solution to Equation (10) can be approximately given by

$${\mathit{d}}_i^{\ell }\approx \sum _{j\in \mathcal{N}}\mathbb{E}\left[{\mathit{Z}}_{ij}^{\ell }\right]=\sum _{j\in \mathcal{N}}{p}_{ij}^{\ell },\forall i\in \mathcal{N}.$$

(14)

In the A-step, the aggregate topology from topologies of all layers in the multiplex network is obtained using the OR aggregation mechanism [22]. Specifically, the respective entries in the aggregate adjacency matrix should be equal to those calculated from the observed part of the adjacency matrix for each layer, which is expressed as

$${\mathit{A}}_{ij}=\phi ({\mathit{Z}}_{ij}^1,\cdots ,{\mathit{Z}}_{ij}^{\left|\mathcal{L}\right|})=1-\prod _{\ell \in \mathcal{L}}(1-{\mathit{Z}}_{ij}^{\ell })$$

(15)

where $\phi (·)$ is the aggregation function. Apparently, ${\mathit{A}}_{ij}=1$ when there exists at least one ℓ such that ${\mathit{Z}}_{ij}^{\ell }=1$ and ${\mathit{A}}_{ij}=0$ only when ${\mathit{Z}}_{ij}^{\ell }=0$ for all ℓ. The aggregate topology (where nodes are connected) is used to update the estimated complete topology of each layer as follows. According to the Bayes’ Theorem, we have

$$\begin{array}{}\mathrm{(16a)}& \hfill P({\mathit{Z}}_{ij}^{\ell }=1\mid {\mathit{A}}_{ij}=1)& =\frac{P({\mathit{A}}_{ij}=1\mid {\mathit{Z}}_{ij}^{\ell }=1)·P({\mathit{Z}}_{ij}^{\ell }=1)}{P({\mathit{A}}_{ij}=1)}\hfill \mathrm{(16b)}& \hfill & =\frac{1·p_{ij}^{\ell }}{1-\prod _{h\in \mathcal{L}}(1-p_{ij}^h)}\hfill \end{array}$$

In the implementation, this means that the latest $p_{ij}^{\ell }$ for all $i,j$, which is equivalent to $P({\mathit{Z}}_{ij}^{\ell }=1\mid {\mathit{A}}_{ij}=1)$ at this point, is updated according to Equation (16b).

The complete EMA algorithm is summarized in Algorithm 1. Note that we can treat the aggregate topology as part of the parameters for generating the multiplex network, i.e., $\mathsf{\Theta }=(\mathit{A},\mathit{d})$. From this perspective, the EMA framework is reduced to an EM framework and the A-step is regarded as a component of the M-step for estimating the optimal model parameters. This is because given the current estimate of the unobserved part of the multiplex network, the aggregate topology obtained by Equation (15) is the only solution and hence the optimal solution to the aggregate topology.

Algorithm 1: EMA for multiplex network reconstruction

Input: Error tolerance ${ϵ}_T$, maximum number of iterations $Ite{r}_{max}$, number of unique nodes in the multiplex network $\left|\mathcal{N}\right|$, and partially observed topology $\mathit{X}$. Output: Reconstructed complete topology $\mathit{Z}$ for the multiplex network. 1: Initialize the degree sequences and predicted adjacency matrices: ${\mathit{d}}_i^{\ell }\sim U(1,|\mathcal{N}\left|\right),\forall i\in \{1,\dots ,|\mathcal{N}\left|\right\}$ and ${\mathit{Z}}_{ij}^{\ell }\sim U(0,1),\forall i\in \{1,\dots ,|\mathcal{N}\left|\right\}$. 2: for$iter$ in 1 to $Ite{r}_{max}$ do 3:     E-step: Update $q\left(\mathit{Z}\right)$ according to Equation (6). 4:     A-step: Update $\mathit{Z}$ using the aggregate topology according to Equation (16). 5:     M-step: Update $\mathit{d}$ according to Equation (14). 6:     Calculate the error $ϵ$ 7:     if $ϵ>{ϵ}_T$ then Continue 8:     else Break 9:     end if 10: end for

4. Datasets

This section presents the details about the real-world multiplex network we use in validating the EMA framework. Covert multiplex networks include the drug trafficking multiplex network and the Sicilian Mafia multiplex network. To demonstrate that the EMA framework is general, we validate it against the other two models on overt multiplex networks, including the London transportation multiplex network and the C. elegans neural multiplex network. This set of networks provide an opportunity to test the EMA framework under varying fractions of observed components of multiplex networks and permit the comparison between social networks, biological, and physical networks.

4.1. Drug Trafficking Multiplex Network

The dataset on drug trafficking networks in this study was obtained from the 2007 ‘E’ Division Provincial Threat Assessment (PTA) report generated by Criminal Intelligence Service Canada and the Royal Canadian Mounted Police (RCMP) [39]. This report contains information about all individuals associated with the criminal operations of 129 different groups. This extensive data collection effort maps a community of drug trafficking operations, revealing how individuals interact within and between groups to sustain regional illicit drug supply and distribution. Individual-level information includes demographic characteristics, a description of three years of drug trafficking activity, their role in drug trafficking, and their current and historical relationships among each other, including co-offenders (individuals who commit crimes together), legitimate business partners (including friends), their kinship, and involvement in a formal criminal organization (including enemies) (Figure 2). This multiplex network integrates information obtained from official records management systems, surveillance, and non-criminal justice information sources. The layers for co-offenders and legitimate business partners are used to build a two-layer multiplex network (representing surveillance data capturing current activity, criminal or otherwise), the layers for co-offenders, legitimate business partners, and formal criminal organization are used to build a three-layer multiplex network (adding social embedding in a community of criminal actors), and the four-layer multiplex network adds familial relations (adding trusted relations that bolster network resiliency). It is obvious from the figure that both layers do not have giant hub nodes that are incident with over 50% of the edges of a layer. Some key network statistics of each of the four layers are presented in

Table 1. Network statistics of each layer in the drug trafficking multiplex network. FCO is the shorthand for formal criminal organization.

Layer	$\left|\mathcal{N}\right|$	$\left|\mathcal{E}\right|$	$\mathit{\rho }(\times {10}^{-3})$	$\langle \mathit{d}\rangle$	$\langle \left|\mathit{CC}\right|\rangle$	$\left|\mathit{GCC}\right|$	CoV of $\left|\mathit{CC}\right|$
Co-offender	1645	1808	1.32	2.19	19.00	1024	5.71
Legitimate	1022	1041	2.00	2.04	10.99	462	4.43
FCO	560	597	3.81	2.13	14.36	150	2.24
Kinship	399	308	3.88	1.54	3.22	25	0.99

. We can see from the table that each layer has a low link density ($\rho$), low average degree ($\langle d\rangle$), and the average size of the connected component ($\langle \left|CC\right|\rangle$). Although the size of the greatest connected component is not small, particularly the one in the co-offender network, the coefficient of variance (CoV, i.e., ratio of standard deviation and mean) of the size of the connected component is small, indicating that members in the drug trafficking community are generally isolated.

4.2. Sicilian Mafia Multiplex Network

The Sicilian Mafia multiplex network (Figure 3) was derived from the pre-trial detention order issued by the local court in 2007 in response to a major investigation [44]. The Italian civil law system is substantively different than federalist case law systems (such as the U.S. system), and as such, the court documents filed in support of the detention order provide sufficient information to map individual-level involvement in criminal enterprise activities [15]. Network information captures five years of observations (2003 to 2007) for two Mafia clans, known as the Mistretta family and the Batanesi clan. These familial clans are embedded within a community of illicit entrepreneurs collectively referred to as the Sicilian Mafia. This data set includes two types of layers: (i) meeting data among suspected (101 nodes and 256 edges) and (ii) phone calls among individuals (100 nodes and 124 edges). These layers represent different modes of communication commonly observed during surveillance activity in support of major investigations. Some key statistics of each layer are presented in

Table 2. Network statistics of each layer in the Sicilian Mafia multiplex network.

Layer	$\left|\mathcal{N}\right|$	$\left|\mathcal{E}\right|$	$\mathit{\rho }(\times {10}^{-3})$	$\langle \mathit{d}\rangle$	$\langle \left|\mathit{CC}\right|\rangle$	$\left|\mathit{GCC}\right|$	CoV of $\left|\mathit{CC}\right|$
Meeting	101	256	55.77	5.07	19.02	86	1.76
Call	100	124	27.45	2.48	23.50	85	1.51

. It can be observed from the figure that both layers do not have giant hub nodes as well.

4.3. London Transportation Multiplex Network

The London transportation multiplex network [45] (Figure 4) consists of underground (260 nodes and 225 edges), overground (81 nodes and 62 edges), and lightrail (44 nodes and 34 edges). The overground and light rail layers are selected to construct a two-layer multiplex network. The third layer adds the underground edges. While the modalities (layers) of transit systems are designed to intersect forming a fully connected giant component, single layers will not always form a connected component. Moreover, this network is intentionally incomplete in that buses that may link different modalities (e.g., underground to light rail) were not included. Some key statistics of each layer are presented in

Table 3. Network statistics of each layer in the London transportation multiplex network.

Layer	$\left|\mathcal{N}\right|$	$\left|\mathcal{E}\right|$	$\mathit{\rho }(\times {10}^{-3})$	$\langle \mathit{d}\rangle$	$\langle \left|\mathit{CC}\right|\rangle$	$\left|\mathit{GCC}\right|$	CoV of $\left|\mathit{CC}\right|$
Underground	260	225	6.68	1.73	5.65	98	2.48
Overground	81	62	19.14	1.53	4.26	17	0.83
Lightrail	44	34	35.94	1.55	4.00	8	0.54

. Each layer of this transportation multiplex network has many connected components and the size of the greatest connected component is not very large compared to the number of nodes in each layer.

4.4. C. elegans Multiplex Network

The C. (Caenorhabditis) elegans neural multiplex network [46] (Figure 5) consists of three layers for different types of synaptic junctions, including electric (242 nodes and 450 edges), chemical (260 nodes and 869 edges), and polyadic (277 nodes and 1666 edges). Layers for chemical and polyadic junctions are selected to construct a two-layer multiplex network. Some key statistics of each layer are presented in

Table 4. Network statistics of each layer in the C. elegans multiplex network.

Layer	$\left|\mathcal{N}\right|$	$\left|\mathcal{E}\right|$	$\mathit{\rho }(\times {10}^{-3})$	$\langle \mathit{d}\rangle$	$\langle \left|\mathit{CC}\right|\rangle$	$\left|\mathit{GCC}\right|$	CoV of $\left|\mathit{CC}\right|$
Electrical	242	450	15.43	3.72	80.67	237	1.37
Chemical	260	869	25.81	6.68	130.00	258	0.98
Polyadic	277	1666	43.58	12.03	277.00	277	0.00

. Compared to the covert multiplex networks introduced earlier, all layers of the C. elegans multiplex network are much denser because each layer has a higher density. In each layer, the size of the greatest connected component is large. In particular, all the nodes in the polyadic layer are connected because the size of the greatest connected component is equal to the number of nodes.

5. Experiments

5.1. Experiment Setup

For each multiplex network, we use 50 random simulations at each fraction of missing nodes to obtain the average predictive performance (we find that 50 repetitions can already achieve a stable mean value). When a node is not observed, the links connected to this node are also considered unobserved. For both EMA and EM, the error tolerance ${ϵ}_T$ is set to $1\times {10}^{-5}$, i.e., the solution is considered converged if the difference in the mean absolute error (MAE) of the predicted adjacency matrices (unobserved entries only) at two consecutive iterations for all unobserved links is lower than the tolerance.

Because covert networks are sparse, the associated datasets are typically imbalanced toward negatives (no links between nodes), we employ the Matthews Correlation Coefficient (MCC) and G-mean to measure the prediction (classification) accuracy [47,48]. MCC can be computed as

$$\mathrm{MCC}=\frac{\mathrm{TP}\times \mathrm{TN}-\mathrm{FP}\times \mathrm{FN}}{\sqrt{(\mathrm{TP}+\mathrm{FP})\times (\mathrm{TP}+\mathrm{FN})\times (\mathrm{TN}+\mathrm{FP})\times (\mathrm{TN}+\mathrm{FN})}}$$

(17)

where TP represents true positives, TN represents true negatives, FP represents false positives, and FN represents false negatives. Note that the value of MCC can be negative. G-mean is given by

$$\mathrm{G}-\mathrm{mean}=\sqrt{\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}\times \frac{\mathrm{TN}}{\mathrm{TN}+\mathrm{FP}}}=\sqrt{\mathrm{Recall}\times \mathrm{Specificity}}$$

(18)

We also include the random model (RM), i.e., a binary uniformly random classifier, as a baseline model. Among all possible links, except for the links among observed nodes, the RM selects a number (i.e., the number of links left) of links in a uniformly random manner as the predicted links.

5.2. Results

We first check the performance of EMA under different fractions of missing components. The predictive performances measured by G-mean and MCC for different multiplex networks under 10% to 90% of missing components are presented in Figure 6, Figure 7, Figure 8 and Figure 9 where c represents the fraction of observed components. We can see that EMA achieves the best MCC regardless of the multiplex networks and the fraction of missing components. In terms of G-mean, EMA outperforms EM in all cases and random model when the observed fraction of components is large. Note that the trends of G-mean and MCC over different fractions of observed components obtained by EMA and EM are close to each other because the values of FP × FN are much smaller than TP × TN.

Next, we examine the impact of the number of layers on the EMA and EM approaches. As the number of layers increases from two to three in drug trafficking (Figure 6), C. elegans (Figure 8), and London transport networks (Figure 9), the performance improvement of EMA when compared to EM in both G-mean and MCC decreases, indicating that the benefits of adding the aggregation step drops. This is because when there are more layers, there are more choices of values for individual entries in the layers, making a positive entry in the aggregate topology less informative about the positive value of each layer.

We also compare the convergence of EMA and EM algorithms on the drug trafficking multiplex networks (Figure 10) and London transportation multiplex networks (Figure 11). Regardless of the size of multiplex networks and the fraction of components observed, EM converges more quickly than EMA. In the worst-case tested, EMA can reach convergence within 40 iterations given a convergence error of $1\times {10}^{-5}$. However, if the convergence error can be set to a larger value, such as $1\times {10}^{-4}$, then EMA can achieve convergence much faster (within 10 iterations across the networks considered). Therefore, both algorithms can be applied to moderately large-scale multiplex networks.

6. Discussion

From a criminal justice perspective, mapping networks using available data on illicit activities is typically not sufficient to obtain the complete structure of the covert networks under study. As a result, networks of criminal operations are always incomplete, thereby compromising the potential practical contributions network science can make toward developing effective interdiction strategies. Two major implications for Disrupting Illicit Networks arise from this investigation for reconstructing a more complete structure of multiplex networks.

First, with regard to the broad context of disrupting covert networks, the experiments support arguments that criminal enterprise is intimately embedded within a complex social system wherein a pair of individuals can be connected by multiple relations or a criminal actor could tap into different social networks to resolve a pressing issue [35,49,50]. For example, a drug smuggler uses legitimate business assets (helicopter and employed pilot) to transport illicit products (multiple relations) and, in a pinch, the smuggler could launder proceeds through a relative’s real-estate business (temporary use of a different type of relation for criminal ends). The extent to which associates and kin are aware of these illegalities is not necessarily relevant; what is critical is that this interwoven social fabric must be fully understood before effective interdiction strategies can be developed and implemented. Results on the performance of the EMA framework for the drug trafficking network and the Sicilian mafia network show that networks generated from the surveillance of current criminal activity are insufficient. Interpolating latent connections by adding at least one more layer that captures historic criminogenic relations (i.e., formal organized crime membership), business or professional associations, or purely social and familial connections significantly improves network completeness. This finding coincides with applied network criminology [34,39] and crime opportunity theory [51,52].

Second, while inference can improve network completeness, such modeling frameworks do not constitute a silver bullet for all methodological issues. Continued efforts are needed to improve data collection protocols to enhance the capacity to generate more complete information about covert operations and the complex social systems within which they are embedded. To this aim, our experiments show that identifying additional layers composing a multiplex network makes the EMA algorithm more accurate. Supporting prior arguments to establish methodological conventions regarding the use of criminal justice system data in network applications [53,54], this finding confirms that the key to uncovering more of the covert multiplex network is to tap different layers of information sources [34,51]. However, additional tests are needed to understand the nature of information benefits. For example, how does the dependency between layers, particularly in social networks, impact information benefits? When does identifying additional layers of covert networks benefit more than identifying more components within known layers? Finally, given that the covert networks were mapped from observations made over many years, does the model offer more information benefits when used with shorter observations, such as six months?

7. Conclusions

In this study, we develop the EMA framework for reconstructing multiplex networks and validate this approach against the EM approach and the random model on several multiplex networks given different fractions of the observed part of the complete networks. We find that EMA achieves the best predictive performance in terms of G-mean and MCC over different levels of available observed parts of the true networks. The findings are consistent across experiments conducted on different types of multiplex networks (social and non-social). This indicates that integrating the estimation of the aggregate topology, which encodes the interlayer dependency of multiplex networks, can improve the reconstruction accuracy. By inferring a more complete structure of covert multiplex networks, the EMA framework can support interdiction efforts aimed at disrupting criminal operations, especially when the observations are only available for two or three types of interactions. As the fraction of observed actors and the type of interactions among actors increases, the added benefits of the EMA framework when compared to the EM framework decrease. In cases where quick decisions are required for interdicting multiplex with many layers, the EM framework is recommended due to a lower computational load.

Future work will consider the use of graph neural networks to reconstruct sparse multiplex networks. Although GNN has been widely applied to predict missing links [23,55,56], research efforts that apply GNN to reconstruct multiplex networks are limited. As such, it is valuable to compare the reconstruction accuracy of EMA and GNN on multiple datasets and identify where each approach performs the best. Furthermore, since the accuracy of the EMA framework reduces for multiplex networks with many layers, future research can explore how to improve the predictive accuracy using network generative models with a prescribed level of interlayer dependency [57,58,59,60], and use the dynamics of the networks [61].