{'fraction_non_alphanumeric': 0.03195182946075218, 'fraction_numerical': 0.009638963750698035, 'mean_word_length': 4.6468330134357005, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 2, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Validation of community robustness

The large amount of work on community detection and its applications leaves unaddressed one important question: the statistical validation of the results. In this paper we present a methodology able to clearly detect if the community structure found by some algorithms is statistically significant or is a result of chance, merely due to edge positions in the network. Given a community detection method and a network of interest, our proposal examines the stability of the partition recovered against random perturbations of the original graph structure. To address this issue, we specify a perturbation strategy and a null model to build a set of procedures based on a special measure of clustering distance, namely Variation of Information, using tools set up for functional data analysis. The procedures determine whether the obtained clustering departs significantly from the null model. This strongly supports the robustness against perturbation of the algorithm used to identify the community structure. We show the results obtained with the proposed technique on simulated and real datasets.


Introduction
Networks are mathematical representation of interactions among the components of a system and can be modeled by graphs. A graph G=(V,E) consists of a collection of vertices V, corresponding to the individual units of the observed system, and a collection of edges E, indicating some relation between pairs of vertices.
Graphs modeling real systems, i.e. social, biological, and technological networks, display non trivial topological features. Indeed they present the properties that define a complex network: big inhomogeneities, a broad degree distribution and distribution of edges locally inhomogeneous. In the study of complex networks, a network is said to have a community structure if the vertices can be divided in g groups, such that nodes belonging to the same group are densely connected and the number of edges between nodes of different groups is minimal.
The problem of community detection (graph partitioning) has been widely studied by researchers in a variety of fields, including statistics, physics, biology, social and computer science in the last 15 years. Finding communities within an arbitrary complex network can be a computationally difficult task. The number of communities, if any, within the network is typically unknown and the communities are often of unequal size and/or density. Despite these difficulties, however, several methods for community finding have been developed and employed with varying levels of success, see [7], [11], [13], [15], [18] and [29] for reviews.
Our work focusses on the problem of testing the robustness of the recovered partition of a given community detection method. In the following we provide a brief review of the state of the art of the literature addressing this problem. Although the huge work developed for community detection and its applications, the question of the significance of results still remains open. Our proposal represents a first attempt to statistically define the robustness of a clustering and hence can not be directly compared to any of the following described methodologies.
State of the art The modularity Q of Newman and Girvan [24] was the first attempt to give an answer to this question. It is defined as the fraction of the edges that fall within the given groups minus the expected such fraction if edges were distributed at random and is based on the idea that a random graph is not expected to have a cluster structure. However, as pointed out in [11] and [17], there is an important limit. Precisely, networks with a strong community structure have high modularity, on the contrary high modularity does not imply networks with a community structure. Other authors, see [14], [31], suggested the use of a z-score to compare the maximum modularity of a graph to the maximum attainable modularity in purely random graphs of the same size and expected degree sequence. The problem is that the distribution of the null model, though peaked, is not Gaussian, causing false positives and false negatives.
A different approach was developed in [21], where the authors studied how canonical ensembles of network partitions depend on temperature to assess the significance and nature of the community structure obtained by algorithms that optimize the modularity. In this case −Q plays the role of energy, i.e. at temperature T , the statistical weight of a given partition in the ensemble is proportional to exp(Q/T ). Typically, as the temperature increases, there is a transition from low entropy/high Q partitions (significant cluster structure) to high entropy/low Q (random partitions). If there is strong community structure, the transition is sharp. The peak is broader for networks with weaker community structure, as there are more reasonable alternative partitions with intermediate values of Q, and so the transition occurs over a broader range of temperature. They also introduced an order parameter to measure the similarity of the sampled partitions at a given temperature, i.e. whether there is just a single partition with high Q or a number of competing partitions. Therefore, it is a useful tool to detect false positives. However the methodology is computationally onerous and cannot be easily generalized to other optimization methods.
In [4] the authors introduced the notion of entropy of graph ensembles to assess the relevance of additional information about the nodes of a network using the information that comes from the topology of the network itself. The indicator of clustering significance Θ they introduce can also reveal statistical regularities that shed light on possible mechanisms underlying the network stability and formation.
In [20] the authors presented the Order Statistics Local Optimisation Method (OSLOM ), a technique based on the local optimization of a fitness function, the C-score [19], expressing the statistical significance of a cluster with respect to random fluctuations. Given a subgraph C in a graph G, the C-score measures the probability that the number of links connecting a node to nodes in C, where C is embedded within a random graph, is higher than or equal to the value seen in the original graph G. This score permits to rank all the vertices external to C (in increasing order of the C-score), having at least one connection with C, and to calculate its order statistic distribution Ω. The minimum of Ω is the random variable whose cumulative is the score of the community C. To asses its significance a threshold parameter P is fixed. The procedure is iterated to analyse the full network. The novelty of this approach is the local estimate of the significance, i.e. of single communities, not of partitions; on the contrary a serious limit is due to the lack of a data driven procedure to estimate P , indeed the authors fix its value to 0.1.
Recently, [38] proposed a testing based community detection procedure called Extrac-tion of Statistically Significant Communities (ESSC). The ESSC procedure measures the statistical significance of connections between a single vertex and a set of vertices in undirected networks under a null distribution derived from the configuration model [2]. Given an observed network G 0 with n vertices and a vertex set B, they introduce the statisticŝ d(u : B), measuring the number of edges between a vertex u and B in the random modelĜ, and show thatd n (u : B) is approximately binomial as n → ∞ in the total variation distance between two probability mass functions. This permits to obtain the p-values of the null distribution using the binomial approximation and gives origin to an iterative deterministic procedure that recovers robust communities. The technique has some similarities with OSLOM, indeed both are extraction methods and use the configuration model as reference distribution, but differentiates because the probabilities have a closed form. Another group of techniques was proposed in [12], [17], [32], and their conceiving was completely different from previous described methodologies. Indeed they introduce a stochastic component in the network by perturbing the graph structure, measure the effect of the perturbation and compare it with the corresponding value for a null model graph. The basic idea is that a significant partition should not be altered by small modifications, as long as the modification is not too extensive. An interesting feature of these methods is their independence from the community detection technique adopted.
In this paper we present a methodology able to clearly detect if the community structure found by some algorithms is statistically significant or is a result of chance, merely due to edge positions in the network. Given a community detection method and a network of interest, our proposal examines the stability of the partition recovered against random perturbations of the original graph structure. To address this issue, following ideas from [17], we specify a perturbation strategy and a null model to build some procedures based on Variation of Information as stability measure. Given a this measure we address the question of evaluating its the significance. This permits to build the Variation of Information curve as a function of the perturbation percentage and to compare it with the corresponding null model curve using analysis tools set up for functional data analysis.
The rest of the paper is organised as follows. In Section 2 we introduce the proposed procedures based on Variation of Information and the functional data analysis techniques, including their detailed description. Section 3 shows the results achieved applying our methodology on simulated and real datasets. Conclusions and ideas for future research are drawn in Section 4.

Variation of Information hypothesis testing
Variation of Information (VI) is an information theoretic criterion for comparing two partitions, or clusterings, of the same data set [23]. It is a metric and measures the amount of information lost and gained in changing from clustering C to clustering C . The criterion makes no assumptions about how the clusterings were generated and applies to both soft and hard clusterings.
Given a dataset D and two clusterings C and C of D, with K and K non empty clusters, respectively, VI is defined as where H (C) is the entropy associated with clustering C and I (C, C ) is the mutual information between C and C , i.e the information that one clustering has about the other P (k) is the probability of a point being in cluster C k and P (k, k ) is the probability that a point belongs to C k in clustering C and to C k in C . Another equivalent expression for VI is The first term measures the amount of information about C that we loose, while the second measures the amount of information about C that we gain, when going from clustering C to clustering C . VI metric is the basis of the hypothesis testing procedures we propose to establish the statistical significance of a recovered community structure in a complex network. Our original idea is to generate two different curves based on the VI measure and to statistically test their difference. The first curve V Ic is obtained computing VI between the partition of our original network and the partition of different perturbed version of our original network. The second curve V Ic random is obtained computing VI between the partition of a null random network and the partition of different perturbed version of such null network. The comparison between the two VI curves turns the question about the significance of the retrieved community structure into the study of the stability/robustness of the recovered partition against perturbations. We expect that it must be robust to small perturbations, because if "small changes" in the network imply a completely different partition of the data, it means that the found communities are not trustworthy, and this cannot be due to the failure of the chosen algorithm for the community detection. Indeed the proposed testing procedure is independent from the clustering algorithm and it is easy to check if such a behaviour is due to it.
To understand well this point we must consider the behaviour of the VI curve for networks having a real community structure and those having a very poor community structure. In the first case the VI curve starts at 0, when the perturbation level p is 0% (unperturbed graph), rises rapidly (perturbation level between 0% and 40%), then levels off when 50% < p < 100%; in the second case the VI curve immediately grows up to a certain value and levels off that value, meaning that whatever partition has been found, at each level of perturbation, it has the same robustness of a random graph. Obviously the set up of a testing procedure is more necessary for all cases where the community structure is moderate or weak and the behaviour of the VI curve could be similar to that of a random graph.
The choice of the null random network is more delicate, because we expect that it has the same structure of our original graph but with completely random edges. This is why our choice relapses on the Configuration Model [2] associated with the degree sequence of the observed graph d = {d(1), . . . , d(n)} with vertex set V = {1, . . . , n}, i.e. CM(d). The CM(d) is a probability measure on the family of multi-graphs with vertex set V and degree sequence d that reflects, within the constraints of the degree sequence, a random assignment of edges between vertices. The generative form is simple: one can simply cut all the edges in the network, so every node still retains its degree by the number of half-edges or stubs emanating from it. The result will be an even number of half-edges. To create new networks with the same degree, one simply needs to randomly pair all the half-edges, creating the new edges in the network. The Configuration Model generates every possible graph with the given degree distribution with equal probability. Note that it naturally creates networks with multiple edges between nodes and self-connections between nodes. If such networks are unacceptable, one can reject those samples and try the algorithm again, repeating until one obtains a network without multiple or self-connections. The perturbation strategy adopted is described in Section 2.4. The basic steps of our method are the following: Overall Procedure 1. find a partition C of the given network, with vertex set V = {1, . . . , n} and degree sequence d = {d (1), . . . , d(n)}, by some chosen method M ; 2. given a perturbation level p ∈ [0, 1], perturb the network shuffling its edges by the percentage p preserving the original graphs degree distribution; 3. using the same method M find a partition C for the perturbed network and compute the V I(C, C ) to compare C and C ; 4. repeat steps 2-3 at different perturbation levels p ∈ [0, 1] so to obtain the V Ic curve. ( Note that p = 0 corresponds to the original unperturbed graph, p = 1 corresponds to the maximal perturbation level (random graph)); 5. repeat steps 1-4 starting in step 1 from the CM(d) (null model) so to get the V Ic random curve; 6. compare the V Ic and V Ic random curves as functions of p, statistically testing "the difference" between the two curves.
Note that the variation of p from 0 to 1 induces an intrinsic order to the data structure as in temporal data and can be treated as time point. Moreover, as it will be described in Section 2.4, we generate many perturbed graphs (i.e. 10) for each different level of p and these are considered as replicates per time points in our strategy.
Step 6 of the above procedure is achieved by a functional data analysis approaches aiming to test if the two groups of curves represent "the same process" or "different processes". The testing procedure we rely on is based on a tool set up for time course microarray data, namely Gaussian Process (GP) regression [16]. Aim of the GP regression in the context of gene expression data is to identify differentially expressed genes in a one-sample time course microarray experiment, i.e. to detect if the profile has a significant underlying signal or the observations are just random fluctuations. In this case we reformulate the testing problem working on log 2 V Ic V Ic random , as described in Section 2.1.
In order to show that our approach is robust with respect to the testing procedure used in Step 6, we also display the overall results achieved when using other two approaches in Step 6, described respectively in Section 2.2 and Section 2.3. Indeed, we can look at the two measured V I curves as independent realisations of two underlying processes say X 1 and X 2 observed with noise on a finite grid of points p ∈ [0, 1] and to test the null hypothesis versus the alternative hypothesis where d = means that the processes on either side have the same distribution. Then, as described in Section 2.2, taking advantage of the Karhunen-Loève expansion we explore the methodology developed in [28] based on Functional Principal Components Analysis (FPCA) to test (5).
On the contrary, the approach described in Section 2.3 addresses a domain-selective inferential procedure, providing an interval-wise non parametric functional testing [27], able not only to assess (5), but also to point out specific differences.
We will briefly describe GP regression, FPCA and interval-wise functional testing in the following sections. We would like to point out that our overall procedure provides a workflow to validate a community structure under different perspectives that can be investigated in dependence of the specific real problem dealt with. The description of the three different testing procedure is functional to the understanding of our overall procedure and in particular how we exploit the theory underling each single methodology to compare the curves V Ic and V Ic random . Hence we will summarise the three testing procedures to highlight the key connection to our testing problem. We choose to provide a review of each single methodology to provide awareness of the differences between the three procedures and their link to our testing problem. Even if addressing the same problem they are not equivalent. We refer to the original papers for any theoretical property of such testing procedures, including type I/II error study. We want to stress that the original contribution of our proposal is summarised in the six steps procedure depicted in the above frame.

GP regression
In this section we briefly summarise the methodology proposed in [16], where the authors present an approach to estimate the continuous trajectory of gene expression time-series from microarray through GP regression.
Briefly we recall that a gaussian process is the natural generalisation of a multivariate Gaussian distribution to a Gaussian distribution over a specific family of functions. More precisely, as defined in [30], a gaussian process is a collection of random variables, any finite number of which have a joint Gaussian distribution and is completely specified by its mean function and its covariance function. If we define the mean function m(x) and the covariance function k(x, x ) of a real process f (x) as: then we can write the GP as The Since observations are noisy, i.e.
, assuming that the noise ε ε ε ∼ N (0, σ 2 n I) and using Eq. (7)-(8), the marginal likelihood with K y = σ 2 w ΦΦ T + σ 2 n I. In this framework the hypothesis testing problem can be reformulated, over the perturbation interval [0, 1], as: The marginal likelihood derived from Eq. (10), enables then to compare or rank different models by calculating the Bayes Factor (BF). More specifically the BF is approximated with a log-ratio of marginal likelihoods of two GPs, each one representing the hypothesis of differential (the profile has a significant underlying signal) and non differential expression (there is no underlying signal in the profile, just random noise). The significance of the profiles is then assessed based on the BF.

Functional Principal Component testing
In this section we briefly summarise the approach proposed in [28] to test the hypothesis (5) when the observed data are realisations of the curves at finite grids and possibly corrupted by noise. Their motivating application is a diffusion tensor imaging study, where the objective is to compare white matter track profiles between healthy individuals and multiple sclerosis patients. They introduce a novel framework based on functional principal component analysis (F P CA) of an appropriate mixture process, referred to as marginal F P CA. The statistical framework for this problem assumes to observe data arising from two groups, namely , where t 1ij , t 2ij ∈ T , a compact interval that in our case is T =[0,1] (time plays the role of perturbation level p). It is assumed that the Y 1ij and the Y 2ij are independent realisations of two underlying processes observed with noise on a finite grid of points: where X 1ij ∼ IID X 1 (·) and X 2ij ∼ IID X 2 (·) are independent and square integrable random functions over T , for some underlying (latent) random processes X 1 and X 2 . It is assumed that X 1 and X 2 are second-order stochastic processes with mean functions assumed to be continuous and covariance functions assumed to be continuous and positive semidefinite, both being unknown. The measurement errors { 1ij } and { 2ij } are independent and identically distributed (IID), with zero mean and variances σ 2 1 and σ 2 2 , respectively. The authors exploit the truncated Karhunen-Loève expansion of the mixture process X(·) of X 1 (·) and X 2 (·) with mixture probabilities p and 1 − p. Let Z a binary random variable taking values in {1, 2} with P (Z = 1) = p, then X 1 (·) = E [(·)/Z = 1] and X 2 (·) = E [(·)/Z = 2]. Let us consider the truncated Karhunen-Loève expansion of X(·) and define X K Z (t) = µ(t) + K k=1 ξ Zk Φ k (t), Z = 1, 2, testing hypothesis (5) reduce to testing if the F P C scores {ξ k 1 } K k=1 and {ξ k 2 } K k=1 have the same distribution: In practice the authors consider K null hypothesis given the finite truncation level and propose a multiple two-sample univariate test, the Anderson-Darling (AD) statistic [26], combined with a multiple-comparison adjustment. The authors propose a Bonferroni correction, a procedure which controls the probability of erroneously rejecting even one of the true null hypotheses, the Family Wise Error Rate (FWER). In this case hypothesis (12) is rejected if min where p k is the p-value that is obtained by using the chosen univariate two-sample test for each H k 0 . The false discovery rate (FDR), suggested in [3] is a different point of view for how the errors in multiple testing could be considered. The FDR is the expected proportion of erroneous rejections among all rejections. If all tested hypotheses are true, controlling the FDR controls the traditional FWER. But when many of the tested hypotheses are rejected, indicating that many hypotheses are not true, the error from a single erroneous rejection is not always as crucial for drawing conclusions from the family tested, and the proportion of errors is controlled instead. Using the individual testing statistics proposed in [28] we will therefore adopt this FDR approach to adjust our tests for multiplicity.
Note that this procedure is designed for a more general framework in which the two curves V I and V Ic can be observed at different time points (i.e. p ∈ [0, 1]).

Interval-wise Functional testing
In the following we will briefly review the Interval-wise Functional testing procedure (IT P ) proposed by [27], where the authors develop a non-parametric domain-selective inferential methodology for functional data embedded in the L 2 (T ) space (where T is any limited open interval of R) to test (5). Their technique is not only able to assess the equality in distribution between functional populations, but also to point out specific differences. They propose a procedure based on the following three steps: More in detail, let us assume to observe two independent samples of sizes n 1 and n 2 of independent random functions on a separable Hilbert space y ij (t), i = 1, . . . , n j , j = 1, 2.
In the first step, data are projected on a finite-dimension subspace generated by a reduced basis y ij (t) = p k=1 c ij Φ (k) (t), where integer p represents the dimension. It follows that each of the n = n 1 + n 2 units can be represented by means of the corresponding p coefficients {c In the second step, the authors build a family of multivariate tests for . . , p and k is a vector of successive indexes in {1, . . . , p}. In addition they add the multivariate tests on the complementary sets of each interval, i.e., they do also test each hypothesis H The tests are performed by the Nonparametric Combination Procedure (NPC), see [25], that constructs multivariate permutation tests by means of combining univariate-synchronized permutation tests.
In the third step they obtain the adjusted p-value for the kth component λ

Perturbation strategy
Mimicking the approach proposed by [17] and [8], we restrict our perturbed networks to having the same numbers of vertices and edges as the original unperturbed network, hence only the positions of the edges change. In other words we apply a Degree Preserving Randomization. Our perturbation strategy relies on the rewire function belonging to the R package igraph, using the option keeping degseq.
Moreover, we expect that a network perturbed by only a small amount has just a few edges moved in different communities, while a maximally perturbed network produces completely random clusters. In [17] the perturbation strategy is achieved by removing each edge with a certain probability α and replacing it with another edge between a pair of vertex (i, j) chosen at random with a probability proportional to the degree of i and j. Our perturbation strategy consists in randomly rewiring a percentage p of edges while preserving the original graph's degree distribution. The rewiring algorithm indeed chooses two arbitrary edges in each step (e.g. (a, b) and (c, d)) and substitutes them with (a, d) and (c, b), if they do not already exists in the graph. The algorithm does not create multiple edges.
A null percentage of permutation p = 0 corresponds to the original unperturbed graph, while p = 1 corresponds to the maximal perturbation level. Varying the percentage p from 0 (original graph) to 1 (maximal perturbation), many perturbed graph are generated and compared to the partition on the original graph by means of V I. Indeed we generated 10 perturbed graph for each different level of p ∈ [0, 1]. Then, from each of the obtained graphs, we generated other 10 graphs rewiring 1% of edges each time. Hence resulting in 100 graphs for each level of p ∈ [0, 1]. In our setting we chose 20 levels of p.

Results
The overall procedure proposed in the present paper was implemented in R and validated both on simulated and real networks as will be described in the following Sections 3.1 and 3.2. In Figure 1 we provide a flow chart that summarises our procedure. For each of the analysed networks (either simulated or real) we performed the first step of the overall procedure described in Section 2, using some tools embedded in the R package igraph. We chose igraph because it provides an implementation of graph algorithms able to fast identifying community structures in large graphs. In particular we used two community extraction functions, one based on a greedy optimisation of the modularity (cluster f ast greedy) and another based on a multi-level optimisation of the modularity (cluster louvain). More Figure 1: Overall procedure map specifically cluster f ast greedy implements the hierarchical agglomeration algorithm for detecting community structure described in [6] and cluster louvain is based on the hierarchical approach proposed in [5]. Both these methods enables for an automatic definition of the optimal number of communities, are specific for large networks and are based on the optimisation of the modularity. They are briefly summarised in the following.
As regards the testing methodologies we used the bioconductor package gprege available at https://www.bioconductor.org/packages/release/bioc/html/gprege.html for the GP regression, the R code from the professor Staicu's web-site http://www4.stat.ncsu.edu/~staicu/ for the Functional Principal Component test and the R package fdatest available at https://cran.r-project.org/web/packages/fdatest/index.html for the Interval-wise Functional test.

Fast greedy method
F ast greedy is the modularity optimization algorithm introduced by Clauset, Newman and Moore [6]. This method is essentially a fast implementation of a previous technique proposed by Newman [24]. Starting from a set of isolated nodes, the links of the original graph are iteratively added such to produce the largest possible increase of the modularity. Adding a first edge to the set of disconnected vertices reduces the number of groups forming a new partition of the graph. The edge is chosen such that this partition gives the maximum increase (minimum decrease) of modularity with respect to the previous configuration. All other edges are added based on the same principle. At each iteration step, the variation of modularity given by the merger of any two communities of the running partition is computed and the best merger chosen. The fast version of Clauset, Newman and Moore, which uses more efficient data structures, has a complexity of O(N log 2 N ) on sparse graphs.

Louvain method
Louvain method is the fast modularity optimization by Blondel et al. [5]. This technique consists of two steps, executed alternatively. Initially, each node is in its own community. In step 1, nodes are considered one by one, and each one is placed in the neighbouring community (including its own) that maximizes the modularity gain. This is repeated until no node is moved (the obtained decomposition provides therefore a local optimization of Newman-Girvan modularity). After a partition is identified in this way, in step 2 communities are replaced by super-nodes, yielding a smaller weighted network where two super-nodes are connected if there is at least an edge between vertices of the corresponding communities. The two steps of the algorithm are then repeated until modularity (which is always computed with respect to the original graph) does not increase any further.
As pointed out in [11], this method offers a fair compromise between the accuracy of the estimate of the modularity maximum, which is better than that delivered by greedy techniques like the one by Clauset et al. above, and computational complexity, which is essentially linear in the number of links of the graph.

Application to simulated data
In order to show the ability of our method to validate a network clustering, we applied it to modular random network graphs generated using the model implemented in [33]. The model generates undirected, simple, connected graphs with prescribed degree sequences and a specified level of community structure, while maintaining a graph structure that is otherwise as random (uncorrelated) as possible over a broad range of distributions of network degree and community size. The model in [33] is specified by the network size, the average network degree, the number of modules, the modularity, the degree distribution and the module size distribution. The generated graph results also to be as random as possible, to contain no self loops (edges connecting a node to itself), multi-edges (multiple edges between a pair of nodes), isolate nodes (nodes with no edges), or disconnected components (see [33] for details). Specifically, we generated a modular random graph for each level of modularity Q = 0, 0.2, 0.4, 0.6, 0.8, using a power law for degree distribution and for module size distribution, with size=2000, number of modules=10 and average degree=10. For each graph, the corresponding null model was generated using the Configuration Model.
The application of the overall procedure on the simulated datasets is summarised in Tables 1 and 2 and in Figures 2 and 3, respectively.

Gaussian Process results
The application of the Gaussian Processes approach described in Section 2.1 to the simulated networks is summarised in Table 1 The resulting BF show a growing trend from modularity Q = 0 to Q = 0.8 after clustering with either clustering f ast greedy or louvain. This gives strong statistical evidence that networks with a high modularity have a robust clustering structure (significantly different from the random). However, note that louvain method is able to recover a non random clustering also at low modularity (Q ≥ 0.2), the other way around f ast greedy needs a more defined structure (Q ≥ 0.4).

Functional Principal Component testing results
Similarly, Table 2   As we can see the False Discovery rate adjusted p-values decreases drastically when the simulated network modularity grows from Q = 0 to Q = 0.8, after clustering with either f ast greedy or louvain algorithms, also if louvain is able to recover a non random clustering at lower modularity (Q ≥ 0.2) than f ast greedy. This result agree with the previous one obtained by GP.

Interval-wise Functional Testing results
The application of the Interval-wise Testing procedure described in Section 2.3 to the simulated datasets after clustering via f ast greedy or louvain are depicted respectively in Figures 2 and 3. In each figure, panels (a), (c), (e), (g), (i) show the VI curves for the null model (V Ic random ) and for the actual model (V Ic). The two curves appear to be very close for low modularity values and depart from each other as the modularity increases till a maximum of Q = 1. In panels (b), (d), (f ), (h), (j) this is quantified locally by a specific adjusted p-value in each sub-interval. Significant p-values are falling under the horizontal red line corresponding to the critical value of 0.05. As we can see, either using louvain or f ast greedy as clustering methods yields to the similar IT P results conclusion. Of course when there is no perturbation (i.e. at level p = 0) the two curves are coincident and hence not significantly different. When Q ≥ 0.4 the two V I curves are significantly different at any perturbation level, apart from some cases at p ≥ 0.5 and Q = 0.8 where the V Ic is close to the V Ic random , indeed note that if we strongly perturb a network (p ≥ 0.5, i.e. we rewire more than 50% of edges) it approaches a random network. Also in this case Real Networks summary Facebook Nexus 5 Nexus 15 Barabasi  Nodes  4039  2617  4941  1870  Edges  88234  11855  6594  2240  Table 3: Nodes and Edges number for each of the four real analysed datasets louvain is able to recover a non random clustering at lower modularity (Q ≥ 0.2) than f ast greedy, confirming the results obtained by the other two approaches. Moreover note that at Q = 0, when applying louvain, the method detects some interval where there is a significant difference between the two curves but for strong perturbation level (p ≥ 0.5).

Application to real data
In order to provide an example of our analysis work-flow, we selected four different publicly available datasets namely two biological (protein-protein interaction networks, N exus 5 and Barabasi), one representing the western states power grid of United States (N exus 15) and a social dataset (F acebook). Note that N exus is an online repository of networks, with an AP I that allow programmatic queries against it, and programmatic data download as well. These functions can be used to query it and download data from it, directly as an igraph graph. The total number of nodes and edges of these four real networks are summarised in Table 3 and displayed in Figure 4.

Nexus 5
This dataset consists of an undirected protein-protein interaction network in yeast. This data set was compiled by von Mering et al. (see [35]) combining various sources. Only the interactions that have high and medium confidence are included here.

Barabasi
This dataset consists of the protein-protein interaction network in Saccharomyces cerevisiae described and analysed in [1]. It is derived from combined, non-overlapping data, obtained mostly by systematic two-hybrid analyses. Data are available at http://www3.nd.edu/ networks/resources.htm.

Nexus 15
This dataset is an undirected, unweighted network representing the topology of the Western States Power Grid of the United States and was compiled by Duncan Watts and Steven Strogatz. Data are available at http://cdg.columbia.edu/cdg/datasets, [36].

Facebook
This dataset consists of 'circles' (or 'friends lists') from Facebook [22]. The authors obtained profile and network data from 10 ego-networks, consisting of 193 circles and 4,039 users. To do so they developed their own Facebook application and conducted a survey of ten users, who were asked to manually identify all the circles to which their friends belonged. On average, users identified 19 circles in their ego-networks, with an average circle size of 22 friends. Examples of such circles include students of common universities, sports teams, relatives, etc. Data are available at http://snap.stanford.edu/data/egonets-Facebook.html.
The application of the overall procedure on the just described real datasets is summarised in Tables 4 and 5 and in Figures 5 and 6, respectively.

Gaussian Process results
The application of the Gaussian Processes approach described in Section 2.1 to the four real networks is summarised in Table 4. The resulting BF are very high for either f ast greedy or louvain clustering. This gives strong statistical evidence that the four analysed networks have a robust clustering structure hence the recovered community structures are not likely to be random.

Interval-wise Functional Testing results
The application of the Interval-wise Testing procedure described in Section 2.3 to the real datasets after clustering via f ast greedy or louvain are depicted respectively in Figures 5  and 6. In each figure, panels (a), (c), (e), (g) show the VI curves for the null model (V Ic random ) and for the actual model (V Ic). In all the cases the two curves appear to be very close for high perturbation values and depart from each other as perturbation level approaches zero. In panels (b), (d), (f ), (h) this is quantified locally by a specific adjusted p-value in each sub-interval. Also in this case significant p-values are falling under the horizontal red line corresponding to the critical value of 0.05. As expected, either using louvein or f ast greedy as clustering methods yields to similar results conclusion. As already observed for the synthetic datasets, if we strongly perturb a network (p ≥ 0.5, i.e. we rewire more than 50% of edges) it approaches a random network, indeed the two V I curves become very close, and the p-value could survive the threshold.

Conclusions and Discussions
In this paper we propose an effective procedure to evaluate the robustness of a clustering. Given a community detection method and a network of interest, our methodology enables to clearly detect if the community structure found by some algorithms is statistically significant or is a result of chance, permitting to examine the stability of the partition recovered. As suggested in [17], we specify a perturbation strategy and a null model to build a set of procedures based on VI as stability measure. This enables to build the VI curve as a function of the perturbation percentage and to compare it with the corresponding null model curve in the functional data analysis framework.
We point out that our methodology could also be used to compare different clustering methodologies, indeed given two clusterings on the same network, we could test the agreement between the two recovered partitions via the direct comparison of the corresponding VI curves as defined by our procedure in Section 2. For example, all the three procedures we used point out that louvain method is able to recover a non random clustering also at low modularity (Q ≥ 0.2), the other way around f ast greedy needs a more defined structure (Q ≥ 0.4).
However, it is out of the scope of the present paper the comparative evaluation of different community extraction methods. The two methodologies louvain and f ast greedy were indeed only instrumental to the exemplification of our procedure. Both of them were selected at this stage as they both enables for an automatic definition of the optimal number of communities and are based on the optimisation of the modularity, that plays a key role in describing community structures.
An interesting and straightforward extension of the current paper would be using a different clustering stability measure, for example the N ormalized M utual Inf ormation measure proposed in [9]. This would also lead to a comparison of the performance of different measures for community structure comparison.