Hierarchical structure and the prediction of missing links in networks

References

Supplementary Information

Appendix A Hierarchical random graphs

Our model for the hierarchical organization of a network is as follows.Computer code implementing many of the analysis methods described in this paper can be found online at www.santafe.edu/$\sim$aaronc/randomgraphs/. Let $G$ be a graph with $n$ vertices. A dendrogram $D$ is a binary tree with $n$ leaves corresponding to the vertices of $G$. Each of the $n-1$ internal nodes of $D$ corresponds to the group of vertices that are descended from it. We associate a probability $p_{r}$ with each internal node $r$. Then, given two vertices $i,j$ of $G$, the probability $p_{ij}$ that they are connected by an edge is $p_{ij}=p_{r}$ where $r$ is their lowest common ancestor in $D$. The combination $(D,\{p_{r}\})$ of the dendrogram and the set of probabilities then defines a hierarchical random graph.

Note that if a community has, say, three subcommunities, with an equal probability $p$ of connections between them, we can represent this in our model by first splitting one of these subcommunities off, and then splitting the other two. The two internal nodes corresponding to these splits would be given the same probabilities $p_{r}=p$. This yields three possible binary dendrograms, which are all considered equally likely.

We can think of the hierarchical random graph as a variation on the classical Erdős–Rényi random graph $G(n,p)$. As in that model, the presence or absence of an edge between any pair of vertices is independent of the presence or absence of any other edge. However, whereas in $G(n,p)$ every pair of vertices has the same probability $p$ of being connected, in the hierarchical random graph the probabilities are inhomogeneous, with the inhomogeneities controlled by the topological structure of the dendrogram $D$ and the parameters $\{p_{r}\}$. Many other models with inhomogeneous edge probabilities have, of course, been studied in the past. One example is a structured random graph in which there are a finite number of types of vertices with a matrix $p_{kl}$ giving the connection probabilities between them.F. McSherry, “Spectral Partitioning of Random Graphs.” Proc. Foundations of Computer Science (FOCS), pp. 529–537 (2001)

Appendix B Fitting the hierarchical random graph to data

Now we turn to the question of finding the hierarchical random graph or graphs that best fits the observed real-world network $G$. Assuming that all hierarchical random graphs are a priori equally likely, the probability that a given model $(D,\{p_{r}\})$ is the correct explanation of the data is, by Bayes’ theorem, proportional to the posterior probability or likelihood $\mathcal{L}$ with which that model generates the observed network.G. Casella and R. L. Berger, “Statistical Inference.” Duxbury Press, Belmont (2001). Our goal is to maximize $\mathcal{L}$ or, more generally, to sample the space of all models with probability proportional to $\mathcal{L}$.

Let $E_{r}$ be the number of edges in $G$ whose endpoints have $r$ as their lowest common ancestor in $D$, and let $L_{r}$ and $R_{r}$, respectively, be the numbers of leaves in the left and right subtrees rooted at $r$. Then the likelihood of the hierarchical random graph is

If we fix the dendrogram $D$, it is easy to find the probabilities $\{\overline{p}_{r}\}$ that maximize $\mathcal{L}(D,\{p_{r}\})$. For each $r$, they are given by

the fraction of potential edges between the two subtrees of $r$ that actually appear in the graph $G$. The likelihood of the dendrogram evaluated at this maximum is then

Figure S1 shows an illustrative example, consisting of a network with six vertices.

It is often convenient to work with the logarithm of the likelihood,

where $h(p)=-p\log p-(1-p)\log(1-p)$ is the Gibbs-Shannon entropy function. Note that each term $-L_{r}R_{r}h(\overline{p}_{r})$ is maximized when $\overline{p}_{r}$ is close to or to $1$, i.e., when the entropy is minimized. In other words, high-likelihood dendrograms are those that partition the vertices into groups between which connections are either very common or very rare.

