Resolution of ranking hierarchies in directed networks

Abstract

Identifying hierarchies and rankings of nodes in directed graphs is fundamental in many applications such as social network analysis, biology, economics, and finance. A recently proposed method identifies the hierarchy by finding the ordered partition of nodes which minimises a score function, termed agony. This function penalises the links violating the hierarchy in a way depending on the strength of the violation. To investigate the resolution of ranking hierarchies we introduce an ensemble of random graphs, the Ranked Stochastic Block Model. We find that agony may fail to identify hierarchies when the structure is not strong enough and the size of the classes is small with respect to the whole network. We analytically characterise the resolution threshold and we show that an iterated version of agony can partly overcome this resolution limit.

About the Authors:

Elisa Letizia

Roles Conceptualization, Formal analysis, Methodology, Software, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Scuola Normale Superiore, Piazza dei Cavalieri 7, Pisa, 56126, Italy

ORCID http://orcid.org/0000-0002-3134-1137

Paolo Barucca

Roles Conceptualization, Methodology

Affiliations University of Zurich, Schönberggasse 1, Zürich, 8001 Switzerland, LIMS, 35a South Street, London, W1K 2XF, United Kingdom

Fabrizio Lillo

Roles Conceptualization, Methodology, Writing – original draft, Writing – review & editing

Affiliation: Department of Mathematics, University of Bologna, Piazza di Porta San Donato 5, Bologna, 40126, Italy

Introduction

Identifying ranking hierarchies in complex networks is of paramount importance in many disciplines and applications. An exact hierarchical organisation in a directed network means that the set of nodes can be divided in an ordered collection of classes such that links exist only from a node of a low rank class to a node of a higher rank class. Clearly it is equivalent to define exact hierarchical structures when links exist only from upper to lower classes.

Recently the relevance of measuring ranking hierarchy was pointed out in the context of ecosystems [1], in which it was shown how species exhibit a property of trophic coherence, measuring how consistently a species falls into a distinct level of hierarchy within a food web. Other major applications include social network analysis [2, 3], the study of funds flow in financial networks [4, 5], and of corporate cross-ownerships in economics [6].

Since real networks are not necessarily exactly hierarchical, the problem considered here is to find an optimal ordered partition of nodes into classes such that the structure has a maximal level of hierarchy.

Framed in this way, the procedure is to choose a suitable hierarchy metric and to devise feasible algorithms which find the ordered partition of nodes in such a way that the hierarchy metric is maximised. In [7], this maximisation was recognised to be a dual problem of circulation, known to be related to the cost max-flow minimisation [8]. The problem is analogous to the more explored problem of community detection in graphs [9]. In such case a common approach is to choose a metric, for example the modularity, and to look for partitions that maximizes it. It is well known that modularity has resolution limits [10, 11], and the associated optimisation problem might be a hard computational task, even if successful heuristics exist [12].

It is important to stress that the concept of ranking hierarchy we employ in this paper, introduced in [13], and further developed in [7, 14–17], models graphs, representing for example social organisations, as command structure or influential communities.

Related literature sharing a similar definition of hierarchy includes [18–22]. This concept is therefore very different from the more common definition of nested hierarchy in networks, studied for example in [23–26], where low-level communities of nodes are nested into bigger ones, in a way directly associated with hierarchical clustering. The former concept of hierarchy is defined in directed networks and look for rankings of nodes into classes, while the latter makes sense also for undirected networks and look for nested clusters of nodes. Still different, but closer it is the concept of ranking used in [27, 28] and related literature, where the ranking to be inferred is a complete order on a undirected network, not an ordered partition.

In this paper we consider the problem of the inference of hierarchies in directed networks via a class of metrics recently introduced and termed agony. Given a ranking of nodes into classes (i.e. an ordered partition), agony is a metric which penalises those links which are against the ranking, i.e. from a high rank to a low rank node. Different forms of penalisation lead to different types of agony. Once the agony function is chosen, one looks for the ranking of nodes which minimises it. Thus optimisation of agony is a non-parametric approach of hierarchy detection. Similarly to community detection with modularity, agony minimisation might be a challenging computational task, even if for some forms of the agony function exact or heuristic algorithms have been recently proposed (see the next Section for more details).

Here we focus our attention on the problem of resolution limiti when detecting ranking hierarchies with agony minimisation. Specifically, we ask when a given hierarchical structure in a network can be identified with agonies. To investigate the possible resolution limits we proceed in a way similar to what has been done for community detection [29–31]. We introduce a class of random graphs, termed Ranked Stochastic Block Models (RSBM) which is a specific subclass of Stochastic Block Models [32] with a tunable hierarchical structure and we study the resolution limit of hierarchy detection with agony in RSBM. Following the comment above, we stress again that our RSBM are different from those recently proposed in [26], where the nested concept of hierarchy is adopted.

We find, both analytically and empirically, that agony succeeds in identifying hierarchies when the structure is strong. However we prove the existence of resolution thresholds in the model parameters such that beyond these thresholds agony minimisation identifies hierarchical structures which are different from the planted one. Using symmetry arguments we explore analytically alternative rankings, showing that they can have a smaller agony (higher hierarchy) than the planted one. These rankings are obtained by merging, splitting, or inverting the classes of the planted ranking. It is important to notice that, as we show numerically, the resolution limits are not due to the RSBM, but to the method. Once more, this is analogous to what observed in community detection with modularity maximisation. Finally we show that in some cases, iterating the optimisation on each class found in the first run of the agony minimisation, it is possible to improve significantly the recovery of the planted structure.

Methods

Agony

