Network reconstruction and community detection from dynamics

TL;DR

This paper introduces a scalable Bayesian approach that simultaneously reconstructs network structures and detects communities from observed dynamics, improving accuracy through their mutual reinforcement.

Contribution

It presents a novel nonparametric Bayesian method that jointly infers network topology and community structure from functional data, enhancing both tasks.

Findings

Effective on synthetic and empirical networks

Improves reconstruction accuracy via community detection

Applicable to epidemic and Ising model data

Abstract

We present a scalable nonparametric Bayesian method to perform network reconstruction from observed functional behavior that at the same time infers the communities present in the network. We show that the joint reconstruction with community detection has a synergistic effect, where the edge correlations used to inform the existence of communities are also inherently used to improve the accuracy of the reconstruction which, in turn, can better inform the uncovering of communities. We illustrate the use of our method with observations arising from epidemic models and the Ising model, both on synthetic and empirical networks, as well as on data containing only functional information.

Figures

Tables

Table 1. Table 1: Country memberships of twitter users, according to the groups inferred by the reconstruction method.

Group	User countries
1	$\text{Brazil}\times 216$, $\text{Japan}\times 2$, $\text{Germany}\times 2$, $\text{USA}\times 1$, $\text{Russia}\times 1$, $\text{Argentina}\times 1$, $\text{Colombia}\times 1$, $\text{France}\times 1$, $\text{South Africa}\times 1$
2	$\text{Japan}\times 104$, $\text{UK}\times 1$, $\text{China}\times 1$, $\text{Italy}\times 1$, $\text{USA}\times 1$
3	$\text{Indonesia}\times 42$, $\text{India}\times 10$, $\text{UK}\times 8$, $\text{Germany}\times 6$, $\text{USA}\times 6$, $\text{Thailand}\times 1$, $\text{Brazil}\times 1$, $\text{Bali}\times 1$, $\text{Argentina}\times 1$, $\text{Australia}\times 1$
4	$\text{Indonesia}\times 52$, $\text{Germany}\times 6$, $\text{UK}\times 4$, $\text{USA}\times 3$, $\text{Russia}\times 2$, $\text{Australia}\times 2$, $\text{New Zealand}\times 1$, $\text{India}\times 1$, $\text{Botswana}\times 1$, $\text{Brazil}\times 1$
5	$\text{USA}\times 54$, $\text{UK}\times 4$, $\text{Russia}\times 1$, $\text{Japan}\times 1$, $\text{Indonesia}\times 1$
6	$\text{USA}\times 23$, $\text{Brazil}\times 9$, $\text{UK}\times 6$, $\text{Netherlands}\times 3$, $\text{Russia}\times 2$, $\text{Canada}\times 2$, $\text{Italy}\times 2$, $\text{Germany}\times 2$, $\text{Mexico}\times 1$, $\text{South Africa}\times 1$, $\text{Kwait}\times 1$, $\text{Austria}\times 1$, $\text{Romania}\times 1$, $\text{Finland}\times 1$, $\text{Japan}\times 1$, $\text{Philippines}\times 1$, $\text{Egypt}\times 1$, $\text{Argentina}\times 1$, $\text{Chile}\times 1$
7	$\text{USA}\times 18$, $\text{India}\times 8$, $\text{Portugal}\times 2$, $\text{Brazil}\times 2$, $\text{Japan}\times 2$, $\text{UK}\times 2$, $\text{Guernsey}\times 1$, $\text{Malaysia}\times 1$, $\text{Mexico}\times 1$, $\text{Australia}\times 1$, $\text{Russia}\times 1$, $\text{France}\times 1$, $\text{Vietnam}\times 1$, $\text{Spain}\times 1$, $\text{Venezuela}\times 1$, $\text{Philippines}\times 1$
8	$\text{Japan}\times 35$
9	$\text{Brazil}\times 31$, $\text{USA}\times 2$, $\text{UK}\times 1$
10	$\text{Korea}\times 24$, $\text{USA}\times 1$, $\text{Argentina}\times 1$, $\text{Russia}\times 1$, $\text{Japan}\times 1$
11	$\text{Indonesia}\times 22$, $\text{UK}\times 2$, $\text{Australia}\times 1$, $\text{India}\times 1$, $\text{USA}\times 1$
12	$\text{USA}\times 11$, $\text{Japan}\times 2$, $\text{Germany}\times 2$, $\text{France}\times 1$, $\text{Korea}\times 1$, $\text{Thailand}\times 1$, $\text{Brazil}\times 1$, $\text{Chile}\times 1$, $\text{Indonesia}\times 1$
13	$\text{USA}\times 4$, $\text{Indonesia}\times 4$, $\text{India}\times 4$, $\text{UK}\times 3$, $\text{Australia}\times 2$, $\text{Canada}\times 1$
14	$\text{Japan}\times 8$, $\text{Venezuela}\times 6$, $\text{USA}\times 1$, $\text{Chile}\times 1$
15	$\text{Indonesia}\times 14$, $\text{Japan}\times 1$
16	$\text{Indonesia}\times 7$, $\text{USA}\times 3$, $\text{Turkey}\times 1$, $\text{Philippines}\times 1$, $\text{Brazil}\times 1$
17	$\text{USA}\times 9$, $\text{Canada}\times 1$, $\text{France}\times 1$, $\text{Belgium}\times 1$
18	$\text{Germany}\times 5$, $\text{USA}\times 1$, $\text{Russia}\times 1$, $\text{France}\times 1$
19	$\text{Japan}\times 4$