We now use a Markov chain Monte Carlo method to sample dendrograms $D$ with probability proportional to their likelihood $\mathcal{L}(D)$. To create the Markov chain we need to pick a set of transitions between possible dendrograms. The transitions we use consist of rearrangements of subtrees of the dendrogram as follows. First, note that each internal node $r$ of a dendrogram $D$ is associated with three subtrees: the subtrees $s,t$ descended from its two daughters, and the subtree $u$ descended from its sibling. As Figure S2 shows, there are two ways we can reorder these subtrees without disturbing any of their internal relationships. Each step of our Markov chain consists first of choosing an internal node $r$ uniformly at random (other than the root) and then choosing uniformly at random between the two alternate configurations of the subtrees associated with that node and adopting that configuration. The result is a new dendrogram $D^{\prime}$. It is straightforward to show that transitions of this type are ergodic, i.e., that any pair of finite dendrograms can be connected by a finite series of such transitions.

Once we have generated our new dendrogram $D^{\prime}$ we accept or reject that dendrogram according to the standard Metropolis–Hastings rule.M. E. J. Newman and G. T. Barkema, “Monte Carlo Methods in Statistical Physics.” Clarendon Press, Oxford (1999). Specifically, we accept the transition $D\to D^{\prime}$ if $\Delta\log\mathcal{L}=\log\mathcal{L}(D^{\prime})-\log\mathcal{L}(D)$ is nonnegative, so that $D^{\prime}$ is at least as likely as $D$; otherwise we accept the transition with probability $\exp(\log\Delta\mathcal{L})=\mathcal{L}(D^{\prime})/\mathcal{L}(D)$. If the transition is not accepted, the dendrogram remains the same on this step of the chain. The Metropolis-Hastings rule ensures detailed balance and, in combination with the ergodicity of the transitions, guarantees a limiting probability distribution over dendrograms that is proportional to the likelihood, $P(D)\propto\mathcal{L}(D)$. The quantity $\Delta\log\mathcal{L}$ can be calculated easily, since the only terms in Eq. (4) that change from $D$ to $D^{\prime}$ are those involving the subtrees $s$, $t$, and $u$ associated with the chosen node.

The Markov chain appears to converge relatively quickly, with the likelihood reaching a plateau after roughly $O(n^{2})$ steps. This is not a rigorous performance guarantee, however, and indeed there are mathematical results for similar Markov chains that suggest that equilibration could take exponential time in the worst case.E. Mossel and E. Vigoda, “Phylogenetic MCMC Are Misleading on Mixtures of Trees.” Science 309, 2207 (2005) Still, as our results here show, the method seems to work quite well in practice. The algorithm is able to handle networks with up to a few thousand vertices in a reasonable amount of computer time.

We find that there are typically many dendrograms with roughly equal likelihoods, which reinforces our contention that it is important to sample the distribution of dendrograms rather than merely focusing on the most likely one.

Appendix C Resampling from the hierarchical random graph

The procedure for resampling from the hierarchical random graph is as follows.

Initialize the Markov chain by choosing a random starting dendrogram.

Run the Monte Carlo algorithm until equilibrium is reached.

Sample dendrograms at regular intervals thereafter from those generated by the Markov chain.

For each sampled dendrogram $D$, create a resampled graph $G^{\prime}$ with $n$ vertices by placing an edge between each of the $n(n-1)/2$ vertex pairs $(i,j)$ with independent probability $p_{ij}=\overline{p}_{r}$, where $r$ is the lowest common ancestor of $i$ and $j$ in $D$ and $\overline{p}_{r}$ is given by Eq. (2). (In principle, there is nothing to prevent us from generating many resampled graphs from a dendrogram, but in the calculations described in this paper we generate only one from each dendrogram.)

After generating many samples in this way, we can compute averages of network statistics such as the degree distribution, the clustering coefficient, the vertex-vertex distance distribution, and so forth. Thus, in a way similar to Bayesian model averaging,T. Hastie, R. Tibshirani and J. Friedman, “The Elements of Statistical Learning.” Springer, New York (2001). we can estimate the distribution of network statistics defined by the equilibrium ensemble of dendrograms.

For the construction of consensus dendrograms such as the one shown in Fig. 2a, we found it useful to weight the most likely dendrograms more heavily, giving them weight proportional to the square of their likelihood, in order to extract a coherent consensus structure from the equilibrium set of models.