Let G = (V, E) be a binary directed graph of N ≡ |V| nodes and m ≡ |E| links. A rank function r: V → {1, …, R} associates each node to an integer number which indicates the position of the subset (or class) containing the node in the hierarchy. Thus a rank function generates an ordered partition of the nodes into R subsets (i ∈ {1, …, R}) of size . From this point, we will refer to the ordered partition induced by the rank function with the term ranking. Once a ranking has been assigned to the graph G, a link between two nodes is classified as forward if it goes from a node in a class to one in a class with a strictly higher rank and backward otherwise. Identifying the optimal hierarchical structure in a directed graph means to find a ranking where the presence of backward links is suitably penalised. The penalisation will in general depend on the number of backward links as well as on the distance in rank between the connected nodes. The penalisation is of course arbitrary and it is interesting to investigate the ability of different forms of penalisation in identifying hierarchies.

The concept of agony in graphs was first introduced in [16] and it is the weighted cost of all the backward links in a ranking. More specifically, given a graph G and a ranking r the value of agony with respect to r is given by:(1)where f is a penalty function such that it is zero for negative argument and non decreasing otherwise. We will consider here f of the formand we will denote the value of agony of the ranking r on graph G with Ad(G, r). The agony of the graph is defined as the minimum value of agony with respect to all possible rankings on the nodes, i.e.(2)where denotes the set of all rankings. Fig 1 shows two examples of optimal rankings for simple graphs and illustrates the difference between backward and forward links.

[Figure omitted. See PDF.]

Fig 1. Optimal rank and agony (d = 1) for simple graphs.

On the top, the graphs are represented without any ordering, on the bottom nodes are divided according to their ranks. The red links are the backward ones, those contributing to agony, and the black links are the forward ones.

https://doi.org/10.1371/journal.pone.0191604.g001

Remarks 1

1. When the graph is a Directed Acyclical Graph (DAG), one can always find a ranking of the nodes such that there are no backward links (see [33] for a simple routine to solve this problem), hence the value of agony of a DAG is 0, and we say the graph has a perfect hierarchy.

2. The presence of +1 in the cost function f makes sure that same same class links are penalised. Without this, the optimal partition would always be the trivial one, with all the nodes in the same class.

3. Thanks to the minimisation, for the value of (generalised) agony for the optimal partition it holds . For the trivial partition rT, i.e. the one with all the nodes in the same class, it holds which gives the indicated upper bound for the optimum.

4. The exponent d acts as a tuning parameter: when it increases, only rankings with stronger hierarchies are privileged over the trivial one.

5. The optimal ranking may be not unique, however there exists a routine to choose the ranking with the smallest number of classes among those with the optimal value of agony (see [34] for more details).

Finally, one can define the hierarchy of a directed graph as(3)From the previous remark (ii) it follows that where indicates a perfect hierarchy.

Once the penalisation has been chosen, the problem of finding the optimal ranking is quite complex. In its original version, agony was defined with the piecewise linear cost function, i.e d = 1 in our notation. With this choice few exact algorithms to identify the optimal ranking of a graph are known [16, 34]. Ref. [34] considered the computational complexity of algorithms for generic d. The case d = 1 is proven to be solved by an algorithm of polynomial complexity, while the case d = 0 can be reformulated into the minimum Feedback Arc Set problem (FAS, or equivalently into the dual problem: Maximum Acyclic Subgraph) [35] which is known to be NP-hard, but for which some heuristics exist [36]. The intermediate cases, 0 < d < 1, have concave cost functions, which also lead to a NP-hard problem according to [34]. The case d > 1, instead, have a convex cost function which gives a problem of polynomial complexity. However, to the best of our knowledge, no algorithm is available at the moment for these latter cases. One of the objectives of this paper is to investigate how the detected optimal ranking depends on the choice of the penalty function. For this reason we need to introduce a class of graphs which have a hierarchical structure and whose strength can be tuned by a suitable choice of parameters. This is what we do in the next Subsection.

Ranked stochastic block model

Our ensemble of graphs belongs to the class of Stochastic Block Models (SBMs). In this ensemble of graphs, nodes are partitioned into R disjointed subsets and the probability of having a link between two nodes depends on the classes they belong to and it is independent of all the other pairs of nodes, i.e.The R × R matrix C = {cij}ij is called the affinity matrix. For our purpose we consider the directed version of SBMs, and C is not symmetric. We choose a parametrisation of C in order to keep the number of parameters small, which allows to have both analytical tractability, and enough flexibility to model different types of hierarchies.

The ranking r(p), which we will refer to as planted ranking, is defined so that it is consistent with the labelling in the affinity matrix, i.e.

This notion of planted ranking is slightly different from the one used in another stream of literature, as for example in [27]. In the latter case, it refers to the latent fully connected graph (clique) underlying the observed one (which is instead sparse). The former graph should be estimated in order to recover the latent ranking. This can indeed be related to the definition of the Stochastic Block Model: in this case the latent fully connected directed and weighted graph is obtained from the affinity matrix, by using as weights for each edge the probability to observe that specific link in one realisation of the ensemble, i.e the entry of the affinity matrix for the corresponding blocks. Moreover, the term ranking is used here as synonym of ordered partition, which implies that in general the order is not complete.

Note that, given the collection of subsets of nodes, any rank function with a range of values larger than R − 1 would have a larger value of Ad.

Consider

This gives the affinity matrixthat correspond to the schematic representation in Fig 2.

[Figure omitted. See PDF.]

Fig 2. Schematic representation of the structure of a RSBM graph.

https://doi.org/10.1371/journal.pone.0191604.g002

In order to have a true hierarchical structure we require that the parameters p, q, s are such that(4)Define ∀k ∈ 1, ‥, R,For any pair (i, j) (i, j = 1, ….R), the number of links between subset i and j mi, j follows a binomial distribution, mi,j ∼ Binom(ni nj, (C)i,j), therefore the constraint (4) is equivalent toIn the case of uniform cardinality of the subsets, ni = n∀i, which we will consider in the following, the inequality further simplifies to(5)A SBM having the above structure and satisfying the constraint (4) will be termed Ranked Stochastic Block Models RSBM(p, q, s, R, {ni}). In the case of uniform cardinality, we denote briefly RSBM(p, q, s, R).

Remarks 2

