A graph is a pair of sets G = (V, E) where V is a set of n vertices and E is a set of m edges that connect two vertices of V. A random graph g is a family of graphs, where the members of the family are generated by some probability law.

An undirected graph G with n vertices can be represented by its n × n adjacency matrix A where A_ij = A_ji = 1 if vertices i and j are connected, and 0 otherwise. The spectrum of G is the set of eigenvalues (λ₁ ≥ λ₂ ≥ … ≥ λ_n) of the adjacency matrix A. For an undirected graph, the adjacency matrix is symmetric, and thus, its eigenvalues are real (Strang, 2011).

Given a set of random graphs g generated by the same probability law, the set of eigenvalues Λ are random vectors. Let δ be the Dirac delta function and the brackets “<>” indicate the expectation with respect to the probability law of the random graph, the spectral distribution of a random graph g is defined as:

The spectral distribution is directly associated with the structural features of the graphs (Albert and Barabási, 2002) and can be considered as a fingerprint of the random graph, where each random graph model is associated with a specific spectral distribution ρ_g (Van Mieghem, 2010).

Once the spectral distribution of a graph is defined, we can describe a measure of similarity between two spectral distributions. If two spectral distributions are different, then the respective graphs should be different.

Let ρ_g1 and ρ_g2 be the spectral distributions of random graphs g₁ and g₂, respectively. If the support of ρ_g2 contains the support of ρ_g1, the Kullback-Leibler (KL) divergence between two spectral distributions ρ_g1 and ρ_g2 is (Kullback and Leibler, 1951): (2)KL(ρg1|ρg2)=∫-∞+∞ρg1(λ)logρg1ρg2dλ, otherwise, KL(ρ_g1|ρ_g2) = +∞ (we assume 0 log00=0).

For Equation (2), ρ_g2 is called the reference measure. The KL divergence is non-negative and zero if and only if ρ_g1 and ρ_g2 are equal.

To estimate the spectral density (ρ^g), we use the same procedure described by Takahashi et al. (2012). First we compute the eigenvalues of the adjacency matrix of the graph. Then, we apply a Gaussian kernel regression using the Nadaraya-Watson estimator (Nadaraya, 1964) for regularization of the estimator. Finally, we normalize the density to obtain the integral below the curve equal to one. In general, smaller sample and/or graph sizes require larger bandwidth and smaller bin numbers for the smoothing kernels. The opposite holds for larger sample and/or graph sizes. The exact bandwidth size and bin number that maximize the statistical power depends on the data and the alternative hypotheses, but some rules of thumb exist in the literature that have shown good performance in our simulations. The bandwidth of the kernel is chosen as λ1 - λnnumber of bins (Sain, 1996), where the number of bins is selected by using the Sturge's criterion (Sturges, 1926). Type I errors are controlled by our bootstrap procedure discussed below for any choice of bandwidth size and bin number.

We are now able to describe ANOGVA. Given k populations of graphs g₁, g₂, …, g_k, the test consists of verifying if all populations of graphs were generated by the same random graph model. For this, we test if all the spectral distributions are equal.

Let ρ^g1,ρ^g2,…,ρ^gk be the estimated spectral distributions of populations of graphs g₁, g₂, …, g_k, respectively, where ρ^gi (i = 1, …, k) is the average of the graphs spectra in population g_i. Also, set ρ^gM=∑i=1kρ^gik. The support of ρ^gM includes the support of ρ^gi for any i. Formally, we test:

H₀: KL(ρ_g1, ρ_gM) = KL(ρ_g2, ρ_gM) = … = KL(ρ_gk, ρ_gM) = 0, i.e., the graphs from g₁, g₂, …, g_k are generated by the same random graph model (the spectral distributions are equal).

H₁: “At least one population of graphs is generated in a different manner”.

We will use the statistic Δ=∑i=1kKL(ρ^gi,ρ^gM) to build the test statistic. This statistic is the generalization of the Jensen-Shannon divergence (Jensen, 1906; Shannon, 1948) for k > 2. Under the null hypothesis, we expect small Δ, while large Δ suggests a rejection of the null hypothesis. The exact or asymptotic distribution of Δ under the null hypothesis is not known; therefore, we use a computational procedure based on the permutation test to construct the empirical distribution. The steps for the permutation test is as follows:

Construct permuted samples gi*, for i = 1, …, k by resampling (without replacement) |g_i| graphs from the entire dataset {g₁∪g₂∪ … ∪g_k}.

