Thresholding normally distributed data creates complex networks

TL;DR

This paper investigates how thresholding normally distributed, correlated data can produce network properties typical of complex networks, revealing which features are emergent and which are not.

Contribution

It introduces a simple correlated data model and analyzes the network properties resulting from thresholding, highlighting which complex network features are emergent.

Findings

Heavy-tailed degree distributions observed

Large numbers of triangles present

Short path lengths observed

Abstract

Network data sets are often constructed by some kind of thresholding procedure. The resulting networks frequently possess properties such as heavy-tailed degree distributions, clustering, large connected components and short average shortest path lengths. These properties are considered typical of complex networks and appear in many contexts, prompting consideration of their universality. Here we introduce a simple model for correlated relational data and study the network ensemble obtained by thresholding it. We find that some, but not all, of the properties associated with complex networks can be seen after thresholding the correlated data, even though the underlying data are not "complex". In particular, we observe heavy-tailed degree distributions, a large numbers of triangles, and short path lengths, while we do not observe non-vanishing clustering or community structure.

Equations

$$X_{ij}\ge t$$

$$\text{Var}[X_{ij}]$$

$$\text{Cov}[X_{ij},X_{ik}]$$

$$\text{Cov}[X_{ij},X_{kl}]$$

$$P\left(\mathbf{x}\right)=\frac{e^{-\frac{1}{2}(\mathbf{x}-\mathbf{\mu }{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{x}-\mathbf{\mu })}}{\sqrt{(2\pi {)}^N\text{det}(\mathbf{\Sigma })}}.$$

$$\mathbf{X}\sim \mathcal{N}(0,\mathbf{\Sigma }),$$

$${\Sigma }_{(i,j),(i,j)}$$

$${\Sigma }_{(i,j),(i,k)}$$

$${\Sigma }_{(i,j),(k,l)}$$

$$A_{ij}=\{\begin{array}{ll}1& \text{ if }{X}_{ij}\ge t,\\ 0& \text{ otherwise}.\end{array}$$

$$W_{ij}=\sqrt{1-2\rho }{Y}_{ij}+\sqrt{\rho }\left(Z_i+Z_j\right).$$

$$\text{Var}\left[W_{ij}\right]$$

$$\text{Cov}\left[W_{ij},W_{ik}\right]$$

$$\text{Cov}\left[W_{ij},W_{kl}\right]$$

$$y>\frac{t-\sqrt{\rho }\left(z_i+z_j\right)}{\sqrt{1-2\rho }}$$

$$E[A_{ij}]=P[A_{ij}=1]=P[X_{ij}\ge t]=1-\Phi (t),$$

$$E[k_i]=\sum _{j\ne i}E[A_{ij}]=(n-1)(1-\Phi (t)).$$

$$t={\Phi }^{-1}\left(1-\frac{⟨k⟩}{n-1}\right).$$

$$C=P[A_{ik}=1|A_{ij},A_{jk}=1]=\frac{P[A_{ik},A_{ij},A_{jk}=1]}{P[A_{ij},A_{jk}=1]}.$$

$$P[X_{ij},X_{ik}\ge t]=\sum _{N=0}^{\infty }\frac{{\rho }^N}{N!}{\left[\varphi (t)H_{N-1}\left(t\right)\right]}^2$$

$$P[X_{ij},X_{ik},X_{jk}\ge t]=$$

$$\sum _{N=0}^{\infty }\sum _{i=0}^N\sum _{j=0}^{N-i}\frac{{\rho }^N\varphi {\left(t\right)}^3}{i!j!(N-i-j)!}{H}_{N-1-i}\left(t\right)H_{N-1-j}\left(t\right)H_{i+j-1}\left(t\right)$$

$$T=\left(\genfrac{}{}{0ex}{}{n-1}{2}\right)P[X_{ij},X_{ik},X_{jk}\ge t].$$

$$\frac{1}{2}\left(⟨k^2⟩-⟨k⟩\right)=\left(\genfrac{}{}{0ex}{}{n-1}{2}\right)P[X_{ij},X_{ik}\ge t].$$

$$\text{Var}[k]$$

$$+(n-1)(n-2)\sum _{N=1}^{\infty }\frac{{\rho }^N}{N!}{\left[\varphi (t)H_{N-1}\left(t\right)\right]}^2.$$

$$p_k=\left(\genfrac{}{}{0ex}{}{n-1}{k}\right)\sqrt{\frac{1-\rho }{2\pi \rho }}{\int }_{-\infty }^{\infty }{e}^{f_k\left(y\right)}dy$$

$$f_k(y)=k$$

$$-\frac{1}{2}{\left(\frac{t-\sqrt{1-\rho }y}{\sqrt{\rho }}\right)}^2.$$

$$y_{0,k}={\Phi }^{-1}\left(1-\frac{k}{n-1}\right).$$

$$p_k\sim \frac{1}{n-1}\sqrt{\frac{1-\rho }{\rho }}\exp \left[-\left(\frac{1-2\rho }{2\rho }\right)y_{0,k}^2+\left(\frac{t\sqrt{1-\rho }}{\rho }\right)y_{0,k}-\frac{t^2}{2\rho }\right].$$

$$P[X_{ij},X_{ik}\ge t]=\frac{1}{2\pi \sqrt{1-{\rho }^2}}{\int }_t^{\infty }{\int }_t^{\infty }{e}^{-\frac{1}{2}\left(\frac{x^2-2\rho xy+y^2}{1-{\rho }^2}\right)}dxdy.$$

$$H_n(x)=(-1{)}^n{e}^{\frac{x^2}{2}}\frac{d^n}{d{x}^n}{e}^{-\frac{x^2}{2}}.$$

$$H_{-1}(x)=\frac{1-\Phi \left(x\right)}{\varphi \left(x\right)}.$$

$${\Sigma }_{jk}^{(0)}={\Sigma }_{(0,j),(0,k)}=\{\begin{array}{ll}1& \mathrm{for}j=k,\\ \rho & \mathrm{otherwise}.\end{array}$$

Full text

I Introduction

Networks are a popular tool for representing and analyzing real-world systems consisting of entities and their relationships. They provide a simple yet intuitive representation for many complex systems. In the most basic incarnation, networks are simple graphs—undirected and unweighted with only one type of node and one type of edge. The usual picture is that nodes represent some group of objects (people, neurons, proteins, etc.), and edges represent some kind of interaction between them (friendship, synapses, binding, etc.) Newman18 ; berg_structure_2004 ; bullmore_brain_2011 ; stopczynski_measuring_2014 ; apeltsin_improving_2011 ; bullmore_complex_2009 ; rubinov_complex_2010 ; sekara_strength_2014 .

In many real-world settings interactions are indicated by real-valued data, and so creating a simple network requires thresholding, which may take several forms apeltsin_improving_2011 ; bullmore_complex_2009 ; langer_problem_2013 ; rubinov_complex_2010 ; huang2009network ; sekara_strength_2014 ; stopczynski_measuring_2014 ; serrano2009extracting ; radicchi2011information ; dianati2016unwinding . The most obvious case of thresholding is when a continuous valued data set is explicitly thresholded by deciding what level of interaction is sufficiently strong to count as an edge in the network. A more subtle case is that of experimental limitation: interactions that exist but are very weak or rare may not be observed. Even for binary valued data sets the sampling method may hide an implicit thresholding mechanism. For example, one commonly uses a combination of a yeast-two-hybrid screen and biochemical assays to detect and verify edges in protein-protein interaction networks. These methods typically do not detect weak protein-protein interactions VAYNBERG200622 and are thus equivalent to applying a threshold on the edge strength in protein-protein interaction networks. For another example, consider friendship networks. Most everyday interactions between people are presumably not strong enough to constitute friendship. At what point does a casual acquaintance cross over to the category of friend? When people list their friends, in a survey for instance, they will implicitly apply some criteria to filter the friends from the acquaintances. Nevertheless, an understanding of the properties one should expect to observe from thresholded relational data is currently lacking.

In this paper we examine the properties of networks created by thresholding relational data. To do this we introduce a basic model of the underlying relational data which is then thresholded to produce edges in the network. The model is derived from three assumptions:

all nodes are statistically identical; 2. 2.

any correlations are local; 3. 3.

the underlying relational data are normally distributed.

All three of these assumptions—which are no doubt violated in real-world systems—are quite natural for a null model. Assumption 1, that all nodes are identical, severely constrains what correlation structures are admissible. In fact, only two free parameters remain in the covariance matrix once this assumption is made: a local correlation strength between edges that share nodes, and a global correlation strength between edges that do not share nodes. Assumption 2 sets the second of these to zero—edges that do not have a node in common are uncorrelated. The other free parameter, the local correlation strength, we call $\rho$. Our remaining freedom is to pick a distribution that is consistent with the required correlation matrix. The most obvious and simple choice is assumption 3, the multivariate normal (Gaussian) distribution. We believe this to be the simplest non-trivial model for relational data.

The thresholding procedure will also be very simple: any of the relational data that falls above some threshold, $t$, will be said to constitute an edge in the network, and any that falls below will not. The threshold value $t$ is a parameter of the model.

Sophisticated methods to extract networks from weighted data have been developed, for example in serrano2009extracting ; radicchi2011information ; dianati2016unwinding . These more complicated methods lead to different networks, but in this paper we do not consider the relative merits of more advanced procedures. We favor the simplistic approach since it contains only one parameter and allows us to derive equations for several network properties. Nevertheless, in Appendix E we present some similar results for the so-called disparity filter of serrano2009extracting .

Our network ensemble on $n$ nodes is thus defined by two parameters: the threshold, $t$, and a local correlation coefficient, $\rho$. Despite the simplicity of the model—the underlying relational data are normally distributed—we nonetheless find a number of the behaviors typically observed in complex networks, such as heavy-tailed degree distributions, short average path lengths, and large numbers of triangles. It does not, however, yield non-vanishing clustering or community structure in the large $n$ limit and so cannot account for this observation in real-world data sets. Finally, the model we study is not constrained to produce positive definite matrices. As a result it is not immediately applicable to the study of thresholded correlation matrices (as for example in huang2009network ).

This paper has two main parts. In Sec. II we define and justify the network model. Then, in Sec. III, we study the properties of the network ensemble. We look at the density of edges, triangles and clustering, the degree distributions, shortest path lengths, and the giant component.

II Model specification

II.1 Thresholding locally correlated data

A network can be represented by its adjacency matrix, $\bm{A}$, where $A_{ij}=1$ if node $i$ and $j$ are connected and $A_{ij}=0$ otherwise. We consider networks created by thresholding underlying relational data, $\bm{X}$, adding an edge between $i$ and $j$ if

[TABLE]

(see Fig. 1a). To fully specify the model we need to pick a distribution for $\bm{X}$. Assuming that all nodes are statistically identical—exchangeable in the parlance of statistics—constrains our choice of distribution.

If nodes are identical then the marginal distribution for $X_{ij}$ must be the same for all (distinct) pairs $i$ and $j$. Further, by a linear transform we can always set $E[X_{ij}]=0$ and $\text{Var}[X_{ij}]=1$. So long as the appropriate transformation is made to $t$, this shift will have no effect on the thresholded network. For this reason we will always assume $X_{ij}$ has mean [math] and variance $1$. Exchangeability puts further constraints on the covariance matrix, whose entries can take only three values. For $i,j,k,l$ all distinct, these are

[TABLE]

where $\mathrm{Cov}[X,Y]$ denotes covariance. We will assume that $\gamma=0$ since this quantifies the correlation between two edges that do not share a node, i.e. two edges that do not “touch” (see Fig. 1b). This leaves us with two free parameters, $t$ and $\rho$. The remaining task is to pick a distribution with the required covariance matrix, $\bm{\Sigma}$.

In principle any distribution could be used, but the obvious choice is a multivariate normal distribution. In standard notation a multivariate normal distribution (MVN) is denoted $\mathcal{N}(\bm{\mu},\bm{\Sigma})$. The probability density function of an $N$-dimensional MVN is

[TABLE]

The normal distribution has many points in its favor. Famously it arises in the central limit theorem, which makes it a plausible model for many random processes. If the relational data $\bm{X}$ arises due to the aggregation of many independent processes then the central limit theorem implies $\bm{X}$ will be multivariate normally distributed. Further, the normal distribution is the maximum entropy distribution with the required covariance matrix, Eq. (2), and so could be justified as the “least informative distribution”—the model that makes the fewest extra assumptions beyond the correlation structure. We can also appeal to simple pragmatism: the multivariate normal distribution is well-studied and has convenient mathematical properties.

A concise statement of the model is as follows: given the freely chosen parameters $t\in\mathbb{R},\rho\in{\left[0,\frac{1}{2}\right]}$, and the number of nodes $n$, draw a random variable $\bm{X}$ with

[TABLE]

where

[TABLE]

Then create the network by thresholding $\bm{X}$,

[TABLE]

Note, we constrain $0\leq\rho\leq\frac{1}{2}$ so that $\bm{\Sigma}$ is positive semi-definite 111To see why $\rho>\frac{1}{2}$ is problematic, consider the marginal distribution for four edges, say $X_{ij},X_{jk},X_{kl},X_{il}$. A simple calculation shows that the covariance matrix has a negative eigenvalue for $\rho>\frac{1}{2}$. Similarly, since $\text{Var}\left[\Sigma_{j}X_{ij}\right]$ must be greater than zero, $\rho$ must be greater than $-1/(n-2)$ and so negative correlations can be vanishingly weak at most..

Even if we have good reason to believe that the marginal distributions for $X_{ij}$ are not normal, the model may still be applicable. Consider an arbitrary cumulative distribution function, $F(x)$, and let $\Phi(x)$ denote the standard normal cumulative distribution function. If we sample $\bm{X}$ from a multivariate normal distribution, and then apply the function $F^{-1}(\Phi(x))$ to each $X_{ij}$ we will have transformed the edge weights to the arbitrary distribution $F$. So long as we apply the same transformation to $t$, however, the resulting network after thresholding will be identical.

The upshot is that our model can be adapted for any marginal distribution, and no network properties change—the assumption that the edge weights have normally distributed marginals is of no real consequence. What is important is the assumption that there is some transformation of the data such that the joint distribution is multivariate normal. While this assumption is a limitation, the above procedure is actually one of the standard methods for creating multivariate distributions with arbitrary marginals.

II.2 Sampling from the model

We now describe a simple algorithm to sample from the model. This algorithm also provides an intuitive model interpretation.

Let $Z_{i}$ be $n$ i.i.d. variables, $\mathcal{N}(0,1)$. Let $Y_{ij}$ be ${n\choose 2}$ i.i.d. variables, $\mathcal{N}(0,1)$. Then let

[TABLE]

Note that $W_{ij}$ is normally distributed with mean zero and further

[TABLE]

Hence, $\bm{W}$ is distributed identically to $\bm{X}$. So, to sample from the model:

Sample $\bm{z}$, a length n vector of i.i.d. standard normal variables. 2. 2.

For $i<j$, generate $y\sim\mathcal{N}(0,1)$, and if

[TABLE]

add edge $(i,j)$ to the network.

If $\rho=\frac{1}{2}$, generating $y$ is unnecessary and one can simply add edge $(i,j)$ if $\sqrt{1/2}(z_{i}+z_{j})\geq t$.

A Python package to generate networks along with scripts for the figures in this paper is publicly available maier_thredgecorr_2018 .

In order to achieve the required correlations, the algorithm above separates $X_{ij}$ into node and edge effects. Each node is given a value $Z_{i}$ and $X_{ij}$ is created by a linear combination of $Z_{i}$ and $Z_{j}$ plus i.i.d. random noise $Y_{ij}$. We can interpret the $Z$’s as latent variables that control the propensity for individual nodes to have edges and $\rho$ controls the relative strength of the noise process. When $\rho=1/2$ edges are entirely determined by the values of $Z$, while at $\rho=0$ edges are entirely random and independent.

Despite this equivalent formulation, our model should not be primarily understood as a latent variable model since it was not constructed as one. Rather, the equivalent latent variable model is derived and used for algorithmic convenience. In fact, the existence of this latent variable interpretation is not surprising. As $n\to\infty$ our model is in a class of models known as exchangeable random graphs diaconis2007graph ; orbanz2015bayesian . The Aldous-Hoover theorem implies that all exchangeable random graphs have an equivalent latent variable model diaconis2007graph ; orbanz2015bayesian ; hoover1979relations .

III Network properties

We now turn our attention to properties of the networks created by the model.

III.1 Edge density

Edges in the network exist whenever the corresponding weight $X_{ij}$ is greater than $t$. The marginal distribution for $X_{ij}$ is simply a standard normal distribution. Thus,

[TABLE]

where $\Phi(x)$ is the cumulative distribution function for the standard normal distribution $\mathcal{N}(0,1)$. When $\rho=0$ all edges exist independently and the model is equivalent to the random graph, $G_{n,p}$, with $p=1-\Phi(t)$.

The mean degree is equally simple to compute. For all $\rho$

[TABLE]

If we want to pick $t$ for a desired mean degree $\left\langle k\right\rangle$, it is easy to invert this to obtain

[TABLE]

III.2 Triangles, clustering, and degree variance

Many complex networks are observed to have large numbers of triangles. The clustering coefficient or transitivity is one way to quantify this. We can quantify the clustering with the probability that a triangle is closed, given that two of its edges already exist,

[TABLE]

The numerator of this equation corresponds to the density of triangles while the denominator corresponds to the density of two-stars (which also determines the variance of the degree distribution). Note that for simplicity we shorten the logical connective “and” (or “$\wedge$”) using commas, e.g. $P[A_{ij}=1\ \wedge\ A_{jk}=1]\equiv P[A_{ij},A_{jk}=1]$.

The marginal distributions of a MVN are themselves MVN, and are found by simply dropping the unwanted rows and columns in the correlation matrix $\bm{\Sigma}$. Thus, $(X_{ij},X_{ik})^{T}$ will be bivariate normally distributed and $(X_{ij},X_{ik},X_{jk})^{T}$ will be trivariate normally distributed, both with correlation coefficient $\rho$. Introducing the Hermite polynomials $H_{N}(x)$ as defined in Appendix A, one finds that

[TABLE]

for the density of two-stars and

[TABLE]

for triangles. Both sums converge for $\rho\leq 0.5$, and we can estimate them accurately with a finite number of terms harris_use_1980 . Noting that there are ${n-1\choose 2}$ potential triangles for each node, the expected number of triangles per node is simply ${n-1\choose 2}$ times their density

[TABLE]

Plots of these functions are shown in Fig. 2. We find that $T$ is much larger in these networks than in the random graph $G_{n,p}$—larger by multiple orders of magnitude. In fact, while $T$ goes to zero in the large $n$ limit for the random graph, in this model we find that $T$ increases with $n$ for large values of $\rho$. On the other hand, the clustering coefficient $C$ decreases with growing number of nodes for all parameter values. This leads to a slightly paradoxical result for large $\rho$: in the limit $n\rightarrow\infty$ the expected number of triangles at each node goes to infinity, and the clustering coefficient still goes to zero! The reason for this is that the number of two-stars diverges faster than the number of triangles.

Equation (14) can also be used to compute the variance of the degree distribution. To see this note that a node of degree $k$ has ${k\choose 2}$ two-stars. Further, noting that there are ${n-1\choose 2}$ potential two-stars (the same number of potential triangles) we find

[TABLE]

Combining this with Eq. (11) the variance of the node degree $k$ can be written

[TABLE]

The first term is simply the variance of a binomial distribution. For $\rho=0$ the second term vanishes and we recover the correct result for the random graph $G_{n,p}$. For $\rho>0$ the sum is positive and monotonically increases with $\rho$ as illustrated in Fig. 3.

III.3 Degree distribution

In the previous two subsections we gave expressions for the mean and variance of the degrees. Here we give expressions for the full distribution of degrees.

The degree distribution $p_{k}$ is the probability that a node has $k$ edges. For this model the degree distribution can be written

[TABLE]

where

[TABLE]

This result is derived in Appendix B.

The integral in Eq. (19) can be computed numerically to high precision using Gauss-Hermite quadrature, centered at the maximum of $f_{k}(y)$. Increasing the order of Gauss-Hermite quadrature (i.e. incorporating more points) increases the accuracy. The full details are in Appendix B.

We can also approximate the integral using Laplace’s method miller2006applied , an asymptotic approximation for integrals of this form (equivalent to a first order Gauss-Hermite quadrature). The idea of the method is to replace the function $f_{k}(y)$ by a second order Taylor series around its maximum. For large $n$, the last (quadratic) term in $f_{k}$ will be negligible and for $0<k<n-1$, the maximum will be at

[TABLE]

Combining this with Stirling’s approximation for the binomial coefficient, we find

[TABLE]

Together with the closed form approximation for $\Phi^{-1}$, given in Appendix C, Eq. (22) provides a closed form approximation for the degree distribution.

Figure 4 shows some example degree distributions, computed to high precision using Eq. (19) along with the asymptotic approximation, Eq. (22), where we chose $n=100\,000$ and $\left\langle k\right\rangle=100$.

To illustrate how these degree distributions compare to the degree distributions of real networks, we chose three data sets from different domains, and fit the model. The first data set is a network of friendships between students at a U.S. high school ($n=2587$) AddHealth , the second data set is a co-authorship network of researchers ($n=16726$) newman_structure_2001 , and the third network describes interactions between proteins ($n=6327$) joshi-tope_reactome:_2005 .

Given a number of nodes $n$ the model studied in this paper has two free parameters, $t$ and $\rho$. A simple procedure to fit the model to the data is to choose $t$ and $\rho$ so that the mean and variance of the model’s degree distribution match the observed values. We use Eq. (12) to fix $t$ and subsequently Newton’s method to solve Eq. (18) for $\rho$.

The results of this exercise are shown in Fig. 5. The networks were chosen for their different degree distributions—note the different scales on the axes: linear, log-linear, and log-log—and the threshold model can qualitatively ape these distributions. Nevertheless, the similarity of the degree distribution should not be over-emphasized. As discussed, this model has vanishing clustering so cannot account for this observation of real-world networks. Comparisons for clustering are shown in Appendix D.

While the degrees in the thresholded networks, $k_{i}=\sum_{j}A_{ij}$, in general follow a complicated distribution, the underlying degrees $d_{i}=\sum_{j}X_{ij}$ are always normally distributed. When $\rho=0$, $d_{i}$ is Gaussian and $k_{i}$ is binomial, or Poisson in the sparse limit. When $\rho>0$, $d_{i}$ is still Gaussian, but $k_{i}$ now follows a heavy-tailed distribution. Thus, the heavy-tailed distribution observed in the model is due to the combination of correlation and thresholding. Without positive correlation we observe Poisson distributions; without thresholding we observe Gaussian distributions.

III.4 Giant component

A well studied problem in the theory of random graphs is the formation of a large connected (giant) component. At very low densities only a handful of nodes can be reached from any other node but at some critical point a macroscopic number of nodes will be connected. For the random graph this transition occurs at a mean degree of $\left\langle k\right\rangle=1$ Newman18 ; ER60 ; bollobas_2001 .

To explore the effects of $\rho>0$ we sampled from the model as described in Sec. II.2 and measured the size of the second largest component as a susceptibility parameter for the phase transition. The maximum of this susceptibility parameter is used to find the transition lines in Fig. 6a.

We find that as $\rho$ or $n$ increases, the transition occurs at lower values of the mean degree. This result is in line with the configuration model for which the transition point decreases with increasing variance in the degree distribution. For $\rho=0$ we recover the standard result for the random graph.

For the other limit case, $\rho=1/2$, recall that all edge weights can be considered to arise from node “propensities”, $Z_{i}$, with $X_{ij}=\sqrt{1/2}(Z_{i}+Z_{j})$. This implies that all nodes that are connected to any other nodes must also be connected to the node with maximum propensity $Z_{\max}$. The size of the largest component is then given by this node’s degree plus $1$, $k_{\max}+1$. The second largest component is then always of size 1. We therefore omit $\rho=1/2$ in the numerical analysis.

III.5 Shortest path lengths

Another phenomenon well established in the complex networks literature is that randomly chosen nodes often have surprisingly short paths between them. This is often referred to as the “six degrees of separation” or “small-world” phenomenon watts_collective_1998 ; Newman18 . By a common definition, network models are considered to demonstrate this property if the average shortest path length $\left\langle d_{ij}\right\rangle$ between nodes grows logarithmically (or slower) as the number of nodes increases Newman18 .

Using the method described in Sec. II.2, we sampled from the threshold model to verify that it displays this property. We looked at networks with between $100$ and $30\,000$ nodes, with mean degree $\left\langle k\right\rangle=5$, and investigated the influence of increasing edge weight correlation $\rho$. After sampling a network from the model we computed the average shortest path length $\left\langle d_{ij}\right\rangle$ on the largest (giant) component. For each parameter combination we computed the mean by averaging 200 sampled networks.

The results are shown in Fig. 7a. Since it is well known that the random graph $G_{n,p}$ has short shortest paths FR07 it is unsurprising that the threshold model does also (recall, for $\rho=0$ they are equivalent, and we see the standard $\left\langle d_{ij}\right\rangle\propto\log n$ scaling behavior). For $\rho>0$ we see that average shortest path lengths grow significantly slower than logarithmically, a behavior sometimes referred to as “ultra small-world” and often related to networks with power-law degree distribution cohen_scale-free_2003 ; bornholdt_structural_2004 . In our model, the effect appears despite the fact that the degree distribution does not follow a power-law.

As discussed, when $\rho=1/2$ all edge weights can be considered to arise from node propensities $Z_{i}$, such that $X_{ij}=\sqrt{1/2}(Z_{i}+Z_{j})$. All nodes are then either disconnected or part of the giant component, and the node with maximum propensity $Z_{\max}$ is connected to all nodes in the giant component. Hence, all nodes in the giant component are either directly connected or can reach each other in two steps through the maximum-degree node. So, when $\rho=1/2$ the average shortest path length must be $1\leq\left\langle d_{ij}\right\rangle<2$.

IV Discussion

In this paper we studied the effects of thresholding relational data. We started with a simple model of multivariate normally distributed data, with only one free parameter, $\rho$, controlling local correlations. We then demonstrated that thresholding this normally distributed correlated relational data reproduces many of the properties commonly associated with complex networks. In particular, we find that the combined effects of correlation and thresholding leads to heavy-tailed degree distributions, relatively large numbers of triangles, and short average path lengths.

The underlying relational data $\bm{X}$ in the model we introduce would not usually be considered complex. It is generated from a highly symmetric multivariate normal distribution with only one free parameter. Since every pair of nodes has some level of interaction, the graphical interpretation for $\bm{X}$ would be a weighted complete graph, with all edge weights (and linear combinations thereof) normally distributed. For example, the “degrees”, $d_{i}=\sum_{j}X_{ij}$, are normally distributed. And yet, after thresholding the networks show several properties commonly associated with complex networks.

One way to think about these results is in the context of the central limit theorem. Whenever interaction strengths are the aggregate result of a large number of processes then we expect $\bm{X}$ to be normally distributed. Constructing a simple graph from these data can lead to complex networks. This provides one simple explanation for the ubiquity of complex networks—they can arise as a consequence of the central limit theorem.

Of course, for most scientific questions of interest the exact details of the mechanisms and structure are what matter. In a social network, for example, answering the question "who influences whom, and why?" is far from trivial and the fact that the network has certain commonly observed properties is usually incidental.

In summary, straightforward assumptions lead to several of the properties associated with complex networks. If a network arises by a simple thresholding procedure then finding that it is “complex” need be no more surprising than finding a bell-shaped curve in a regular data set.

Acknowledgements.

The authors thank the students, staff, and faculty of the Santa Fe Institute Complex Systems Summer School for bringing them together and providing a stimulating research environment. All authors made significant contributions to the project. G.T.C and B.F.M. wrote the final manuscript. B.F.M. was supported as an Add-On Fellow for Interdisciplinary Life Science by the Joachim Herz Stiftung. A.C.S. was supported by the Clarendon Fund, e-Therapeutics plc, and the Engineering and Physical Sciences Research Council under grant number EP/L016044/1. J.S. was supported by the U.S. Army Research Office under Multidisciplinary University Research Initiative Award No. W911NF-13-1-0340. G.S. was supported by the Natural Sciences and Engineering Research Council of Canada and the Sentinel North program.

Appendix A Multivariate normal integrals and Hermite polynomials

The probability of a two-star existing with nodes $i,j$ and $k$ as constituents is given by

[TABLE]

Direct computation of the integral is not straightforward but we can compute it quickly using the Hermite polynomials harris_use_1980 . A quick outline of this method: for $n\geq 0$, define the Hermite polynomials as

[TABLE]

As the name suggests, the Hermite polynomials are in fact polynomials, for example $H_{0}(x)=1$, $H_{1}(x)=x$, $H_{2}(x)=x^{2}-1$, and so on. For notational convenience also define

[TABLE]

Using the Hermite polynomials, we can expand Eq. (23) as an infinite sum and integrate term by term. The final result is given by Eq. (14). The same trick is used for the 3-dimensional integral to give Eq. (15).

Appendix B Degree distribution

Since, by assumption, all nodes in this model are equivalent, we will simply consider the one-node marginal to compute the degree distribution. Let $\bm{U}$ be all the terms in $\bm{X}$ that are associated with node [math], i.e. $U_{j}=X_{0j}$. Then, $\bm{U}$ is multivariate normally distributed, $\mathcal{N}(\bm{0},\bm{\Sigma}^{(0)})$, where $\bm{\Sigma}^{(0)}$ has ones along the diagonal and $\rho$ everywhere else

[TABLE]

The focal node will have degree $k$ when exactly $k$ terms in $\bm{U}$ are larger than the threshold $t$. There are ${n-1\choose k}$ different ways this can happen and each is equally likely. So, to compute $p_{k}$ we can compute the probability that the first $k$ terms in $\bm{U}$ are larger than $t$ and all others are smaller, and then multiply by ${n-1\choose k}$ to obtain

[TABLE]

To solve this integral we use a standard trick tong_multivariate_1990 . First, we note that if $Z_{0},Z_{1},\dots,Z_{n-1}$ are i.i.d. $\mathcal{N}(0,1)$ then

[TABLE]

will be distributed identically to $\bm{U}$. Further, once we know the value of $Z_{0}$ then all the terms are independent, and the probability that any one of them is greater than $t$ is the probability that $Z_{1}\geq\frac{t-\sqrt{\rho}z}{\sqrt{1-\rho}}$. Given $Z_{0}=z$, the probability that exactly $k$ values will greater than $t$ and the rest less than $t$ is

[TABLE]

Averaging this quantity over $z$ then provides us with the correct expression,

[TABLE]

where $I_{n,k}$ is the integral. A change of variables allows us to write

[TABLE]

where

[TABLE]

A standard approach to approximate such an integral is to use Laplace’s method. In this approach one expands $f$ about its maximum and then neglects higher order terms, $f(y)\approx f(y_{0})-\frac{|f^{\prime\prime}(y_{0})|}{2}(y-y_{0})^{2}$. Having done this, the integral reduces to a standard Gaussian integral. While this approach is asymptotically correct (in the large $n$ and $k$ limit), we can improve the approximation by including more terms using Gauss-Hermite quadrature. Re-writing the integral again, and making another change of variables:

[TABLE]

where $R_{k}$ is the remaining terms of $f_{k}$ after expansion:

[TABLE]

Now we can approximate the integral using Gauss-Hermite quadrature:

[TABLE]

where $x_{i}$ are the points for which $H_{N}\left(x_{i}\right)=0$ and the weights $w_{i}$ are

[TABLE]

Note that $I_{n,k}(1)$ is Laplace’s approximation, i.e. Laplace’s approximation is a first order Gauss-Hermite quadrature at the maximum of $f_{k}$, while $I_{n,k}(N)$ approximates the remainder terms with increasingly high order polynomials and so we expect $I_{n,k}(N)\rightarrow I_{n,k}$ as $N$ increases.

Appendix C Approximation of inverse cumulative distribution function

The normal distribution’s inverse cumulative distribution function, $\Phi^{-1}\left(x\right)$, can be approximated abramowitz_handbook_2013 for $0<x\leq 0.5$ as

[TABLE]

with

[TABLE]

For $0.5<x\leq 1$ we use $\Phi^{-1}(x)=-\Phi^{-1}(1-x)$.

Appendix D Other properties of real networks

In Fig. 8 we compare simulations from the thresholded normal model to the real networks from Sec. III.3. Other degree properties, such as the average neighbor-degree seem to be modelled well, while the local clustering coefficient is generally smaller in the simulations than real data, as expected.

Appendix E Disparity filtering

The simplistic style of thresholding that we have considered is not the only method to extract a network from relational data. Indeed, more sophisticated algorithms have been developed serrano2009extracting ; radicchi2011information ; dianati2016unwinding . By limiting our analysis to a simple thresholding procedure, rather than a more sophisticated algorithm, we have been able to derive several of the basic properties. In this appendix we repeat some of the analysis for the more complex “disparity filter” algorithm serrano2009extracting .

The disparity filter assumes that the relational data are positive and vary by orders of magnitude. In contrast, we have considered data that are distributed according to a standard normal distribution—$X_{ij}$ is negative half of the time and will virtually never be larger in magnitude than $10$. To match the assumptions of the disparity filter, we apply the algorithm to $e^{X_{ij}}$.

Since $X_{ij}$ is normally distributed, $e^{X_{ij}}$ is log-normally distributed. And, since the exponential function is monotonic, a simple thresholding procedure at the value of $e^{t}$ is equivalent to our previous analyses. However, since $e^{X_{ij}}$ is non-negative and will vary by orders of magnitude, these data also match the assumptions of the disparity filter.

The disparity filter algorithm proceeds by assigning a local significance score $\alpha_{ij}\neq\alpha_{ji}$ to each potential edge. In the normal model this score will be

[TABLE]

An edge is considered present between nodes $i$ and $j$ if either $\alpha_{ij}$ or $\alpha_{ji}$ exceeds a predetermined significance threshold, $\alpha$.

Making use of the decomposition from Eq. (7), we can write

[TABLE]

Further, we can make a mean-field style approximation and replace $\sum_{k\neq j}e^{\sqrt{1-2\rho}Y_{ik}+\sqrt{\rho}Z_{k}}$ by its expected value to arrive at

[TABLE]

Defining the constant

[TABLE]

a quick calculation establishes that either $\alpha_{ij}$ or $\alpha_{ji}$ is larger than $\alpha$ with probability

[TABLE]

where $\phi$ and $\Phi$ are again the density function and cumulative distribution function for the standard normal distribution. This integral can easily be computed using Gauss-Hermite quadrature. Its derivatives with respect to $\alpha$ are also simple to compute, and so root finding algorithms such as Newton’s method can find the correct choice of $\alpha$ for a desired final density.

While the above derivation for the network density in this model is reasonably straight-forward, other properties such as the degree distribution are more challenging. Instead of re-deriving the results of this paper for a sophisticated filtering algorithm we present the results of simulations in Fig. 9. Increasing the correlation $\rho$ similarly increases the variance on the node degree $k$. This demonstrates that the qualitative properties derived for the naive thresholding still apply for at least some more sophisticated thresholding methods.