1. A possible interpretation of the RSBM is that p, q give the backbone of the hierarchical structure, while s represents the noise.

2. As mentioned in Remarks 1.2, the definition of the cost function implies that links between nodes of the same rank have a positive cost. This means that those links are classified as backward, and for this reason they are assigned a probability s as the other backward links in the model.

Since RSBMs are random graphs, different realisations of the model give different values of agony and hierarchy. We will compute below the expected value of these quantities.

We estimate the expected value of hd(G, r(p)), the hierarchy of the planted ranking of RSBM graphs. Note that we make a little abuse of notation indicating with hd the value 1 − Ad/m, i.e. we do not consider the minimisation of agony. For this reason h is not necessarily bounded between 0 and 1 as .

Indicating with the ensemble average of hd(G, r(p)), we obtainwhere In order to have closed form expressions we need to estimate the terms eij. We consider a second order Taylor expansion:(6)

If we assume that ni = O(N)∀i, then the last two terms in Eq (6) vanish when N → ∞, henceThis gives the first order estimate for It is possible to compute higher order estimates or estimates based on exact expected values. The expressions are however less transparent and we find in simulations that first order estimates are quite accurate, thus in the following we use them.

Results

Looking for optimal hierarchies in RSBM

RSBMs are constructed with a specific ranking, the planted one, which is determined by the choice of the classes and the model parameters. When minimising a generalised agony Ad on realisations of such graphs, it is not a priori obvious that the optimal ranking is the planted one. We therefore ask the following question:

Given a RSBM(p, q, s, R, {ni}), find the ranking r which minimises the generalised agony Ad. In particular check when the planted ranking r(p) is optimal.

This is in general a complicated problem and we do not have a complete answer to this question, despite the fact that it is possible, at least for d = 1, to find numerically the optimal ranking of a specific realisation of a RSBM. In order to simplify the problem, in this paper we will restrict our attention to the homogeneous case ni = N/R, ∀i. Given the form of the affinity matrix and the homogeneity assumption, we expect that the optimal solution, when different from the planted one, preserves the homogeneity of the planted ranking. Possible boundary effects (for example the first and last class have different size from the other ones) are not considered and we expect to play a minor role when the number of planted classes is large. In any case in the subsection Simulation we use numerical simulations to test our intuition.

For this reason we shall compute the generalised agony of the following alternative rankings:

1. the number of classes changes either by merging adjacent classes or by splitting each class; due to homogeneity, merged or split classes have all the same size;

2. the rank is inverted, , ∀j = 1, ‥, R, i.e. nodes in highest ranks of the planted ranking are given lowest ranks in the alternative. Moreover we consider also the case when the number of classes is arbitrary, but again their size is assumed to be uniform.

To distinguish between the two families of ranking, we will denote the former as direct, in contrast with inverted for the latter. For each of these alternative rankings we compute the value of as a function of the number of classes and we look for the optimal one among these alternatives and the planted ranking. Clearly there is no guarantee that this will be the global optimum over all the possible rankings. To maintain this distinction, we will denote optimal the ranking with highest value of within the subset of alternatives just described, while we will always refer to the best among all the rankings, i.e. that which gives , as the global optimum. We will see for example that numerical simulations of some RSBM indicate that the globally optimal ranking is a partial inversion of the planted hierarchy. However this analysis serves to show that planted ranking might not be globally optimal for some generalised agony and to provide an upper bound for the resolution threshold as well as getting intuition on the characteristics of the optimal ranking in a RSBM.

In the following we will focus on two regions of the parameter space of RSBMs:

* p ≥ q > s, termed a twitter-like hierarchy;

* q = 0, p ≠ 0 termed a military-like hierarchy.

In the former hierarchy forward links can connect low rank nodes with nodes of any higher rank, while in the latter the forward links can connect a node only with nodes in a direct superior class. In both cases backward links can exists with a probability s. As we will see the global optimal ranking of the two hierarchies is quite different.

Finally we consider the casewhere is the number of classes after splitting (b < 0) or merging (b > 0). The parameters a > 1 and b < a are such that . We denote the direct and inverted rankings with 2a−b classes as r(b) and r(i,b) respectively.

We will focus our attention on the case d = 1, d = 0, and d = 2. Results for other values of d are left for a future paper.

Agony with d = 1.

In this case exact algorithms for its optimisation are known, allowing the comparison of calculations with numerical simulations.

Provided that the constraints in (5) are satisfied, one can easily verify that ∀b < 0i.e. splitting is never optimal, neither in the direct nor in the inverted ranking.

As for merging (b > 0), the first order estimate of is given by(7)Similarly, one can write the estimate for the value of hierarchy of the inverted rankingIn this notation p, q, s, a are the parameters of the RSBM, while b refers to the modified ranking r(b) or r(i,b). Moreover it is clearly .

In the twitter-like hierarchy (p ≥ q > s) it is , i.e. the inverted ranking is never optimal. Merging, instead, can give rankings with higher hierarchy than the planted ranking.

To show this, in the left panel of Fig 3 we plot the behaviour of as a function of the number of classes, , after merging. Each line is associated to a RSBM(p, q, s, R). The parameters p = q = 0.5, R = 32 are fixed, while different curves refer to different values of s. We plot the variable as a continuous variable to help the interpretation of the observed behaviour. When s is small the maximum value of is correctly identified at . Above a critical value sm of the parameter describing the probability of a backward link, the planted ranking is no longer optimal and merging classes gives a ranking with higher hierarchy. Notice that for , the hierarchy of the planted ranking becomes negative. This might seem counterintuitive since we showed before that h* ∈ [0, 1]. The condition simply means that putting all the nodes in the same class has a higher hierarchy than the one of the planted ranking when .

[Figure omitted. See PDF.]

Fig 3.

Panel (a) shows the value of the estimate of h1 for different values of s as a function of the number of classes, , for twitter-like graphs with parameters p = q = 0.5, R = 32. Panel (b) gives a schematic representation of the estimated optimal number of classes as s varies.