Figure 1 illustrates the idea behind ANOGVA. In summary, the spectral distribution of each population is compared to the reference distribution (the average of the spectral distributions, ρ_gM). If the sum of the distances (KL divergence) is large, it means that at least one of the spectral distributions is different when compared to the reference (ρ_gM). In other words, at least one of the populations of graphs was generated by a different random graph model and/or set of parameters.

For the simulation studies, we analyzed five random graph models, namely Erdös-Rényi (Erdös and Rényi, 1960), geometric (Penrose, 2003), k-regular (Meringer, 1999), Watts-Strogatz (Watts and Strogatz, 1998), and Barabási-Albert (Barabasi and Albert, 1999). All simulations were carried out in R using the package igraph. The name of the functions and respective parameters analyzed in our study are: function erdos.renyi.game - parameter p; function grg.game - parameter radius; function k.regular.game - parameter k; function watts.strogatz.game - parameter p; and function barabasi.game - parameter power. For the Watts-Strogatz random graph model, we selected the parameter p as the varying parameter in our simulations because it changes the graph structure without altering the number of edges.

We evaluated the controls of types I and II errors of ANOGVA in verifying whether graphs are generated by the same random graph model.

Figures 2A–D illustrate, respectively, the ROC curves obtained by simulating scenarios 1, 2, 3, and 4 described in Section 2.6.1. The x-axis represents the p-value threshold and the y-axis represents the proportion of rejected null hypothesis given a p-value threshold. The ROC curve under the null hypothesis lies in the diagonal. Under the alternative hypothesis, we expect to obtain a curve above the diagonal. In our case, the nominal p-value is on the x-axis and the proportion of rejected null hypotheses is on the y-axis.

The further is the curve above the diagonal, the higher is the power of the test. By analyzing the ROC curves, it is possible to notice that (i) Figure 2A—the ROC curves under the null hypothesis lie in the diagonal as expected, i.e., the statistical test is effectively controlling the rate of false positives (the proportion of rejected null hypothesis is as expected by the p-value threshold); (ii) Figure 2B—the power of the test increases as the number of graphs increases under the alternative hypothesis (the parameter of one of the populations is generated in a different manner); (iii) Figure 2C—the power of the test increases as the datasets are more balanced; and (iv) Figure 2D—the case the graph samples are generated by different models, the power of the test increases as the number of misclassified graphs decreases. When the number of misclassification is zero, it is possible to notice that the power of the method is higher than when only the parameter is different (Figure 2C).

We compared ANOGVA with the application of ANOVA on other features, such as the number of edges, betweenness centrality, closeness centrality, assortativity, and transitivity in five random graph models namely, Erdös-Rényi, geometric, k-regular, Watts-Strogatz, and Barabási-Albert. Figure 3 describes the proportion of rejected null hypotheses in these graph models. By analyzing Figure 3, we notice that only ANOGVA followed by transitivity are able to discriminate all random graph models with different parameters, including the Watts-Strogatz random graph model. Betweenness centrality presented low power for geometric, k-regular, Watts-Strogatz, and Barabási-Albert random graph models. Closeness centrality was not able to identify differences between Watts-Strogatz random graph models. Assortativity coefficient presented low power for Erdös-Rényi and Watts-Strogatz random graph models.

Functional brain networks were constructed as described in Section 2.7.2. The five sub-networks depicted in Figure 4 are based on sub-networks defined by Sato et al. (2015), namely somatomotor, visual, default-mode, cerebellar, and fronto-parietal.

Here, we focused on the identification of which sub-network is associated with autism and Asperger by using ANOGVA. The number of permutations was set to 1,000. First, we compared the three groups (controls vs. autism vs. Asperger) for each sub-network to verify if there is at least one population that differs from the others. The test indicated a significant difference at a p-value threshold of 0.05 for cerebellar (p = 0.04), and not for somatomotor (p = 0.13), visual (p = 0.95), default-mode (p = 0.17), and fronto-parietal (p = 0.21). These results suggest that the structure of the cerebellar functional sub-network is different at least in one of the populations among controls, autistic and Asperger patients. To identify which population is not equally generated in the cerebellar cluster, we carried out pairwise comparisons among the groups. Results indicate that there is no statistical evidence to discriminate controls vs. Asperger (p = 0.99) and autism vs. Asperger (p = 0.13), but there is significant difference between controls and autism (p < 0.001). This result indicates that the random process underlying the group control and autism are different.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.