Appendix D Predicting missing connections

Our algorithm for using hierarchical random graphs to predict missing connections is as follows.

Initialize the Markov chain by choosing a random starting dendrogram.

Run the Monte Carlo algorithm until equilibrium is reached.

Sample dendrograms at regular intervals thereafter from those generated by the Markov chain.

For each pair of vertices $i,j$ for which there is not already a known connection, calculate the mean probability $\langle p_{ij}\rangle$ that they are connected by averaging over the corresponding probabilities $p_{ij}$ in each of the sampled dendrograms $D$.

Sort these pairs $i,j$ in decreasing order of $\langle p_{ij}\rangle$ and predict that the highest-ranked ones have missing connections.

In general, we find that the top 1% of such predictions are highly accurate. However, for large networks, even the top 1% can be an unreasonably large number of candidates to check experimentally. In many contexts, researchers may want to consider using the procedure interactively, i.e., predicting a small number of missing connections, checking them experimentally, adding the results to the network, and running the algorithm again to predict additional connections.

The alternative prediction methods we compared against, which were previously investigated inD. Liben-Nowell and J. Kleinberg, “The link prediction problem for social networks.” Proc. Internat. Conf. on Info. and Know. Manage. (2003)., consist of giving each pair $i,j$ of vertices a score, sorting pairs in decreasing order of their score, and predicting that those with the highest scores are the most likely to be connected. Several different types of scores were investigated, defined as follows, where $\Gamma(j)$ is the set of vertices connected to $j$.

Common neighbors: score$(i,j)=|\Gamma(i)\,\cap\,\Gamma(j)|$, the number of common neighbors of vertices $i$ and $j$.

Jaccard coefficient: score$(i,j)=|\Gamma(i)\,\cap\,\Gamma(j)|\,/\,|\Gamma(i)\,\cup\,\Gamma(j)|$, the fraction of all neighbors of $i$ and $j$ that are neighbors of both.

Degree product: score$(i,j)=|\Gamma(i)|\,|\Gamma(j)|$, the product of the degrees of $i$ and $j$.

Short paths: score$(i,j)$ is $1$ divided by the length of the shortest path through the network from $i$ to $j$ (or zero for vertex pairs that are not connected by any path).

One way to quantify the success of a prediction method, used by previous authors who have studied link prediction problems7, is the ratio between the probability that the top-ranked pair is connected and the probability that a randomly chosen pair of vertices, which do not have an observed connection between them, are connected. Figure S4 shows the average value of this ratio as a function of the percentage of the network shown to the algorithm, for each of our three networks. Even when fully $50\%$ of the network is missing, our method predicts missing connections about ten times better than chance for all three networks. In practical terms, this means that the amount of work required of the experimenter to discover a new connection is reduced by a factor of $10$, an enormous improvement by any standard. If a greater fraction of the network is known, the accuracy becomes even greater, rising as high as $200$ times better than chance when only a few connections are missing.

We note, however, that using this ratio to judge prediction algorithms has an important disadvantage. Some missing connections are much easier to predict than others: for instance, if a network has a heavy-tailed degree distribution and we remove a randomly chosen subset of the edges, the chances are excellent that two high-degree vertices will have a missing connection and such a connection can be easily predicted by simple heuristics such as those discussed above. The AUC statistic used in the text, by contrast, looks at an algorithm’s overall ability to rank all the missing connections over nonexistent ones, not just those that are easiest to predict.

Finally, we have investigated the performance of each of the prediction algorithms on purely random (i.e., Erdős–Rényi) graphs. As expected, no method performs better than chance in this case, since the connections are completely independent random events and there is no structure to discover. We also tested each algorithm on a graph with a power-law degree distribution generated according to the configuration model.M. Molloy and B. Reed, “A critical point for random graphs with a given degree sequence.” Random Structures and Algorithms 6, 161–179 (1995) In this case, guessing that high-degree vertices are likely to be connected performs quite well, whereas the method based on the hierarchical random graph performs poorly since these graphs have no hierarchical structure to discover.