https://doi.org/10.1371/journal.pone.0191604.g003

The right panel of Fig 3 shows the optimal number of classes as a function of s. As explained, when s < sm it is , while after this value the optimal number of classes decreases and in the limit s = smax it is . Therefore the value sm sets a resolution threshold, since twitter-like graphs with a probability of backward links larger than sm will not be correctly identified by agony with d = 1. More precisely sm is an upper bound of the resolution threshold, since other rankings, not considered here, could have higher hierarchy than the planted and the merged ones when s < sm.

Interestingly for large number of classes R, as we prove in the following Proposition 1, the resolution threshold scales as sm ∼ (6p − 3q)/R2, i.e. the more communities are present the more it is difficult to detect them. The same happens for large networks (N → +∞). Taking the number of classes constant and letting p and q scale as 1/N to keep the connectivity fixed, one immediately sees that sm = O(N−1), i.e. for large networks and fixed number of classes the detectable structures are those with very strong hierarchical structure. Thus agony with d = 1 has strong resolution limits for large graphs, similarly to what happens with modularity and community detection.

The situation is more complex in the military-like hierarchy (q = 0) because for large s inverted rankings become better than direct ones. To show this, we refer to the left panel of Fig 4, which is the analogous of left panel of Fig 3. In this case, alongside we also plot , with matching line colours to distinguish those associated to the same values of s, and circles to identify . In all cases we chose p = 0.5 and R = 32. For small values of s (solid blue lines), is convex in and has its maximum at , whereas is negative for inverted rankings different from the trivial one. Thus in this regime the planted ranking is optimal. When s reaches the critical value si (dashed red lines), the optimal choices for both the direct and inverted rankings give the same value of hierarchy. For higher s (dotted green lines) the only direct ranking with non negative hierarchy is the trivial one, i.e , while the inverted rankings are (strictly) positive for a suitable choice of b. Therefore in this regime inverted rankings outperfom the planted one.

[Figure omitted. See PDF.]

Fig 4.

Panel (a) shows how depending on the value of s the inverted rank can give a higher value of than the planted rank in military-like graph with parameters p = 0.5, q = 0, R = 32. Panel (b) gives a schematic representation of the estimated optimal number of classes as s varies, dashed lines are associated to the inverted rank.

https://doi.org/10.1371/journal.pone.0191604.g004

The right panel of Fig 4 shows the optimal number of classes as a function of s together with an indication of the sign of the hierarchy of the optimal direct and inverted ranking. For the hierarchy of the optimal direct ranking is positive and the one of the optimal negative ranking is negative, for they are both positive, while for s1 < s < smax the inverted optimal hierarchy is positive and the optimal direct one is negative. Thus for s < si the optimal ranking is direct and coincides with the planted one, while after this value the inverted ranking with two classes becomes optimal. This is true in the region after which the inverted ranking with three classes becomes optimal. By increasing s further, the optimal ranking is always inverted with an increasing number of classes up to a value smaller or equal to for s = smax. Therefore for the military-like hierarchy the resolution threshold is si which for large R scales as 6p/R2, displaying a resolution limit similar to the twitter-like hierarchy, both for large number of classes R and for large graphs (N → ∞).

We summarise the results for d = 1 in the following proposition.

Proposition 1 When d = 1 and p ≥ q > s, (Twitter hierarchy) the first order estimate for the optimal value of h (8) where Furthermore, when q = 0, (Military hierarchy) where The proof and the extended expression for are given in S1 Appendix.

In conclusion, we explicitly showed that for RSBMs there exist alternative rankings with a smaller agony (d = 1) than the planted one. The merging of the classes for the twitter hierarchy is due to fact that for a large number of classes it might be more convenient to aggregate classes paying a penalty equal to one than to leave them separate but paying a higher penalty for the distant backward links. Similarly, for the military hierarchy, when the number of backward links is relatively large, it is more convenient (in terms of agony) to invert the ranking because forward links do not enter the cost minimisation. Thus even if p is much larger than s and the number of forward links is much larger than the number of the backward links, it is more convenient to invert the ranking to avoid to pay large penalties of backward links between very distant classes.

Thus our results depend on the choice of the penalisation function and on the choice of the affinity matrix. In the next Subsection we show indeed that a very different result is obtained for d = 0. Changing the affinity matrix, for example introducing a probability of backward links which depends on the distance between classes, and changing the penalty function by including the negative cost of forward links is left for a future study.

Agony with d = 0.

This case corresponds to the FAS problem. The optimal ranking is obtained when each node is in a different class, , and the inverted ranking is never optimal as stated by the following:

Proposition 2 When d = 0, ∀RSBM(p, q, s, R = 2a) the optimal value for the first order estimate of h is given by (for both Twitter and Military hierarchy)See S1 Appendix for the proof. The reason for this result is that backward links are weighted in the same way irrespectively from the distance between the ranks of the nodes connected by the link. Thus, for example, the naive ranking with all nodes in one class has a agony equal to the number of links, while the ranking where each node is in one class has an agony equal to the number of backward links, which is smaller than the total number of links.

Finally we note that the value of increases very slowly when approaches N, so in specific realisations of the RSBM the optimal ranking can have a number of classes smaller than N.

Agony with d = 2.

Finally, we consider the case of d = 2. Similarly to the case d = 1, splitting is never optimal, both for the direct and inverted rankings, while merging can give rankings with higher value of than the planted one. One can proceed as before, considering the expressions for the alternative rankings when b > 0:and