Table 2. Table 2: Party affiliation of deputies of the lower house of the Brazilian congress, according to the groups they were classified by the reconstruction method. Parties in red belong to the center-left government coalition, and in blue to the right-wing opposition.

Group	Parties
1	$\text{PMDB}\times 83$, $\text{PT}\times 70$, $\text{PP}\times 38$, $\text{PR}\times 36$, $\text{PSB}\times 27$, $\text{PDT}\times 22$, $\text{PTB}\times 17$, $\text{PV}\times 14$, $\text{PCdoB}\times 12$, $\text{PSC}\times 10$, $\text{DEM}\times 5$, $\text{PMN}\times 5$, $\text{PRB}\times 3$, $\text{PSOL}\times 3$, $\text{PHS}\times 2$, $\text{PTdoB}\times 1$, $\text{PTC}\times 1$
2	$\text{PSDB}\times 54$, $\text{PPS}\times 8$, $\text{PFL}\times 2$
3	$\text{DEM}\times 41$
4	$\text{DEM}\times 8$, $\text{PMDB}\times 6$, $\text{PPS}\times 4$, $\text{PP}\times 2$, $\text{PSB}\times 1$

Equations

$$P(\mathbf{A}|\mathcal{D})=\frac{P(\mathcal{D}|\mathbf{A})P(\mathbf{A})}{P(\mathcal{D})},$$

$$P(\mathbf{A},\mathbf{b}|\mathcal{D})=\frac{P(\mathcal{D}|\mathbf{A})P(\mathbf{A}|\mathbf{b})P(\mathbf{b})}{P(\mathcal{D})}.$$

$$P(\mathbf{A}|\mathbf{\lambda },\mathbf{\kappa },\mathbf{b})=\prod _{i<j}\frac{e^{-{\kappa }_i{\kappa }_j{\lambda }_{b_i,b_j}}({\kappa }_i{\kappa }_j{\lambda }_{b_i,b_j}{)}^{A_{ij}}}{A_{ij}!},$$

$$P(\mathbf{A}|\mathbf{b})=\int P(\mathbf{A}|\mathbf{\lambda },\mathbf{\kappa },\mathbf{b})P(\mathbf{\kappa }|\mathbf{b})P(\mathbf{\lambda }|\mathbf{b})\mathrm{d}\mathbf{\kappa }\mathrm{d}\mathbf{\lambda }.$$

$$P(\mathbf{\sigma }|\mathbf{A},\mathbf{\tau },\gamma )=\prod _t\prod _{i}P({\sigma }_i(t)|\mathbf{\sigma }(t-1)),$$

$$\begin{gathered} P({\sigma }_i(t)|\mathbf{\sigma }(t-1))= \\ f({\mathrm{e}}^{m_i(t-1)},{\sigma }_i(t){)}^{1-{\sigma }_i(t-1)}\times f(\gamma ,{\sigma }_i(t){)}^{{\sigma }_i(t-1)} \end{gathered}$$

$$P(\mathbf{A},\mathbf{b},\mathbf{\tau }|\mathbf{\sigma })=\frac{P(\mathbf{\sigma }|\mathbf{A},\mathbf{\tau },\gamma )P(\mathbf{A}|\mathbf{b})P(\mathbf{b})P(\mathbf{\tau })}{P(\mathbf{\sigma }|\gamma )}.$$

$$P(\mathbf{s}|\mathbf{A},\beta ,\mathbf{J},\mathbf{h})=\frac{\exp \left(\beta {\sum }_{i<j}{J}_{ij}{A}_{ij}{s}_i{s}_j+{\sum }_i{h}_i{s}_i\right)}{Z(\mathbf{A},\beta ,\mathbf{J},\mathbf{h})},$$

$$P(\mathbf{s}|\mathbf{A},\beta ,\mathbf{J},\mathbf{h})=\prod _i\frac{\exp (\beta s_i{\sum }_j{J}_{ij}{A}_{ij}{s}_j+h_i{s}_i)}{2\cosh (\beta {\sum }_j{J}_{ij}{A}_{ij}{s}_j+h_i)}.$$

$$\begin{gathered} P(\mathbf{A},\mathbf{b},\beta ,\mathbf{J},\mathbf{h}|\bar{\mathbf{s}})= \\ \frac{P(\bar{\mathbf{s}}|\mathbf{A},\beta ,\mathbf{J},\mathbf{h})P(\beta )P(\mathbf{h})P(\mathbf{J}|\mathbf{A})P(\mathbf{A}|\mathbf{b})P(\mathbf{b})}{P(\bar{\mathbf{s}})}. \end{gathered}$$

$$\min \left(1,\frac{P({\mathbf{A}}^{\prime },\mathbf{b},\mathbf{\theta }|\mathcal{D})}{P(\mathbf{A},\mathbf{b},\mathbf{\theta }|\mathcal{D})}\frac{P({\mathbf{A}}^{\prime }\to \mathbf{A})}{P(\mathbf{A}\to {\mathbf{A}}^{\prime })}\right)$$

$$S({\mathbf{A}}^{\ast },\mathbf{\pi })=1-\frac{{\sum }_{i<j}|A_{ij}^{\ast }-{\pi }_{ij}|}{{\sum }_{i<j}|A_{ij}^{\ast }+{\pi }_{ij}|},$$