As before we describe the behaviour for the two considered hierarchies and then we state the proposition summarising our results. For the twitter-like hierarchy (p ≥ q > s), the behaviour is similar to the d = 1 case. Since , ∀b, inverted rankings are never optimal. The planted ranking is optimal up to the critical value s2,m for the probability of backward links. After that, merged rankings outperform the planted one, and the number of classes decreases with s. When s2,1 < s ≤ smax the optimal choice is the trivial ranking, i.e. . Despite the similarity with the d = 1 case, the resolution threshold is now higher, since it can be shown that s2,m ≤ sm. Moreover, while, as noted before, in the d = 1 case sm = O(R−2), in the d = 2 case the resolution threshold is not only stricter but also it decreases faster as the number of classes increases, since it scales as . Finally, when d = 2 the large s case has the trivial ranking as the optimal one, whereas in the d = 1 case the optimal ranking has two classes.

For the military-like hierarchy (q = 0), the planted ranking is proven to be optimal with respect to the direct rankings up to the critical value . After this value the optimal choice is the trivial ranking. Then when it becomes optimal to merge inverted rankings and the optimal number of classes increases with s, starting from . Differently from the case d = 1, in this case it holds , hence for the optimal rank is the trivial one, and the resolution threshold is given by , which scales as 12p/R3, while inverted rankings are to be preferred for any .

We summarise the results for d = 2 in the following proposition.

Proposition 3 When d = 2 and p ≥ q > s (Twitter hierarchy), the first order estimate for the optimal value of h where and is given S1 Appendix.

Furthermore, when q = 0 (Military hierarchy), where

With this last proposition we showed that hierarchy detection with quadratic cost function has a behaviour very similar to the linear case. However the resolution limits we highlighted before escalates in this case, and, as a result, only very strong hierarchies are detected correctly when the number of class is large. The same computations can be done also for greater integers d, for which the sums in the estimates of agony have a closed formula. Intuitively as increases, backward links to distant classes are given a larger penalisation, hence rankings with merged classes become more convenient than the planted one even for smaller values of s. In other words agonies with d > 1 are strongly suboptimal and are able to identify very strong structures.

Following this remark, better candidates as penalty functions are likely those with 0 < d < 1. For at least some of those d one can expect to soften the resolution limits associated to integer d. However the approach to study the regime cannot rely on analytical formulae.

Numerical simulations

In this Subsection we show the results of numerical simulations to test the propositions we presented before. This is important for two reasons. First, to show that the guessed rankings, obtained by merging, splitting, or inverting the planted one, are indeed the optimal ones or have a hierarchy close to the optimal one. Second, to prove that the first order approximation and other simplifying assumptions give analytic expressions close to numerical simulations.

We use igraph [37] to sample a graph from the RSBM ensemble. For computing agony we use the algorithm described in [34], which we will refer to it as agony (in italics) for brevity, and which gives the exact solution for the optimisation problem when d = 1. Finally, we use the MCMC algorithm in the GraphTool [38] package for the inference of the SBM (without constraint on the structure of the affinity matrix).

We perform the same analysis with different choices for the parameters p, q, s, R, N and the results are consistent, thus in the following we present only representative cases. We use the adjusted Rand Index (RI) [39] to measure the similarity between the planted and the inferred ranking. The RI is 0 between independent rankings and 1 when each pair of elements that are in the same class in one ranking are also in the same class in the other. Ordering of classes does not matter in computing RI, thus the RI between a ranking and its inverted version is 1. Nevertheless we checked that high values of RI do not correspond to inverted rankings.

Twitter-like hierarchy.

We generate RSBM with parameters N = 212 = 4096, and we vary the value of s. Fig 5 shows the heat maps of the classes found by agony for different values of s. The heat-maps are constructed as follow: a square in position (i, j) refers to the number of nodes that belong to class i in the planted rank and are placed in class j by agony: the darker the colour, the higher the number. For small s (almost DAG structures) the algorithm recovers faithfully the planted ranking and the RI is high. When the hierarchical structure becomes weaker, the ranking obtained by agony is the merging of contiguous classes in the hierarchy, as postulated in the theoretical part above. For this choice of p, q, R the resolution threshold for s is sm = 0.00151 consistently with our simulations. As we predicted, classes merge more and more when s increases. The inferred rankings are close to uniform, and the main exception is the first and last class which are smaller than the other ones.

[Figure omitted. See PDF.]

Fig 5. Heat maps comparing the planted ranking with the ranking inferred with agony for twitter-like hierarchy.

In each panel a square in position (i, j) contains the number of nodes that belong to class i in the planted rank and are placed in class j by agony: the darker the colour, the higher the number. The parameters are p = q = 0.5, R = 32 and 9 values of s. Each plot refers to a single realisation from the ensemble.

https://doi.org/10.1371/journal.pone.0191604.g005

We show numerically that the ranking we proposed as optimal in the previous Section has indeed a value of hierarchy very close to the one obtained from simulations. In Fig 6 we show a scatter plot of the true value of computed with agony on the simulated graphs against the hierarchy of the planted rank (circles), and against , the hierarchy computed with Eq (8) (stars). To evaluate the latter we use the coefficients of the RSBM estimated from the sample graph with GraphTool. We estimate p = q and s as the average elements of the inferred affinity matrix on the corresponding classes and we leave free the number of classes. For s < sm (red symbols) the two methods agree and give a value of hierarchy consistent with the real one. When s > sm (green and blue symbols depending on whether s is smaller or larger of ) the hierarchy of the planted ranking is significantly smaller than , showing that another ranking is optimal. This has a value of hierarchy which is very close to the one computed from Eq (8), even when the coefficients of the RSBM are estimated from data. It is interesting to note that this is true also for s very close to smax where the number of classes detected by GraphTool is significantly smaller than R. This is due to the fact that the analytical expression in Eq (8) of the value of hierarchy of the merged ranking depends weakly on the number of classes. This is a strong indication that the ranking we suggested, and obtained by merging the classes, has a value of hierarchy which is indeed very close to the globally optimal one. In conclusion, the planted hierarchy is optimal for a very small range of values of s and, as we expected, it gives negative values of h1 when s is large enough. On the other side, our estimate for optimal h1 is accurate for all the value of s considered.

[Figure omitted. See PDF.]

Fig 6. Comparison of hierarchies for twitter-like RSBMs.

The parameters are p = q = 0.5, R = 32, s varies in [0, smax], with smax = 0.448. Each point refers to a single realisation of the ensemble. The circles represent the pairs , i.e. the optimal hierarchy computed with agony and the one of the planted hierarchy . The stars represent where is the theoretical hierarchy of Eq (8) with the parameters of the SBM estimated via GraphTool. Finally, sm is the theoretical resolution threshold and is the theoretical value of s for which the estimate for the planted hierarchy is zero.

https://doi.org/10.1371/journal.pone.0191604.g006

Finally in Fig 7 we show that the resolution problem is due to the choice of the method (agony with d = 1) and not necessarily to the model itself. In fact it is well known that SBM have a resolution threshold both when inference is done using Maximum Likelihood methods [40] and spectral methods [41]. To this end we infer a SBM on the adjacency matrix, keeping free the number of classes (see [38] for the model selection adopted by GraphTool) and we compute the RI of the planted ranking versus the one obtained with agony and the SBM fit. The result is shown in Fig 7 for different values of s. We see that the SBM fit outperforms agony. It is clear that, since we are using SBM for generating the graphs, its fitting will be better. However what we want to stress is that there is remarkably wide interval of values of s for which agony is not able to detect a hierarchical structure even if it is strong enough to be detected by another method. Hence the limit in resolution is not embedded in the RSBM but in the objective function associated to agony.

[Figure omitted. See PDF.]

Fig 7. The figure shows the value of the Rand Index between the planted ranking and the inferred ones.

The blue squares considers the ranking obtained with agony (hence d = 1), while the red triangles considers the ranking obtained with a RSBM fit via GraphTool. The parameters of the twitter-like hierarchy are p = q = 0.5, R = 32, s varies in [0, smax], with smax = 0.448, and each point refers to a single realisation of the ensemble.

https://doi.org/10.1371/journal.pone.0191604.g007

Military-like hierarchy.

For the military-like hierarchy things are more complicated. Fig 8 shows the heat map of the classes for p = 0.5 and nine values of s. With these parameters our formulas give si = 0.00280 and s1 = 0.00284. We see that for strong hierarchical structures (small s) agony recovers well the classes. However when s increases a partial inversion of the hierarchy is observed and only for large s we recover the fully inverted ranking we studied in the previous Section. Thus simulations show that the latter is not always the optimal ranking but rather there are partially inverted rankings with a larger hierarchy. The purpose of the above analysis on the military-like hierarchy is to show that there exist values of the parameters for which the planted ranking is not optimal and to demonstrate that partial inversion can outperform the planted one. Moreover the partial inversion is observed for s = 0.002 < s1, hence our computations provide a upper bound of the true resolution threshold.

[Figure omitted. See PDF.]

Fig 8. Heat maps comparing the ranking inferred using agony with the planted ranking for military-like hierarchy.

In each panel a square in position (i, j) contains the number of nodes that belong to class i in the planted rank and are placed in class j by agony: the darker the color, the higher the number. The parameters are p = 0.5, q = 0, R = 32, s varies in [0, smax], with smax = 0.0294, and each plot refers to a single realisation of the ensemble.

https://doi.org/10.1371/journal.pone.0191604.g008

Fig 9 shows, similarly to Fig 6, the scatter plot of the true value of computed via agony on the simulated graphs against the hierarchy of the planted rank (circles), and against , the hierarchy computed with Eq (8) using the coefficients of the SBM estimated from the sample graph with GraphTool. The main message of the Fig is that, despite the fact the symmetrically inverted ranking is not the optimal one according to numerical simulations, its value of hierarchy is very close to the one of the optimal ranking, while the planted one strongly mis-estimates the value of h. Thus our computation in the previous Section can be used to reliably estimate the hierarchy of a military-like ranking. This is obviously a partial answer and analytical calculations of the hierarchy of partially inverted rankings are left for a future study.

[Figure omitted. See PDF.]

Fig 9. Comparison of hierarchies of military-like RSBMs.

The parameters are p = 0.5, q = 0, R = 32, s varies in [0, smax], with smax = 0.0294, and each point refers to a single realisation of the ensemble. The circles represent the pairs , i.e. the optimal hierarchy computed with agony and the one of the planted hierarchy . The stars represent where is the theoretical hierarchy with the parameters of the SBM estimated via GraphTool.

https://doi.org/10.1371/journal.pone.0191604.g009

Beyond the resolution limit: Iterated agony.

In the previous Sections we have shown theoretically and numerically that inference of ranking hierarchies based on agony suffers from significant resolution limit. In twitter-like hierarchies, the identified classes are merging of adjacent classes and thus small classes are not identified. In military-like hierarchies inversions start to play a significant role.

An heuristic method to overcome this problem is to iterate the application of agony. As done with modularity, one can apply agony to each class found in the first iteration of the algorithm, in order to find subclasses. In principle one could continue to iterate, even if the fact that agony finds two classes in an Erdös-Renyi graph suggests a careful design of the stopping criterion. The purpose of this Section is not to propose a full criterion for the improvement of agony via iteration, but to show that indeed improvement is possible, both considering model graphs and real networks.

We first consider the model graphs with twitter-like hierarchy we presented in the previous Section. Fig 10 shows the RI between the planted ranking and the one inferred with one (as in the previous Section) and two iterations of agony with d = 1. For small values of s the second iteration does not improve the inference because one iteration already recovers the planted structure. For larger values of s, i.e. weaker structures, the second iteration dramatically outperforms the result of the first one, indicating that iterated applications of agony can significantly improve the hierarchies detection. Table 1 shows some details of the obtained results. It is worth noticing that the value of h after the second run is actually smaller than the one from the first run, despite the fact that the RI follows the opposite pattern. This is expected since agony finds the optimal value of h, while the RI looks at the similarity with the planted ranking. A closer look to the results of the two iterations (see S1 Table) highlights that high number of classes after the second iteration and high hierarchy in each subclass are associated to the cases for which there is no significant improvement in the RI, hence a successful routine would rely on the control of these two quantities to stop the iterations.