$${\hat{A}}_{ij}=\{\begin{array}{ll}1& \text{if }{\pi }_{ij}>1/2,\\ 0& \text{if }{\pi }_{ij}<1/2.\end{array}$$

$$P(\mathbf{A}|\mathbf{b})$$

$$=P(\mathbf{A}|\mathbf{k},\mathbf{e},\mathbf{b})P(\mathbf{k}|\mathbf{e},\mathbf{b})P(\mathbf{e}|\mathbf{b}),$$

$$P(\mathbf{A}|\mathbf{k},\mathbf{e},\mathbf{b})$$

$$P(\mathbf{k}|\mathbf{e},\mathbf{b})$$

$$P(\mathbf{e}|\mathbf{b})$$

$$P(\mathbf{k}|\mathbf{e},\mathbf{b})=P(\mathbf{k}|\mathbf{\eta })P(\mathbf{\eta }|\mathbf{e},\mathbf{b})$$

$$P(\mathbf{k}|\mathbf{\eta })=\prod _r\frac{{\prod }_k{\eta }_k^r!}{n_r!}$$

$$P(\mathbf{\eta }|\mathbf{e},\mathbf{b})=\prod _{r}q(e_r,n_r{)}^{-1}$$

$$q(m,n)=q(m,n-1)+q(m-n,n).$$

$$P(\mathbf{b})=P(\mathbf{b}|\mathbf{n})P(\mathbf{n}|B)P(B)=\frac{{\prod }_r{n}_r!}{N!}{\left(\genfrac{}{}{0ex}{}{N-1}{B-1}\right)}^{-1}{N}^{-1}.$$

$$P({\mathbf{e}}_l|{\mathbf{e}}_{l+1},{\mathbf{b}}_l)=\prod _{r<s}{\left(\left(\genfrac{}{}{0ex}{}{n_r^l{n}_s^l}{e_{rs}^{l+1}}\right)\right)}^{-1}\prod _r{\left(\left(\genfrac{}{}{0ex}{}{n_r^l(n_r^l+1)/2}{e_{rr}^{l+1}/2}\right)\right)}^{-1},$$

$$A_{ij}(G_{ij})=\{\begin{array}{ll}1& \text{ if }{G}_{ij}>0\text{ and }i\ne j,\\ 0& \text{ otherwise.}\end{array}$$

$$P(\mathbf{A}|\theta )$$

$$=\sum _{\mathbf{G}}P(\mathbf{A}|\mathbf{G})P(\mathbf{G}|\theta ),$$

$$P(\mathbf{A}|\mathbf{G})=\{\begin{array}{ll}1& \text{ if }\mathbf{A}=\mathbf{A}(\mathbf{G}),\\ 0& \text{ otherwise.}\end{array}$$

$$P(\mathbf{A},\mathbf{b}|\mathcal{D})=\frac{P(\mathcal{D}|\mathbf{A})P(\mathbf{A},\mathbf{b})}{P(\mathcal{D})},$$

$$P(\mathbf{A},\mathbf{G},\mathbf{b}|\mathcal{D})=\frac{P(\mathcal{D}|\mathbf{A})P(\mathbf{A}|\mathbf{G})P(\mathbf{G},\mathbf{b})}{P(\mathcal{D})},$$

$$P(\mathbf{A},\mathbf{b}|\mathcal{D})=\sum _{\mathbf{G}}P(\mathbf{A},\mathbf{G},\mathbf{b}|\mathcal{D}),$$

$$P(\mathbf{A},\mathbf{b}|\mathcal{D})=\frac{P(\mathcal{D}|\mathbf{A})P(\mathbf{A}|\mathbf{b})P(\mathbf{b})}{P(\mathcal{D})},$$

Full text

The observed functional behavior of a wide variety large-scale system is often the result of a network of pairwise interactions. However, in many cases these interactions are hidden from us, either because they are impossible to measure directly, or because their measurement can be done only at significant experimental cost. Examples include the mechanisms of gene and metabolic regulation Wang et al. (2006), brain connectivity Breakspear (2017), the spread of epidemics Keeling and Rohani (2002), systemic risk in financial institutions Musmeci et al. (2013), and influence in social media Bakshy et al. (2012). In such situations, we are required to infer the network of interactions from the observed functional behavior. Researchers have approached this reconstruction task from a variety of angles, resulting in many different methods, including thresholding the correlation between time-series Kramer et al. (2009), inversion of deterministic dynamics Timme (2007); Shandilya and Timme (2011); Nitzan et al. (2017), statistical inference of graphical models Abbeel et al. (2006); Bresler et al. (2008); Montanari and Pereira (2009); Höfling and Tibshirani (2009); Nguyen et al. (2017) and of models of epidemic spreading Gomez Rodriguez et al. (2010); Myers and Leskovec (2010); Netrapalli and Sanghavi (2012); Ma et al. (2018); Prasse and Van Mieghem (2018); Braunstein Alfredo et al. (2019), as well as approaches that avoid explicit modeling, such as those based on transfer entropy Runge et al. (2012), Granger causality Sun et al. (2015), compressed sensing Shen et al. (2014); Ma et al. (2015); Han et al. (2015), generalized linearization Li et al. (2017), and matching of pairwise correlations Ching et al. (2015); Lai (2017).

In this work, we approach the problem of network reconstruction in a manner that is different from the aforementioned methods in two important ways. First, we employ a nonparametric Bayesian formulation of the problem, which yields a full posterior distribution of possible networks that are compatible with the observed dynamical behavior. Second, we perform network reconstruction jointly with community detection Fortunato and Hric (2016), where at the same time as we infer the edges of the underlying network, we also infer its modular structure Peixoto (2017a). As we will show, while network reconstruction and community detection are desirable goals on their own, joining these two tasks has a synergistic effect, whereby the detection of communities significantly increases the accuracy of the reconstruction, which in turn improves the discovery of the communities, when compared to performing these tasks in isolation.

Some other approaches combine community detection with functional observation. Berthet et al. Berthet et al. (2016) derived necessary conditions for the exact recovery of group assignments for dense weighted networks generated with community structure given observed microstates of an Ising model. Hoffmann et al. Hoffmann et al. (2018) proposed a method to infer community structure from time-series data that bypasses network reconstruction, by employing instead a direct modeling of the dynamics given the group assignments. However, neither of these approaches attempt to perform network reconstruction together with community detection. Furthermore, they are tied down to one particular inverse problem, and as we will show, our general approach can be easily extended to an open-ended variety of functional models.

Bayesian network reconstruction — We approach the network reconstruction task similarly to the situation where the network edges are measured directly, but via an uncertain process Newman (2018); Peixoto (2018): If $\bm{\mathcal{D}}$ is the measurement of some process that takes place on a network, we can define a posterior distribution for the underlying adjacency matrix $\bm{A}$ via Bayes’ rule,

[TABLE]

where $P(\bm{\mathcal{D}}|\bm{A})$ is an arbitrary forward model for the dynamics given the network, $P(\bm{A})$ is the prior information on the network structure, and $P(\bm{\mathcal{D}})=\sum_{\bm{A}}P(\bm{\mathcal{D}}|\bm{A})P(\bm{A})$ is a normalization constant comprising the total evidence for the data $\bm{\mathcal{D}}$. We can unite reconstruction with community detection via an, at first, seemingly minor, but ultimately consequential modification of the above equation, where we introduce a structured prior $P(\bm{A}|\bm{b})$ where $\bm{b}$ represents the partition of the network in communities, i.e. $\bm{b}=\{b_{i}\}$, where $b_{i}\in\{1,\dots,B\}$ is group membership of node $i$. This partition is unknown, and is inferred together with the network itself, via the joint posterior distribution

[TABLE]

The prior $P(\bm{A}|\bm{b})$ is an assumed generative model for the network structure. In our work, we will use the degree-corrected stochastic block model (DC-SBM) Karrer and Newman (2011), which assumes that, besides differences in degree, nodes belonging to the same group have statistically equivalent connection patterns, according to the joint probability

[TABLE]

with $\lambda_{rs}$ determining the average number of edges between groups $r$ and $s$ and $\kappa_{i}$ the average degree of node $i$. The marginal prior is obtained by integrating over all remaining parameters weighted by their respective prior distributions,

[TABLE]

which can be computed exactly for standard prior choices, although it can be modified to include hierarchical priors that have an improved explanatory power Peixoto (2017b) (see Appendix A for a concise summary).

The use of the DC-SBM as a prior probability in Eq. 2 is motivated by its ability to inform link prediction in networks where some fraction of edges have not been observed or have been observed erroneously Guimerà and Sales-Pardo (2009); Peixoto (2018). The latent conditional probabilities of edges existing between groups of nodes is learned by the collective observation of many similar edges, and these correlations are leveraged to extrapolate the existence of missing or spurious ones. The same mechanism is expected to aid the reconstruction task, where edges are not observed directly, but the observed functional behavior yields a posterior distribution on them, allowing the same kind of correlations to be used as an additional source of evidence for the reconstruction, going beyond what the dynamics alone says.

Our reconstruction approach is finalized by defining an appropriate model for the functional behavior, determining $P(\bm{\mathcal{D}}|\bm{A})$. Here we will consider two kinds of indirect data. The first comes from a SIS epidemic spreading model Pastor-Satorras et al. (2015), where $\sigma_{i}(t)=1$ means node $i$ is infected at time $t$, [math] otherwise. The likelihood for this model is

[TABLE]

where

[TABLE]

is the transition probability for node $i$ at time $t$, with $f(p,\sigma)=(1-p)^{\sigma}p^{1-\sigma}$, and where $m_{i}(t)=\sum_{j}A_{ij}\ln(1-\tau_{ij})\sigma_{j}(t)$ is the contribution from all neighbors of node $i$ to its infection probability at time $t$. In the equations above the value $\tau_{ij}$ is the probability of an infection via an existing edge $(i,j)$, and $\gamma$ is the $1\to 0$ recovery probability. With these additional parameters, the full posterior distribution for the reconstruction becomes

[TABLE]

Since $\tau_{ij}\in[0,1]$ we use the uniform prior $P(\bm{\tau})=1$. Note also that the recovery probability $\gamma$ plays no role on the reconstruction algorithm, since its term in the likelihood does not involve $\bm{A}$ (and hence, gets cancelled out in the denominator $P(\bm{\sigma}|\gamma)=P(\gamma|\bm{\sigma})P(\bm{\sigma})/P(\gamma)$). This means that the above posterior only depends on the infection events $0\to 1$, and thus is also valid without any modifications to all epidemic variants SI, SIR, SEIR, etc Pastor-Satorras et al. (2015), since the infection events occur with the same probability for all these models.

The second functional model we consider is the Ising model, where spin variables on the nodes $\bm{s}\in\{-1,1\}^{N}$ are sampled according to the joint distribution

[TABLE]

where $\beta$ is the inverse temperature, $J_{ij}$ is the coupling on edge $(i,j)$, $h_{i}$ is a local field on node $i$, and $Z(\bm{A},\beta,\bm{J},\bm{h})=\sum_{\bm{s}}\exp(\beta\sum_{i<j}J_{ij}A_{ij}s_{i}s_{j}+\sum_{i}h_{i}s_{i})$ is the partition function. Note that this is not a dynamical model, as each microstate $\bm{s}$ is sampled independently according to the above distribution. Unlike the SIS model considered before, this distribution cannot be written in closed form since $Z(\bm{A},\beta,\bm{J},\bm{h})$ cannot be computed exactly, rendering the reconstruction problem intractable. Therefore, instead, we make use of the pseudolikelihood approximation Besag (1974), which is very accurate for the purpose at hand Nguyen et al. (2017), where we approximate Eq. 8 as a product of (properly normalized) conditional probabilities for each spin variable $s_{i}$

[TABLE]

With the above likelihood, reconstruction is performed by observing a set of $M$ microstates $\bar{\bm{s}}=\{\bm{s}_{1},\dots,\bm{s}_{M}\}$, sampled according to $P(\bar{\bm{s}}|\bm{A},\beta,\bm{J},\bm{h})=\prod_{l}P(\bm{s}_{l}|\bm{A},\beta,\bm{J},\bm{h})$, which yields the posterior distribution

[TABLE]

In the above we use uniform priors $P(\bm{J}|\bm{A})=\prod_{ij}[-1/2<J_{ij}<1/2]^{A_{ij}}$, thus forcing, without loss of generality, the values of $J_{ij}$ to lie in the shifted unit interval $[-1/2,1/2]$. For the remaining parameters we use uniform priors, $P(\bm{h})\propto 1$ and $P(\beta)\propto 1$, for $\beta\in[-\infty,\infty]$ and $\bm{h}\in[-\infty,\infty]^{N}$.

For any of the above posterior distributions, we perform sampling using Markov chain Monte Carlo (MCMC): For each proposal $\bm{A}\to\bm{A}^{\prime}$, it is accepted with the Metropolis-Hastings probability Metropolis et al. (1953); Hastings (1970)

[TABLE]

and likewise for the node partition $\bm{b}\to\bm{b}^{\prime}$, and any of the remaining parameters $\bm{\theta}\to\bm{\theta}^{\prime}$. Note that the acceptance probability does not require the intractable normalization constant $P(\bm{\mathcal{D}})$ to be computed. For both functional models considered, a whole sweep over $E$ entries of the adjacency matrix and $N$ nodes is done in time $O(EM+N{\left<k\right>})$, where $M$ is the number of data samples per node, allowing the method to be applied for large systems. We summarize and give more details about the technical aspects of the algorithm in Appendix C.

Synthetic networks — We begin by investigating the reconstruction performance of networks sampled from the planted partition model (PP), i.e. a DC-SBM with $\kappa_{i}=1$, $\lambda_{rs}=\lambda_{\text{in}}\delta_{rs}+\lambda_{\text{out}}(1-\delta_{rs})$, with $\lambda_{\text{in}}={\left<k\right>}(1+\epsilon(B-1))/N$ and $\lambda_{\text{out}}={\left<k\right>}(1-\epsilon)/N$, where $\epsilon=N(\lambda_{\text{in}}-\lambda_{\text{out}})/{\left<k\right>}B$ controls the strength of the modular structure. The detectability threshold for this model is given by $\epsilon^{*}=1/\sqrt{{\left<k\right>}}$, below which it is impossible to recover the planted community structure Decelle et al. (2011). Given a network $\bm{A}^{*}$ from this model, we sample $M$ independent Ising microstates $\bm{s}$ according to Eq. 8 using $J_{ij}=1$, $h_{i}=0$ and $\beta=\beta^{*}$ being the critical inverse temperature for the particular network. We compare two inference approaches: In the first we sample both the reconstructed network as well as its community structure form the joint posterior of Eq. 10. In the second approach, we perform reconstruction and community detection separately, by first performing reconstruction in isolation, by replacing the DC-SBM prior $P(\bm{A}|\bm{b})$ by the likelihood of an Erdős-Rényi model. We evaluate the quality of the reconstruction via the posterior similarity $S\in[0,1]$, defined as

[TABLE]

where $\bm{A}^{*}$ is the true network and $\bm{\pi}$ is the marginal posterior probability for each edge, i.e. $\pi_{ij}=\sum_{\bm{A},\bm{b},\bm{\theta}}A_{ij}P(\bm{A},\bm{b},\bm{\theta}|\bm{\mathcal{D}})$. A value $S=1$ means perfect reconstruction. We then perform community detection a posteriori by obtaining the maximum marginal point estimate

[TABLE]

and then sampling from the posterior $P(\bm{b}|\hat{\bm{A}})$. Fig. 1 contains the comparison between both approaches for networks sampled from the PP model, which shows how sampling from the joint posterior improves both the reconstruction as well as community detection. For the latter, the joint inference allows the detection all the way down to the detectability threshold, for the examples considered, which, otherwise, is not possible with the separate method.

Real networks with synthetic dynamics — Now, we investigate the reconstruction of networks not generated by the DC-SBM. We take two empirical networks, the worldwide network of $N=3\,286$ airports 111Retrieved from openflights.org. with $E=39\,430$ edges, and a food web from Little Rock Lake Martinez (1991), containing $N=183$ nodes and $E=2\,434$ edges, and we sample from the SIS (mimicking the spread of a pandemic) and Ising model (representing simplified inter-species interactions) on them, respectively, and evaluate the reconstruction obtained via the joint and separate inference with community detection, with results shown in Fig. 2. As is also the case for synthetic networks, the reconstruction quality is significantly improved by performing joint community detection 222Note that in this case our method also exploits the heterogeneous degrees in the network via the DC-SBM, which can in principle also aid the reconstruction, in addition to the community structure itself.. The quality of the reconstruction peaks at the critical threshold for each model, at which the sensitivity to perturbations is the largest. As the number of observed samples increases, so does the quality of the reconstruction, and the relative advantage of the joint reconstruction diminishes, as the data eventually “washes out” the contribution from the prior. For the Ising model, we compare the results of our method with the mean-field inverse correlations method Nguyen et al. (2017), i.e. $\beta A_{ij}J_{ij}=[\bm{C}^{-1}]_{ij}$, where $C_{ij}={\left<\sigma_{i}\sigma_{j}\right>}-{\left<\sigma_{i}\right>}{\left<\sigma_{j}\right>}$ is the connected correlation matrix. As seen in Fig. 2, this simpler reconstruction method can be just as accurate as our separate reconstruction approach, but only close to the critical point. For higher inverse temperatures the reconstruction deteriorates rapidly, and breaks down completely as the system becomes locally magnetized, with whole rows and columns of the matrix $\bm{C}$ being equal to zero, causing it to be singular 333Refinements of this approach including TAP and BP corrections Nguyen et al. (2017) yield the same performance for this example.. In such situations this kind of approach requires explicit regularization techniques Decelle and Ricci-Tersenghi (2014), which become unnecessary with our Bayesian method. The joint inference with community structure improves the reconstruction even further, beyond what is obtainable with typical inverse Ising methods, since it incorporates a different source of evidence.

In Fig. 3 we show a comparison of the reconstruction of the food web network from a simulated Ising model, using different approaches. Optimal thresholding corresponds to the naive approach of imputing the existence of an edge to the connected correlation between two nodes exceeding a threshold $c^{*}$, i.e. $\pi_{ij}=\{1\text{ if }C_{ij}>c^{*},\;0\text{ otherwise}\}$. The value of $c^{*}$ was chosen to maximize the posterior similarity, which represents the best possible reconstruction achievable with this method. Nevertheless, the network thus obtained is severely distorted. The inverse correlation method comes much closer to the true network, but is superseded by the joint inference with community detection.

Empirical dynamics — We turn to the reconstruction from observed empirical dynamics with unknown underlying interactions. The first example is the sequence of $M=619$ votes of $N=575$ deputies in the 2007 to 2011 session of the lower chamber of the Brazilian congress. Each deputy voted Yes, No, or abstained for each legislation, which we represent as $\{1,-1,0\}$, respectively. Since the temporal ordering of the voting sessions is likely to be of secondary importance to the voting outcomes, we assume the votes are sampled from an Ising model (the addition of zero-valued spins changes Eq. 9 only slightly by replacing $2\cosh(x)\to 1+2\cosh(x)$). Fig 4 shows the result of the reconstruction, where the division of the nodes uncovers a cohesive government and a split opposition, as well as a marginal center group, which correlates very well with the known party memberships and can be use to predict unseen voting behavior (see Appendix D). In Fig 5 we show the result of the reconstruction of the directed network of influence between $N=1\,833$ twitter users from $58\,224$ re-tweets Hodas and Lerman (2014) using a SI epidemic model (the act of “re-tweeting” is modelled as an infection event, using Eqs. 5 and 6 with $\gamma=0$) and the nested DC-SBM. The reconstruction uncovers isolated groups with varying propensities to re-tweet, as well as groups that tend to be influence a large fraction of users. By inspecting the geo-location metadata on the users, we see that the inferred groups amount to a large extent do different countries, although clear sub-divisions indicate that this is not the only factor governing the influence among users (see Appendix D.2).

Conclusion — We have presented a scalable Bayesian method to reconstruct networks from functional observations that uses the SBM as a structured prior, and, hence, performs community detection together with reconstruction. The method is nonparametric, and, hence, requires no prior stipulation of aspects of the network and size of the model, such as number of groups. By leveraging inferred correlations between edges, the SBM includes an additional source of evidence, and, thereby, improves the reconstruction accuracy, which in turn also increases the accuracy of the inferred communities. The overall approach is general, requiring only appropriate functional model specifications, and can be coupled with an open ended variety of such models, other than those considered here.

Appendix A Nonparametric DC-SBM model summary

The DC-SBM used in this work is the same derived in detail in Ref. Peixoto (2017b). We give a succinct summary in the following. The marginal likelihood of the DC-SBM can be written as

[TABLE]

where $\bm{e}=\{e_{rs}\}$ is the matrix of edge counts betwen groups, and $\bm{k}$ is the degree sequence of the network, and with

[TABLE]

being the microcanonical likelihood and corresponding noninformative priors. We further increase the explanatory power of this model Peixoto (2017b) by replacing the microcanonical prior for the degrees with

[TABLE]

where $\bm{\eta}=\{\eta_{k}^{r}\}$ are the degree frequencies of each group, with $\eta_{k}^{r}$ being the number of nodes with degree $k$ that belong to group $r$, and

[TABLE]

is a uniform distribution of degree sequences constrained by the overall degree counts, and finally

[TABLE]

is the distribution of the overall degree counts. The quantity $q(m,n)$ is the number of different degree counts with the sum of degrees being exactly $m$ and that have at most $n$ non-zero counts, given by

[TABLE]

For the node partition we use the prior,

[TABLE]

which is agnostic to group sizes.

Finally, the hierarchical degree-corrected SBM (HDC-SBM) is obtained by replacing the uniform prior for $P(\bm{e}|\bm{b})$ by a nested sequence of SBMs, where the edge counts in level $l$ are generated by a SBM at a level above,

[TABLE]

where $\left(\kern-3.00003pt\left(\genfrac{}{}{0.0pt}{}{n}{m}\right)\kern-3.00003pt\right)={n+m-1\choose m}$ is the multiset coefficient. The prior for the hierarchical partition is obtained using Eq. 22 at every level. The entire model above is also easily modified for directed networks. We refer to Ref. Peixoto (2017b) for further details.

Appendix B Adapting multigraph models to simple graphs

The DC-SBM variations considered above generate multigraphs with self-loops, however the functional models presented in the main text operate on simple graphs. We amend this inconsistency in the same manner as in Ref. Peixoto (2018), by adapting the multigraph models to simple graphs in tractable way by generating multigraphs and then collapsing the multiple edges. In other words, if $\bm{G}$ is a multigraph with entries $G_{ij}\in\mathbb{N}$, the collapsed simple graph $\bm{A}(\bm{G})$ has binary entries

[TABLE]

Therefore, if $\bm{G}$ is a multigraph generated by $P(\bm{G}|\theta)$, where $\theta$ are arbitrary parameters, then the corresponding collapsed simple graph $\bm{A}$ is generated by

[TABLE]

with

[TABLE]

Even if $P(\bm{A}|\theta)$ cannot be computed in closed form, the joint distribution $P(\bm{A},\bm{G}|\theta)=P(\bm{A}|\bm{G})P(\bm{G}|\theta)$ is trivial, provided we have $P(\bm{G}|\theta)$ in closed form. Therefore, instead of directly sampling from the posterior distribution

[TABLE]

we sample from the joint posterior

[TABLE]

using MCMC, treating the values $G_{ij}$ as latent variables, and then we marginalize

[TABLE]

which is done simply by sampling from $P(\bm{A},\bm{G},\bm{b}|\bm{\mathcal{D}})$ and ignoring the actual magnitudes of the $G_{ij}$ values, and the diagonal entries.

Appendix C Inference algorithm

The inference algorithm used here is identical to Ref. Peixoto (2018), with the only difference being the likelihoods for the forward model $P(\bm{\mathcal{D}}|\bm{A})$. To summarize, we use MCMC to sample from the joint posterior distribution

[TABLE]

where $\bm{b}$ is the partition of nodes used for the SBM. The MCMC algorithm consists of making proposals of the kind $P(\bm{b}^{\prime}|\bm{A},\bm{b})$ and $P(\bm{A}^{\prime}|\bm{A},\bm{b})$ for the partition and network, respectively (or equivalently for any other remaining model parameter), and accepting them according to the Metropolis-Hastings probability

[TABLE]

which does not require the computation of the intractable normalization constant $P(\bm{\mathcal{D}})$. In practice, at each step in the chain we make either a move proposal for $\bm{A}$ or $\bm{b}$, not both at once. For the node partition, we use the move proposals described in Refs. Peixoto (2014a, 2017b), where for any given node $i$ in group $r$ we propose to move it to group $s$ (which can be previously unoccupied, in which case it is labelled $s=B+1$) according to

[TABLE]

where $P(t|i)=\sum_{j}A_{ij}\delta_{b_{j},t}/k_{i}$ is the fraction of neighbours of $i$ that belong to group $t$, $\epsilon>0$ is a small parameter which guarantees ergodicity, and $d$ is the probability of moving to a previously unoccupied group. (If $k_{i}=0$, we assume $P(b_{i}=r\to s|\bm{A},\bm{b})=d\delta_{s,B+1}+(1-d)(1-\delta_{s,B+1})/B$.) This move proposal attempts to the use the currently known large-scale structure of the network to better inform the possible moves of the node, without biasing with respect to group assortativity. The parameters $d$ and $\epsilon$ do not affect the correctness of the algorithm, only the mixing time, which is typically not very sensitive, provided they are chosen within a reasonable range (we used $d=0.01$ and $\epsilon=1$ throughout). When using the HDC-SBM, we used the variation of the above for hierarchical partitions described in Ref. Peixoto (2017b). The move proposals above require only minimal bookkeeping of the number edges incident on each group, and can be made in time $O(k_{i})$, which is also the time required to compute the ratio in Eq. 32, independent on how many groups are currently occupied.

For the network, we change the values of the latent multigraph $\bm{G}$ with unit proposals

[TABLE]

for $\delta\in\{-1,1\}$. We choose the entries to update with a probability given by the current DC-SBM,

[TABLE]

with

[TABLE]

being the probability of selecting node $i$ from its group $b_{j}$, proportional to its current degree plus one, and

[TABLE]

is the probability of selecting groups $(r,s)$, where $e_{rs}=\sum_{ij}G_{ij}\delta_{b_{i},r}\delta_{b_{j},s}$. The above probabilities guarantee that every entry will be eventually sampled, but it tends to probe denser regions more frequently, which we found to typically lead to faster mixing times. This sampling can be done in time $O(1)$, simply by keeping urns of vertices and edges according to the group memberships. The time required to compute the ratio in Eq. 32 is also $O(1)$ for the DC-SBM and $O(L)$ for the HDC-SBM, where $L$ is the hierarchy depth, again independent of the number of occupied groups.

C.1 Algorithmic complexity

When combining the move proposals defined above for the partition and network, the time required to perform $N$ node move proposals and $E$ edge addition or removal proposals is $O({\left<k\right>}N+EM)$, where ${\left<k\right>}$ is the average degree, and $M$ is the number of samples per node of the functional model (i.e. SIS or Ising). The $O(EM)$ contribution is seen by noting that the addition and removal of an edge requires the re-computation of the likelihood $P(\bm{\mathcal{D}}|\bm{A})$ involving only terms associated with each endpoint over all $M$ samples, each requiring only $O(1)$ computations. For the SIS model we note that we need only to keep track of the summary quantities $m_{i}(t)$ for each node, and update them by adding or subtracting contributions for each added or removed edge, and the same is true for the Ising model with respect to edge contributions to the Hamiltonian. This linear complexity of sweeps allows for the reconstruction of large networks.

For dynamical data where changes of the state of each node are relatively rare (e.g. in a SI or SIR dynamics, a node changes its state only once or twice, respectively, for the whole cascade), it is possible to optimize the inference algorithm by listing for each node only its initial state and the times it changes, together with the new state values. In this way, the contribution to the likelihood of a single node can be computed by only going through the times that its neighbours or the node itself change state, instead of the whole time-span of the dynamics. Therefore, the complexity for a whole MCMC sweep changes to $O({\left<k\right>}N+Ea)$, where $a$ is the average number of times a single node changes its state during the whole dynamics. For very active dynamics we have $a=O(M)$, and hence this algorithm has the same complexity as the version above, but for $a\ll O(M)$ it gives noticeable speed-ups.

In addition to the algorithmic complexity of each sweep, the MCMC needs time to converge to the target posterior distribution. This mixing time depends not only on the structure of the network being reconstructed, how easy it is to uncover it from the data, but also on how close the algorithm is initiated to the target distribution. Because of this, it not straightforward to estimate the general algorithmic complexity of the mixing time. In our numerical experiments we found that both starting from random or empty networks lead to reasonable mixing times in most cases, and the results coincide with initializing from the true planted network (which as expected, shows faster equilibration).

A reference implementation of the above algorithm is freely available as part of the graph-tool library Peixoto (2014c).

Appendix D Datasets with empirical dynamics

D.1 Cross-validation

To further evaluate the ability of our reconstruction method to capture the empirical voting behavior of deputies in the lower house of the Brazilian congress, we randomly divided all $M=619$ voting sessions into $M-M_{t}$ “training” sessions, used to fit the model, and $M_{t}$ test sessions, used to compare with the predictions from the model. To evaluate the predicition error, the correlation matrix $C_{ij}={\left<s_{i}s_{j}\right>}-{\left<s_{i}\right>}{\left<s_{j}\right>}$ was computed from the posterior distribution of the fitted model via

[TABLE]

where $\bar{\bm{s}}_{t}$ is the training data, and compared with the correlation matrix obtained for the test data $C_{ij}^{(t)}={\left<s_{i}s_{j}\right>}_{t}-{\left<s_{i}\right>}_{t}{\left<s_{j}\right>}_{t}$, via

[TABLE]

where the sums go over microstates of the test data. The prediction error is then computed as

[TABLE]

We repeated the calculation using $M_{t}\in\{10,40,100\}$, and for each value of $M_{t}$ we averaged the results over $30$ random choices of the test data. We also compared with the reconstruction obtained via inverse correlations Nguyen et al. (2017). The results are shown in Fig. 6. As can be seen, the Bayesian joint reconstruction method outperforms the results based on inverse correlations. The prediction error decreases with larger $M_{t}$, as in this limit the fluctuations in the dynamics are averaged out.

D.2 Comparison with metadata

Here we expand on the comparison of the community structure found both for the Brazilian congress as well as the twitter data, with metadata available in both cases.

In table 2 we list the party affiliations of the deputies according to the groups they were classified by the method. The largest group accounts for all left-wing parties as well as the center parties belonging to the government collation, whereas groups 2 and 3 accounts for the right-wing opposition. Group 4 is composed of a small number deputies who are members of both government and oppositions parties, but vote independently.

In table 1 is shown the country of each twitter user, independently obtained via twitter’s API (not contained in the original dataset of Ref. Hodas and Lerman (2014)), according to each group identified by the reconstruction method. As can be seen, most groups are characterized by a single dominating country, with only a few exceptions. This indicates, plausibly, that the probability of re-tweets is largely shaped by language and cultural barriers. Nevertheless, the method also uncovers subdivisions within the distinct geographical locations, indicating that this is not the only factor determining the influence among users.