[Figure omitted. See PDF.]

Fig 10. Comparison of the Rand Index between the planted ranking and one (blue squares) or two (orange triangles) iterations of agony.

Data refers to simulation of twitter-like HSBM with parameters p = q = 0.5, R = 32, s ∈ [0, smax], with smax = 0.448, and each point refers to a single realisation of the ensemble.

https://doi.org/10.1371/journal.pone.0191604.g010

[Figure omitted. See PDF.]

Table 1. Simulated graphs, output of the two runs of agony.

https://doi.org/10.1371/journal.pone.0191604.t001

We now show that the same phenomenon is relevant also for real networks. We investigate four datasets from SNAP, Stanford Network Analysis Platform [42], which were also used in [34]. Note that these datasets have been updated since they have been used in [34] so our results are slightly different.

The networks are quite different in size (from a minimum of 7K nodes to almost 400K nodes) but they are all quite sparse.

* Wiki vote. The network contains all the Wikipedia voting data from the inception of Wikipedia till January 2008. Nodes in the network represent Wikipedia users and a directed edge from node i to node j represents that user i voted for user j.

* Higgs Reply. The network contains replies to existing tweets: nodes are users and i is linked to j if i replied to a j’s tweet.

* Higgs mention. Similar to the previous case, here links represent mentions: a link from i to j means that user i mentioned user j.

* Amazon. Network was collected by crawling the Amazon website. It is based on Customers Who Bought This Item Also Bought feature of the Amazon website. If a product i is frequently co-purchased with product j, the graph contains a directed edge from i to j.

Table 2 reports some properties of the networks alongside the output of one and two iterations of the agony algorithm. Specifically, for each network the table contains: the number of nodes N, the density (, where m is the number of edges), the percentage of nodes in the largest strongly connected component (SCC), the value of , the number of classes inferred in the first run (R) and the total number of classes after the second run (R′) of agony.

[Figure omitted. See PDF.]

Table 2. Networks summary.

SCC is the percentage of nodes in the largest strongly connected component, is the hierarchy of the ranking obtained with one iteration of agony, R is the number of classes in the globally optimal ranking, and R′ is the number of classes after two iterations of agony.

https://doi.org/10.1371/journal.pone.0191604.t002

It is clear that the second application the algorithm to the classes detected in the first iteration increases significantly the number of classes, suggesting that the classes identified in the first iteration could be aggregation of smaller classes. In S2 Table we report more details on the classes identified in the iteration and on the subclasses identified by the second iteration.

Since agony penalises links among nodes in the same class, the subgraphs in some cases have no links (those with ∗ in S2 Table. Notice this would be the case for any class in a DAG. Thus, a low value of h in each class and a number of sub classes larger than 2 indicate a non trivial and not completely resolved structure of the class.

Conclusion

In this paper we have studied the inference of hierarchical structures in directed networks by introducing an ensemble of random graphs, termed the Hierarchical Stochastic Block Model, and studying how agonies, penalising links contrary to the hierarchy, are able to identify the planted ranking.

Using symmetry arguments we have explored ranking alternative to the planted one and obtained from it by merging, splitting or inverting its classes. We have shown that when the hierarchy is not strong enough some of these alternative rankings of nodes have a value of the hierarchy larger than the planted one. This demonstrates that (generalised) agonies have a resolution limit, being unable to detect small classes in large networks. This is somewhat similar to the well known resolution limit of modularity in community detection. In some cases we have strong numerical indications that the proposed alternative rankings, are actually close to the global optimal one. Finally we have shown that in these cases the iterated application of agony can lead to significant improvement of the hierarchy detection.

There are several directions along which our work can be extended. First, we have investigated in detail the case of uniform cardinality of the classes, even if our formulae can be used to study more complex structures, such as a pyramidal hierarchy with a small top class and larger bottom classes. With a careful choice of the sizes one might be able to maintain analytical tractability, however the study of these structures are left for future investigation. The second direction is to consider, at least theoretically, other values of d (or other agony functions). We have shown some results indicating that the resolution threshold depends on d, however numerical simulations cannot be performed because of the lack of heuristic methods for optimisation of agony with d ≠ 0, 1. Finally, other methods to identify ranking hierarchies could be investigated, for example suitably modifying the agony function or by considering optimisations for a set of functions.

We leave these extensions for future work and we are confident that the results will be of interest in the general problem of hierarchy detection in networks.

Supporting information

[Figure omitted. See PDF.]

S1 Appendix. Detailed proofs.

https://doi.org/10.1371/journal.pone.0191604.s001

(PDF)

S1 Table. Numerical results for simulated graphs.

https://doi.org/10.1371/journal.pone.0191604.s002

(PDF)

S2 Table. Numerical results for real graphs.

https://doi.org/10.1371/journal.pone.0191604.s003

(PDF)

Acknowledgments

We would like to thank Ilaria Bordino, Francesco Gullo, Marcello Paris, and Daniele Tantari for useful discussions. We also thank Nikolai Tatti for providing us the agony code.

Citation: Letizia E, Barucca P, Lillo F (2018) Resolution of ranking hierarchies in directed networks. PLoS ONE 13(2): e0191604. https://doi.org/10.1371/journal.pone.0191604

References

1. Johnson S, Domínguez-García V, Donetti L, Muñoz MA. Trophic coherence determines food-web stability. Proceedings of the National Academy of Sciences. 2014;111(50):17923–17928.

2. Shetty J, Adibi J. Discovering Important Nodes through Graph Entropy: The Case of Enron Email Database. Proceedings of the 3rd International Workshop on Link Discovery. 2005; p. 74–81.

3. Nguyen H, Zheng R. A data-driven study of influences in Twitter communities. In: 2014 IEEE International Conference on Communications (ICC). IEEE; 2014. p. 3938–3944.

4. Fama EF, French KR. Testing trade-off and pecking order predictions about dividends and debt. Review of financial studies. 2002;15(1):1–33.

5. Frank MZ, Goyal VK. Testing the pecking order theory of capital structure. Journal of financial economics. 2003;67(2):217–248.

6. Glattfelder JB, Battiston S. Backbone of complex networks of corporations: The flow of control. Physical Review E. 2009;80(3):036104.

7. Tatti N. Hierarchies in directed networks. In: Data Mining (ICDM), 2015 IEEE International Conference on. IEEE; 2015. p. 991–996.

8. Orlin JB. A faster strongly polynomial minimum cost flow algorithm. Operations research. 1993;41(2):338–350.

9. Fortunato S. Community detection in graphs. Physics Reports. 2010;486:75–174.

10. Kumpula JM, Saramäki J, Kaski K, Kertész J. Limited resolution in complex network community detection with Potts model approach. The European Physical Journal B. 2007;56(1):41–45.

11. Fortunato S, Barthelemy M. Resolution limit in community detection. Proceedings of the National Academy of Sciences. 2007;104(1):36–41.

12. Blondel V, Guillaume J, Lambiotte R, Lefebvre E. Fast unfolding of communities in networks. JSTAT. 2008; p. P10008.

13. Simon HA. The architecture of complexity. Facets of systems science. 1991; p. 457–476.

14. Krackhardt D. Graph theoretical dimensions of informal organizations. Computational organization theory. 1994;89(112):123–140.

15. Maiya AS, Berger-Wolf TY. Inferring the maximum likelihood hierarchy in social networks. In: Computational Science and Engineering, 2009. CSE’09. International Conference on. vol. 4. IEEE; 2009. p. 245–250.

16. Gupte M, Shankar P, Li J, Muthukrishnan S, Iftode L. Finding hierarchy in directed online social networks. In: Proceedings of the 20th international conference on World wide web. ACM; 2011. p. 557–566.

17. Romei A, Ruggieri S, Turini F. The layered structure of company share networks. In: Data Science and Advanced Analytics (DSAA), 2015. 36678 2015. IEEE International Conference on. IEEE; 2015. p. 1–10.

18. Trusina A, Maslov S, Minnhagen P, Sneppen K. Hierarchy measures in complex networks. Physical review letters. 2004;92(17):178702. pmid:15169201

19. Tibély G, Pollner P, Vicsek T, Palla G. Extracting tag hierarchies. PloS one. 2013;8(12):e84133. pmid:24391901

20. Nepusz T, Vicsek T. Hierarchical self-organization of non-cooperating individuals. PloS one. 2013;8(12):e81449. pmid:24349070

21. Corominas-Murtra B, Goñi J, Solé RV, Rodríguez-Caso C. On the origins of hierarchy in complex networks. Proceedings of the National Academy of Sciences. 2013;110(33):13316–13321.

22. Mones E. Hierarchy in directed random networks. Physical Review E. 2013;87(2):022817.

23. Clauset A, Moore C, Newman MEJ. Structural Inference of Hierarchies in Networks. Statistical Network Analysis: Models, Issues, and New Directions Lecture Notes in Computer Science, Vol 4503. 2007; p. 1–13.

24. Clauset A, Moore C, Newman MEJ. Hierarchical structure and the prediction of missing links in networks. Nature. 2008;453(7191):98–101. pmid:18451861

25. Newman MEJ. Communities, modules and large-scale structure in networks. Nature Physics. 2012;8(1):25–31.

26. Peixoto TP. Hierarchical Block Structures and High-Resolution Model Selection in Large Networks. Phys Rev X. 2014;4:011047.

27. Cucuringu M. Sync-Rank: Robust Ranking, Constrained Ranking and Rank Aggregation via Eigenvector and SDP Synchronization. IEEE Transactions on Network Science and Engineering. 2016;3(1):58–79.

28. Langville AN, Meyer CD. Who’s# 1?: the science of rating and ranking. Princeton University Press; 2012.

29. Bickel PJ, Chen A. A nonparametric view of network models and Newman–Girvan and other modularities. Proceedings of the National Academy of Sciences. 2009;106(50):21068–21073.

30. Zhang P, Moore C. Scalable detection of statistically significant communities and hierarchies, using message passing for modularity. Proceedings of the National Academy of Sciences. 2014;111(51):18144–18149.

31. Newman MEJ. Community detection in networks: Modularity optimization and maximum likelihood are equivalent. arXiv preprint arXiv:160602319. 2016;.

32. Holland PW, Laskey KB, Leinhardt S. Stochastic blockmodels: First steps. Social networks. 1983;5(2):109–137.

33. Newman MEJ. Networks: An introduction; 2014.

34. Tatti N. Tiers for peers: a practical algorithm for discovering hierarchy in weighted networks. Data Mining and Knowledge Discovery. 2017;31(3):702–738.

35. Slater P. Inconsistencies in a schedule of paired comparisons. Biometrika. 1961;48(3/4):303–312.

36. Eades P, Lin X, Smyth WF. A fast and effective heuristic for the feedback arc set problem. Information Processing Letters. 1993;47(6):319–323.

37. Csardi G, Nepusz T. The igraph software package for complex network research. InterJournal. 2006;Complex Systems:1695.

38. Peixoto TP. The graph-tool python library. figshare. 2014;

39. Hubert L, Arabie P. Comparing partitions. Journal of Classification. 1985;2(1):193–218.

40. Decelle A, Krzakala F, Moore C, Zdeborová L. Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications. Physical Review E. 2011;84(6):066106.

41. Nadakuditi RR, Newman MEJ. Graph spectra and the detectability of community structure in networks. Physical Review Letters. 2012;108(18):1–5.

42. Leskovec J, Krevl A. SNAP Datasets: Stanford Large Network Dataset Collection; 2014. http://snap.stanford.edu/data.