# A spectral-based framework for hypothesis testing in populations of networks

Li Chen1 ,  Nathaniel Josephs2 ,  Lizhen Lin3 ,  Jie Zhou1 and  Eric D.
Kolaczyk2 College of Mathematics, Sichuan University, China  Department of
Mathematics and Statistics, Boston University, USA  Department of Applied and
Computational Mathematics and Statistics, The University of Notre Dame, USA

###### Abstract.

In this paper, we propose a new spectral-based approach to hypothesis testing
for populations of networks. The primary goal is to develop a test to
determine whether two given samples of networks come from the same random
model or distribution. Our test statistic is based on the trace of the third
order for a centered and scaled adjacency matrix, which we prove converges to
the standard normal distribution as the number of nodes tends to infinity. The
asymptotic power guarantee of the test is also provided. The proper interplay
between the number of networks and the number of nodes for each network is
explored in characterizing the theoretical properties of the proposed testing
statistics. Our tests are applicable to both binary and weighted networks,
operate under a very general framework where the networks are allowed to be
large and sparse, and can be extended to multiple-sample testing. We provide
an extensive simulation study to demonstrate the superior performance of our
test over existing methods and apply our test to three real datasets.

##  1\. Introduction

In this work, we consider an inference problem related to populations of
networks in which each sample or data point is a network. The statistical
network analysis literature has been largely focused on proposing models and
algorithms for analyzing a single network. However, the increasing prevalence
of multiple network datasets, in which the network is the fundamental data
object, along with the need to extract useful scientific information from
them, have motivated the demand for developing statistical methods of
inference for populations of networks.

For example, in brain network data (see the COBRE data of Section 5), one may
be interested in testing whether the brain network structure from a group of
individuals with schizophrenia is different from that of a group of healthy
controls. Given a collection or sample of such networks, one might also be
interested in estimating some mean network feature, which could provide a
notion either of averaging networks or of clustering networks into different
groups [23]. All of these cases are inference tasks for one or two samples of
network objects, both of which have been recently explored in the literature.

In [15], the authors consider two-sample testing for networks with
applications to functional neuroimaging. Very recently, this work was extended
in [18] through a geometric and statistical framework for inference on
populations of unlabeled networks by providing a geometric characterization of
the space of unlabeled networks and deriving a central limit theorem for the
sample Fréchet mean. Supervised and unsupervised learning such as clustering,
regression, and classification for network objects have also been considered
in the literature. See, e.g., [2] and [17], with the former considering
network classification in neuroimaging and the latter employing Bayesian
methods for classification, anomaly detection, and survival analysis.

Herein, we focus on the problem of two-sample hypothesis testing for
populations of networks. Compared to much of the work in the literature, such
as [15], in which the number of nodes is fixed, we consider a general
framework which allows both the number of nodes and the sample size (the
number of networks) to grow. In another related work, [13] study two-sample
problems from a minimax testing perspective on testing whether two samples of
binary networks of $n$ nodes are generated from the same link probability
matrix against an alternative that says the two link probability matrices are
$\rho$ apart with respect to some matrix norm. Their work focuses on the
theoretical characterization of minimax separation with respect to the number
of networks $m$, the number of nodes $n$, and different matrix norms. [28]
study whether two random dot product networks ($m=1$) defined on different
vertex sets are generated from the model or not.

Our test statistics are spectral based and not restricted to a given network
structure. We utilize the trace of the third order for a centered and scaled
adjacency matrix, which is proven to converge to the standard normal
distribution as the number of nodes tends to infinity. In addition, we show
that the asymptotic power tends to one as the number of nodes increases. Since
we also want to understand the limiting behavior as the sample size increases,
we explore the proper interplay between the asymptotics in the number of
networks and in the number of nodes for each network when characterizing the
theoretical properties of our proposed testing statistics. These statistics
are conceptually simple and computational friendly and we provide an extensive
simulation study under various models to demonstrate the superior performance
of our test over existing methods. In almost all the cases examined in our
study, the proposed test statistics achieve the nominal rejection rate under
the null and a power close to one under the alternative. We also apply our
test to three real datasets, based on both weighted and binary networks.

The idea of applying a spectral method based on random matrix theory to
network data is a natural one, as the network data (e.g., the adjacency or
Laplacian matrix) can be naturally viewed as a random matrix. Spectral-based
hypothesis tests, in particular, have been proposed in [4] and [9] for testing
the community structure and determining the number of clusters within a single
network. A spectral-based test based on a Tracy-Widom law for hypothesis
testing of populations of networks and change point detection in networks can
be found in [7] and [6]. Compared to those two works, our spectral-based test
has asymptotic standard normal distribution with much faster convergence rate
under the null compared to the slow convergence of a test that has a Tracy-
Widom law. Furthermore, our testing statistics require much milder conditions
for the theoretical performance guarantees: an error estimate of the link
probability estimates with $o_{p}(1)$ is needed in comparing to an error
condition of $o_{p}(n^{-2/3})$ required in [6]. Note that although the
statistic in [9] also has an asymptotically normal law, it is limited to
testing the presence of community structure in a single network versus the
null Erdös-Rényi model, whereas our statistic can test the difference between
arbitrary network models and can be applied to either binary or weighted
networks in both two-sample and multiple-sample frameworks.

The remainder of the paper is organized as follows. In Section 2, we describe
our proposed spectral-based testing statistics and derive their asymptotic
null distributions as well as asymptotic power results. We extend our test for
weighted networks and multiple-sample testing in Section 3. Results of
extensive simulation studies are reported in Section 4 and analysis on three
real network datasets are given in Section 5. We conclude in Section 6 with a
few final remarks and possible future directions of this work.

##  2\. A new spectral-based test for binary networks

In this section, we first propose a new spectral-based test for testing the
difference between distributions of two samples of binary networks.
Specifically, we consider two samples of networks on the same $n$ nodes with
possibly different sample sizes $m_{1}$ and $m_{2}$. We assume one observes
the independent and identically distributed symmetric binary adjacency
matrices $A_{1}^{(1)},\ldots,A_{1}^{(m_{1})}$, with conditionally independent
entries generated from a symmetric link probability matrix $P_{1}$, i.e.

$A_{1,ij}^{(k)}\sim\text{Bernoulli}(P_{1,ij})\enskip,$

for $k=1,2,\ldots,m_{1}$, $i,j=1,2,\ldots,n$. Similarly, one observes a second
sample of adjacency matrices $A_{2}^{(1)},\ldots,A_{2}^{(m_{2})}$ with

$A_{2,ij}^{(k)}\sim\text{Bernoulli}(P_{2,ij})\enskip,$

generated from the same model with link probability matrix $P_{2}$. Our goal
is to test whether the two samples of networks have the same graph structure
or not, which is equivalent to testing

(2.1) $H_{0}:P_{1}=P_{2}~{}\text{against}~{}H_{1}:P_{1}\neq P_{2}\enskip.$

To address this, we propose a new statistic that utilizes results from random
matrix theory. For necessary background on spectral properties of
inhomogeneous networks, which are used heavily in this work, see Appendix A.

###  2.1. New spectral test for binary networks

Given two samples of networks $\\{A_{1}^{(k)}\\}_{k=1}^{m_{1}}$ and
$\\{A_{2}^{(k)}\\}_{k=1}^{m_{2}}$ sampled from the link probability matrices
$P_{1}$ and $P_{2}$, respectively, we introduce the normalized matrix with
elements as follows:

(2.2)
$Z_{ij}=\begin{cases}\frac{\bar{A}_{1,ij}-\bar{A}_{2,ij}}{\sqrt{n\big{(}\frac{1}{m_{1}}P_{1,ij}(1-P_{1,ij})+\frac{1}{m_{2}}P_{2,ij}(1-P_{2,ij})\big{)}}}&\text{
if }i\neq j\\\ B_{ij}&\text{ if }i=j\end{cases}\enskip,$

where $\bar{A}_{u}$ is the sample average of adjacency matrices in the $u$th
group, for $u=1,2$,

(2.3) $\bar{A}_{u}=\frac{1}{m_{u}}\sum_{k=1}^{m_{u}}A_{u}^{(k)}\enskip,$

and $B$ is an $n\times n$ diagonal matrix with $B_{ii}$ given by i.i.d random
variables such that

(2.4) $P(B_{ii}=-1/\sqrt{n})=P(B_{ii}=1/\sqrt{n})=1/2\enskip,$

for $i=1,\ldots,n$.

Consider the test statistic

(2.5) $\theta=\frac{1}{\sqrt{15}}\text{Tr}(Z^{3})\enskip,$

where $\text{Tr}(\cdot)$ represents the trace operator. We have the following
theorem on the asymptotic distribution of $\theta$ under the null hypothesis.

######  Theorem 1.

Let $Z$ be given as in (2.2). Assume the sample sizes
$m_{1}=O(n^{\alpha_{1}})$ and $m_{2}=O(n^{\alpha_{2}})$ for some
$\alpha_{1},\alpha_{2}\in(0,1)$. Then, under the null hypothesis
$P_{1}=P_{2}$, for the scaled test statistic
$\theta=\frac{1}{\sqrt{15}}\text{Tr}(Z^{3})$, we have

(2.6) $\theta\overset{d}{\to}\mathcal{N}(0,1)\quad\text{as}\quad
n\to\infty\enskip.$

We defer the details of the proof to Appendix B.1. However, the overview of
the argument is as follows. First, one can see that under the null hypothesis
of $P_{1}=P_{2}$, $Z$ is a Wigner matrix satisfying $\mathrm{E}(Z_{ij})=0$ and
$\text{Var}(Z_{ij})=1/n$. Then the remainder of the proof proceeds in three
steps. We begin by showing that the empirical spectral distribution of $Z$
converges to the semicircular law almost surely. Next, we verify that
$X=\sqrt{n}Z$ satisfies conditions $(1)$–$(3)$ of Lemma 2, after which the
asymptotic normality of $\theta$ follows. Finally, the mean and the variance
are obtained from [9].

To formalize a testing framework using $\theta$ in (2.5), we need to account
for the fact that the diagonal matrix $B$ in (2.4) is random. We do so by
employing a Monte Carlo procedure, which we describe in Algorithm 1. Our
output is an empirical confidence level, which is the rejection rate based on
the test statistics computed from the Monte Carlo samples of $B$.

New Spectral-Based Hypothesis Test
$\big{(}\\{A_{1}^{(k)}\\}_{k=1}^{m_{1}},\\{A_{2}^{(k)}\\}_{k=1}^{m_{2}},\alpha,Q\big{)}$;

Input : Adjacency matrices $\\{A_{1}^{(k)}\\}_{k=1}^{m_{1}}$ and
$\\{A_{2}^{(k)}\\}_{k=1}^{m_{2}}$ for groups 1 and 2

Significance level $\alpha$

Number of Monte Carlo samples $Q$

Output : Empirical significance level rej_rate

Compute $\bar{A}_{u}$ for $u=1,2$ using (2.3) ;

for  _$q=1,\ldots,Q$_ do in parallel

Sample $B^{(q)}$ satisfying (2.4) ;

Compute $Z^{(q)}$ in (2.2) using $B^{(q)}$ ;

Compute $\theta^{(q)}$ in (2.5) using $Z^{(q)}$ ;

end forpar

$\text{rej\\_rate}=\frac{1}{Q}\sum_{q=1}^{Q}I\big{(}|\theta^{(q)}|>\mu_{\alpha/2}\big{)}$

Algorithm 1 Procedure for testing using the statistic in (2.5). The output is
an empirical significance level based on Monte Carlo test statistics, where
$I(\cdot)$ is an indicator function and $\mu_{\alpha/2}$ is the $\alpha/2$
upper quantile of $\mathcal{N}(0,1)$.

######  Remark 1.

In Algorithm 1, we deliberately do not output a p-value. For $Q=1$, we could
obtain a p-value using $2P\big{(}\theta>|\theta_{obs}^{(Q=1)}|)$ as in [4] and
[9], where $\theta_{obs}^{(Q=1)}$ is the sample test statistic and $\theta$
follows the null distribution of the testing statistic. In this case, though,
the p-value is implicitly conditional on $B$ and the authors’ simulations
reveal that the randomness of $B$ leads to highly variable p-values. Instead,
for our test, we propose computing many $\theta_{obs}^{(q)}$ in parallel to
reduce the noise induced by $B$. The analogous p-value estimate combining
these Monte Carlo test statistics would be
$2P\big{(}\theta_{Q}>|\bar{\theta}_{obs}|\big{)}$, where
$\bar{\theta}_{obs}=\frac{1}{Q}\sum_{q=1}^{Q}\theta_{obs}^{(q)}$.

######  Remark 2.

The rejection rate from our Monte Carlo estimator has the property that its
expectation under the null is the nominal significance level:

$\mathrm{E}\Big{(}\frac{1}{Q}\sum_{q=1}^{Q}I\big{(}|\theta^{(q)}|>\mu_{\alpha/2}\big{)}\Big{)}=P\big{(}|\theta^{(q)}|>\mu_{\alpha/2}\big{)}=\alpha\enskip.$

###  2.2. Test statistic based on estimated link probability matrices

Theorem 1 assumes that the true link probability matrices $P_{1}$ and $P_{2}$
are known, which is not the case in practice. Therefore, $\theta$ cannot be
used directly as a test statistic. A natural alternative is to plug in some
appropriate estimates of $P_{1}$ and $P_{2}$ with the hope that the plug-in
estimator for the test statistic retains asymptotic normality.

We denote the plug-in estimates of $P_{1}$ and $P_{2}$ by $\hat{P}_{1}$ and
$\hat{P}_{2}$, respectively. Then the empirical version of the normalized
matrix $Z$ in (2.2) can be written as

(2.7)
$\hat{Z}_{ij}=\begin{cases}\frac{\bar{A}_{1,ij}-\bar{A}_{2,ij}}{\sqrt{n\left(\frac{1}{m_{1}}\hat{P}_{1,ij}(1-\hat{P}_{1,ij})+\frac{1}{m_{2}}\hat{P}_{2,ij}(1-\hat{P}_{2,ij})\right)}}&\text{
if }i\neq j\\\ B_{ij}&\text{ if }i=j\end{cases}\enskip.$

The resulting test statistic is thus

(2.8) $\hat{\theta}=\frac{1}{\sqrt{15}}\text{Tr}{(\hat{Z}^{3})}\enskip,$

which has the following limiting law.

######  Theorem 2.

Under the two-sample framework of binary networks, let $\hat{Z}$ be given in
(2.7). As before, assume the sample sizes $m_{1}=O(n^{\alpha_{1}})$ and
$m_{2}=O(n^{\alpha_{2}})$ for some $\alpha_{1},\alpha_{2}\in(0,1)$. Suppose
$\hat{P}_{1}$ and $\hat{P}_{2}$ are some estimates of $P_{1}$ and $P_{2}$,
respectively. If $\max_{ij}|\hat{P}_{u,ij}-P_{u,ij}|=o_{p}(1)$ for $u=1,2$,
then, under the null hypothesis $P_{1}=P_{2}$, we have the following
asymptotic distribution of the scaled test statistic
$\hat{\theta}=\frac{1}{\sqrt{15}}\text{Tr}(\hat{Z}^{3})$:

$\hat{\theta}\overset{d}{\to}\mathcal{N}(0,1)\quad\text{as}\quad
n\to\infty\enskip.$

Again, we defer the proof to Appendix B.2, which relies on rewriting

$\text{Tr}(\hat{Z})^{3}=\text{Tr}(Z^{3})+3\text{Tr}\big{(}Z^{2}(Z\circ
H)\big{)}+3\text{Tr}\big{(}Z(Z\circ H)^{2}\big{)}+\text{Tr}\big{(}(Z\circ
H)^{3}\big{)}\enskip,$

where $\circ$ denotes the Hadamard product and $H$ is an $n\times n$ matrix
with entries $H_{ij}=\max_{u=1,2}O(\hat{P}_{u,ij}-P_{ij})=o_{p}(1)$. Each term
in the right of this equality can be written as an element-wise sum and can be
divided into two parts in terms of the corresponding subscripts, which we
bound to show convergence as $n$ goes to infinity.

###  2.3. Estimating link probability matrices

Theorem 2 requires the sample sizes of observed networks $m_{u}$ to grow with
$n$ at a rate of $n^{\alpha}$ for any $\alpha\in(0,1)$ and
$\max_{ij}|\hat{P}_{u,ij}-P_{u,ij}|=o_{p}(1)$, which are very mild conditions
satisfied by many estimation methods.

The simplest estimator of $P_{u,ij}$ is the sample mean of all the $(i,j)$
elements in the adjacency matrices of group $u$, where $u=1,2$. We refer to
this spectral method based on simple averages as SPE-AVG. It is not difficult
to see that $\max_{i,j}|\hat{P}_{u,ij}-P_{u,ij}|=o_{p}(m_{u}^{-1/2}\log(n))$.
Intuitively, SPE-AVG requires sample sizes to be large enough to achieve good
performance. This is also confirmed empirically by our extensive simulation
study in which SPE-AVG typically yields inferior performance compared to the
next methods we present.

Another possible average estimator of $P_{u,ij}$ is based on the stochastic
block model (SBM). The key idea is to approximate any graph with an SBM,
which, for large networks, is reasonable by Szemerédi’s regularity lemma [20].
The membership vector of nodes can be obtained by community algorithms such as
the method proposed in [24]. After the membership vector has been estimated,
we can simply approximate $P_{u,ij}$ by the sample mean of all the entries in
the submatrix over all $A_{u}^{(k)}$, $k=1,2,\cdots,m_{u}$, restricted to the
corresponding block consisting of the communities of $i$ and $j$. We refer to
this test method based on SBM as SPE-SBM. Assuming the true community number
is $K_{u}$, then the variance error satisfies
$\max_{i,j}|\hat{P}_{u,ij}-P_{u,ij}|=o_{p}(K_{u}m_{u}^{-1/2}n^{-1}\log(n))$.
It can be seen that the rate of SPE-SBM is better than that of SPE-AVG as long
as $K_{u}<n$, which is very easy to be satisfied. However, the property may be
limited by the assumption that the network topologies follow an SBM structure.

Finally, we introduce a new method of estimation based on the modified
neighborhood smoothing (MNBS) proposed in [31]. The idea is to perform
neighborhood smoothing to the matrix $\bar{A}$, which is the weighted average
of $m$ networks and the smoothing procedure is applied to a shrunken
neighborhood size. This results in a better bias-variance tradeoff compared to
the neighborhood smoothing (NBS) method proposed in [30] leading to a better
estimate of the link probability matrix with a smaller error. When the sample
size $m_{u}$ is small compared with $n$, which is more common in practice, and
satisfies $(m_{u}\log n)^{1/2}<n^{1/2}$, then from Lemma 9.3 in [31], the size
of neighborhood is $O_{p}((n\log n/m_{u})^{1/2})$. Thus the estimation error
of the link probability is $|\hat{P}_{u,ij}-P_{u,ij}|=O_{p}((m_{u}n\log
n)^{-1/4})$. We refer to this test method based on MNBS as SPE-MNBS. We also
note that SPE-MNBS puts no structure conditions on the networks. Therefore, we
expect SPE-MNBS to be generally applicable.

###  2.4. Asymptotic power guarantee

Next we consider the power of the test based on $\hat{\theta}$, which we
summarize in the following theorem.

######  Theorem 3.

Consider the alternative model of $P_{1}\neq P_{2}$ under the assumptions of
Theorem 1. Let $Z^{\prime\prime}$ be an $n\times n$ matrix with zero diagonals
and, for any $i\neq j$,

(2.9) $\displaystyle Z^{\prime\prime}_{ij}$
$\displaystyle=\frac{P_{1,ij}-P_{2,ij}}{\sqrt{n\left(\frac{1}{m_{1}}P_{ij}(1-P_{ij})+\frac{1}{m_{2}}P_{ij}(1-P_{ij})\right)}}\enskip.$

Define the partition $\\{1,\cdots,n\\}^{3}=S_{a}\cup S_{b}$, where $(i,k,l)\in
S_{a}$ indicates that
$Z^{\prime\prime}_{ik}Z^{\prime\prime}_{kl}Z^{\prime\prime}_{li}\geq 0$, and
$(i,k,l)\in S_{b}$, that
$Z^{\prime\prime}_{ik}Z^{\prime\prime}_{kl}Z^{\prime\prime}_{li}<0$. Write
$|S_{a}|=an^{3}$ and $|S_{b}|=bn^{3}$, with $a,b\in[0,1]$ satisfying $a+b=1$.
If either of the following conditions are satisfied,

$\displaystyle(i)$ $\displaystyle\quad a\min_{(i,k,l)\in
S_{a}}({Z^{\prime\prime}}_{ik})^{3}+b\min_{(i,k,l)\in
S_{b}}({Z^{\prime\prime}}_{ik})^{3}>0\enskip,$ $\displaystyle(ii)$
$\displaystyle\quad a\max_{(i,k,l)\in
S_{a}}({Z^{\prime\prime}}_{ik})^{3}+b\max_{(i,k,l)\in
S_{b}}({Z^{\prime\prime}}_{ik})^{3}<0\enskip,$

then

$\lim_{n\to\infty}P(|\hat{\theta}|>\mu_{\alpha/2})=1,\ \alpha\in(0,1)\enskip.$

The details of the proof, which are similar to those in Theorem 2, are given
in Appendix B.3.

######  Remark 3.

Note that there is a slight abuse of notation in conditions $(i)$ and $(ii)$
where the minimum operator is taken over all pairs of indices among $(i,j,k).$
Conditions $(i)$ and $(ii)$ above characterize the minimum signal difference
between $P_{1}$ and $P_{2}$ required for Theorem 3 to hold, which implies that
the power is asymptotically one when either of the sets $S_{a}$ or $S_{b}$ is
large enough.

##  3\. Extending our test to other settings

In this section, we extend our test to be used for weighted networks, as well
as for multiple samples in a manner analogous to one-way analysis of variance
(ANOVA).

###  3.1. Extension to weighted networks

We now consider a more general framework that focuses on weighted networks.
Let $F_{1}=\\{F_{1,ij}\\}$ and $F_{2}=\\{F_{2,ij}\\}$ for $i,j=1,\ldots,n$ be
two sequences of distributions defined on bounded intervals and specified by
some parameters. Let
$A_{1}^{(1)},\ldots,A_{1}^{(m_{1})}\overset{i.i.d}{\sim}F_{1}$ and
$A_{2}^{(1)},\ldots,A_{2}^{(m_{2})}\overset{i.i.d}{\sim}F_{2}$ be symmetric
weighted adjacency matrices for networks that are undirected and without self-
loops, i.e. $\ A_{u,ii}^{(k)}=0$ for $u=1,2,\ i=1,\ldots,n$, and
$k=1,\ldots,m_{u}$. Let $\Sigma_{u}$ denote an $n\times n$ matrix in which the
$(i,j)$ element is the variance of $A_{u,ij}^{(k)}$ for $k=1,2,\ldots,m_{u}$.
Note that its diagonal elements are 0 since $A_{u,ii}^{(k)}=0$. Finally, let
$\hat{\Sigma}_{u,ij}$ be an estimate of ${\Sigma}_{u,ij}$.

Our approach for weighted networks is to replace $P_{u,ij}(1-{P}_{u,ij})$ in
(2.2) and $\hat{P}_{u,ij}(1-\hat{P}_{u,ij})$ in (2.7) by ${\Sigma}_{u,ij}$ and
$\hat{\Sigma}_{u,ij}$, respectively. Just as in Section 2.3, estimates
$\hat{\Sigma}_{u,ij}$ can also be obtained using various methods, which will
be discussed later. For simplicity, we use the same notation as in Section
2.3.

For the weighted case, the testing problem in (2.1) is equivalent to

(3.1) $H_{0}:F_{1}=F_{2}~{}\text{against}~{}H_{1}:F_{1}\neq F_{2}\enskip.$

We define the normalized matrix $Z$ as

(3.2)
$Z_{ij}=\begin{cases}\frac{\bar{A}_{1,ij}-\bar{A}_{2,ij}}{\sqrt{n\big{(}\frac{1}{m_{1}}\Sigma_{1,ij}+\frac{1}{m_{2}}\Sigma_{2,ij}\big{)}}}&\text{
if }i\neq j\\\ B_{ij}&\text{ if }i=j\end{cases}\enskip.$

Then the asymptotic distribution of
$\theta=\frac{1}{\sqrt{15}}\text{Tr}(Z^{3})$ follows a standard normal
distribution under the null hypothesis, as stated in the following theorem.

######  Theorem 4.

Under the two-sample framework of weighted networks, let $Z$ be given in
(3.2). Assume sample sizes $m_{1}=O(n^{\alpha_{1}})$ and
$m_{2}=O(n^{\alpha_{2}})$ for some $\alpha_{1},\alpha_{2}\in(0,1)$. Then,
under the null hypothesis $F_{1}=F_{2}$, for scaled test statistic
$\theta=\frac{1}{\sqrt{15}}\text{Tr}(Z^{3})$, we have

(3.3) $\theta\overset{d}{\to}\mathcal{N}(0,1)\quad\text{as}\quad
n\to\infty\enskip.$

The proof is omitted since it is similar to that of Theorem 1.

######  Remark 4.

Although the two-sample testing framework for binary networks is a special
case of that of (3.1), we discuss the two cases separately. In the binary
case, our test statistic is obtained by plugging in an estimate of the link
probability matrix $P$, while our test statistic for the weighted networks
requires a plug-in estimate of the variance of each edge weight. Hence the
estimation methods differ for these two cases.

For practical application, the variance matrices $\Sigma_{1}$ and $\Sigma_{2}$
need to be estimated, with some conditions assumed to ensure that the
asymptotic normality of the new test statistic still holds. For
$\hat{\Sigma}_{1}$ and $\hat{\Sigma}_{2}$, the plug-in estimates of
$\Sigma_{1}$ and $\Sigma_{2}$, respectively, the empirical normalized matrix
of $Z$ in (3.2) can be written with entries as

(3.4)
$\hat{Z}_{ij}=\begin{cases}\frac{\bar{A}_{1,ij}-\bar{A}_{2,ij}}{\sqrt{n\left(\frac{1}{m_{1}}\hat{\Sigma}_{1,ij}+\frac{1}{m_{2}}\hat{\Sigma}_{2,ij}\right)}}&\text{
if }i\neq j\\\ B_{ij}&\text{ if }i=j\end{cases}\enskip.$

Therefore, our test statistic is

(3.5) $\hat{\theta}=\frac{1}{\sqrt{15}}\text{Tr}{(\hat{Z}^{3})}\enskip.$

Then we have the following limiting law.

######  Theorem 5.

In the two-sample framework of weighted networks, let $\hat{Z}$ be given in
(3.4). Assume the sample sizes $m_{1}=O(n^{\alpha_{1}})$ and
$m_{2}=O(n^{\alpha_{2}})$ for some $\alpha_{1},\alpha_{2}\in(0,1)$. Suppose
$\hat{\Sigma}_{1}$ and $\hat{\Sigma}_{2}$ are some estimates of $\Sigma_{1}$
and $\Sigma_{2}$, respectively. If
$\max_{ij}|\hat{\Sigma}_{u,ij}-\Sigma_{u,ij}|=o_{p}(1)$, $u=1,2$, then under
the null hypothesis $F_{1}=F_{2}$, we have the following asymptotic
distribution of the scaled test statistic
$\hat{\theta}=\frac{1}{\sqrt{15}}\text{Tr}{(\hat{Z}^{3})}$:

$\hat{\theta}\overset{d}{\to}\mathcal{N}(0,1)\quad\text{as}\quad
n\to\infty\enskip.$

The proof is similar to that of Theorem 2, so we only include the key
differences in Appendix B.4, as the remainder of the proof can be completed
straightforwardly.

We consider two estimates of $\Sigma_{u,ij}$. The first is obtained simply as
the sample variance of each element over all adjacency matrices in the same
group. For convenience, we still refer to this method as SPE-AVG. Then we have

(3.6)
$\max_{i,j}|\hat{\Sigma}_{u,ij}-\Sigma_{u,ij}|=O_{p}(m_{u}^{-1/2})\enskip.$

The details of the proof can be found in Appendix B.5. The order of the error
is the same as the binary case which implies that SPE-AVG is suitable for
large sample-size cases.

The second estimate of $\Sigma_{u,ij}$ is obtained similarly as SPE-SBM for
unweighted networks: assume each network comes from an SBM, approximate the
community membership vector, and compute the sample covariance within each
community as the sample variance of the nodes corresponding to that community
block (rather than the sample mean). Again, we refer to this method as SPE-SBM
as in the binary case. With a similar argument as the proof in the Appendix
B.5, one has
$\max_{i,j}|\hat{\Sigma}_{u,ij}-\Sigma_{u,ij}|=O_{p}(K_{u}m_{u}^{-1/2}n^{-1})$.
Therefore, the error condition in Theorem 5 is satisfied as long as $K_{u}<n$,
which should hold for most cases.

The power of the test for weighted networks is presented in the following
theorem.

######  Theorem 6.

Under the assumptions of Theorem 4 and the alternative model $F_{1}\neq
F_{2}$, let $Z^{\prime\prime}$ be an $n\times n$ matrix with zero diagonals
and for any $i\neq j$,

$\displaystyle Z^{\prime\prime}_{ij}$
$\displaystyle=\frac{P_{1,ij}-P_{2,ij}}{\sqrt{n\left(\frac{1}{m_{1}}\Sigma_{1,ij}+\frac{1}{m_{2}}\Sigma_{2,ij}\right)}}\enskip,$

Define $S_{a}$ and $S_{b}$ as in Theorem 3 based on the above
$Z^{\prime\prime}$. If either of the following conditions are satisfied,

$\displaystyle(i)$ $\displaystyle\quad a\min_{(i,k,l)\in
S_{a}}({Z^{\prime\prime}}_{ik})^{3}+b\min_{(i,k,l)\in
S_{b}}({Z^{\prime\prime}}_{ik})^{3}>0\enskip,$ $\displaystyle(ii)$
$\displaystyle\quad a\max_{(i,k,l)\in
S_{a}}({Z^{\prime\prime}}_{ik})^{3}+b\max_{(i,k,l)\in
S_{b}}({Z^{\prime\prime}}_{ik})^{3}<0\enskip,$

then

$\lim_{n\to\infty}P(|\hat{\theta}|>\mu_{\alpha/2})=1,\ \alpha\in(0,1)\enskip.$

Again, the proof is omitted as it is similar to that of Theorem 3.

###  3.2. Extension to multiple-sample testing

Finally, we consider the case when $S$ groups are present and $S>2$. Assume
one observes the symmetric binary adjacency matrices
$A_{s}^{(1)},\ldots,A_{s}^{(m_{s})}$ that are generated from a symmetric link
probability matrix $P_{s}$, i.e.

$A_{s,ij}^{(k)}\sim\text{Bernoulli}(P_{s,ij})\enskip,$

for $s=1,\ldots,S$, $k=1,\ldots,m_{s}$, and $i,j=1,\ldots,n$. Our goal is to
test whether there are any differences in the distributions of the $S$ groups,
which is equivalent to testing

(3.7) $H_{0}:P_{1}=P_{2}=\cdots=P_{S}~{}\text{ against }~{}H_{1}:P_{s}\text{
are not all equal}\enskip.$

This is analogous to one-way ANOVA.

We define the pairwise normalized matrices with elements as follows:

(3.8)
$Z^{(s)}_{ij}=\begin{cases}\frac{\bar{A}_{s,ij}-\bar{A}_{ij}}{\sqrt{n\Big{(}\big{(}\frac{1}{m_{s}}-\frac{2}{m}\big{)}P_{s,ij}(1-P_{s,ij})+\frac{1}{m^{2}}\sum_{s=1}^{S}m_{s}P_{s,ij}(1-P_{s,ij})}\Big{)}}&\text{
if }i\neq j\\\ B_{ij}&\text{ if }i=j\end{cases}\enskip,$

where $\bar{A}_{s}$ is the sample average of adjacency matrices in group $s$
as in (2.3), $\bar{A}$ is the overall sample average of all the adjacency
matrices,

$\bar{A}=\frac{1}{m}\sum_{s=1}^{S}\sum_{k=1}^{m_{s}}A_{s}^{(k)}\enskip,$

$m$ is the total sample size,

$m=\sum_{s=1}^{S}m_{s}\enskip,$

and $B$ is defined as in (2.4).

If $\theta^{(s)}=\frac{1}{\sqrt{15}}\text{Tr}\big{(}(Z^{(s)})^{3}\big{)}$,
then, under the null distribution and appropriate conditions on $m_{s}$,
Theorem 2 gives

$\theta^{(s)}\overset{d}{\to}\mathcal{N}(0,1)\quad\text{as}\quad
n\to\infty\enskip,$

and it follows that

$\big{(}\theta^{(s)}\big{)}^{2}\overset{d}{\to}\chi^{2}(1)\quad\text{as}\quad
n\to\infty\enskip.$

Unfortunately, $\theta^{(1)},\ldots,\theta^{(S)}$ are not independent, so the
sum of their squares is not $\chi^{2}(S)$. However, it is shown in [10] that
the sum of dependent $\chi^{2}$ random variables can be approximated with a
gamma distribution. Therefore, we have

(3.9)
$\theta\equiv\sum_{s=1}^{S}\big{(}\theta^{(s)}\big{)}^{2}\overset{d}{\to}\Gamma\Big{(}\frac{S}{u},u\Big{)}\quad\text{as}\quad
n\to\infty\enskip,$

where the scale parameter $u$ is given by

$u=2\Big{(}1+\frac{2\sum_{q\neq r}^{S}\rho_{qr}}{S}\Big{)}\enskip,$

with $\rho_{qr}$ the pairwise correlation between
$\big{(}\theta^{(q)}\big{)}^{2}$ and $\big{(}\theta^{(r)}\big{)}^{2}$.

As before, the true link probability matrices $P_{s}$ are unknown and need to
be estimated. We can estimate each $\hat{P}_{s}$ as in Section 2.3, and then
plug in these estimates to $Z^{(s)}$ in (3.8). Furthermore, although the
pairwise correlations $\rho_{qr}$ are not analytically tractable, they can
easily be estimated using the Monte Carlo simulations in Algorithm 1, which
does not add to the computational complexity.

Moreover, with this setup, it is possible to follow the same development of
Theorem 2 in proving convergence of the plug-in estimator $\hat{\theta}$ that
uses the estimated link probability matrices and estimated pairwise
correlations. Similarly, (3.8) can be extended to weighted networks as in
Section 3.1.

##  4\. Simulation studies

In this section, we illustrate the performance of our proposed tests through
extensive simulation study. For binary networks, we evaluate three plug-in
estimators for the link probability matrices – AVG, SBM, and MNBS – and
compare the results to the test proposed in [12], which involves an estimated
distance between two network distributions based on the Frobenius measure for
binary networks that allows $n$ to go to infinity. We refer to these four
tests as $\text{SPE-AVG},\ \text{SPE-SBM},\ \text{SPE-MNBS}$, and DFRO,
respectively. We do not include the test proposed in [15] for comparison
because their results are asymptotic in the sample size with a fixed number of
nodes and the authors note that their test is expected to decrease power in
larger dimensions, i.e. with more nodes.

In our simulation studies, we evaluate the test performance by estimating the
power when the alternative is true, as well as the null rejection rate
(rejection rate under the null). We also vary the number of nodes,
$n\in\\{100,200,\ldots,1000\\}$, and the sample sizes,
$m_{1}=m_{2}=m\in\\{10,50\\}$. In each example, we set the significance level
to $\alpha=0.05$. We follow the procedure described in Algorithm 1 with $Q=1$
and report the empirical significance level as the average rejection rate on
5000 separate samples of networks from the underlying distributions. Note that
sampling new networks allows us to use $Q=1$, but the results are similar if
we use 5000 separate samples of networks with $Q>1$.

With this design, we consider four types of random graph model for sampling
binary and weighted networks, as well as one multiple-sample setting. The
results are detailed in the remainder of this section, but the conclusions are
as follows. Overall, it appears that SPE-MNBS is the most robust to different
network structures and sample sizes. If the networks are drawn from an SBM,
then unsurprisingly SPE-SBM is suitable. Throughout, we see that SPE-AVG shows
significant improvement as the sample size increases. Finally, all three plug-
in estimates of the link probability matrices yield superior results for our
test compared to DFRO.

###  4.1. Stochastic block model (SBM)

In the first example, we consider an SBM structure with a block matrix given
as

(4.1) $\displaystyle
P_{\text{SBM}}=\left[\begin{matrix}0.5+\varepsilon_{1}&0.25\\\
0.25&0.5\end{matrix}\right]\enskip,$

where $\varepsilon_{1}$ depends on our hypothesis. The membership of the $i$th
node is

$M(i)=I(1\leq i\leq\lfloor n/3\rfloor)+2I(\lfloor n/3\rfloor+1\leq i\leq
n)\enskip.$

The first group of networks, $\\{A_{1}^{(k)}\\}_{k=1}^{m_{1}}$, is generated
from $P_{\text{SBM}}$ with $\varepsilon_{1}=0$. In the null setting, the
second group of networks, $\\{A_{2}^{(k)}\\}_{k=1}^{m_{2}}$, is also generated
from $P_{\text{SBM}}$ with $\varepsilon_{1}=0$, whereas
$\varepsilon_{1}=1/(5\log m)$ in the alternative setting. The results are
shown in the first row of Figure 4.1.

To investigate the performance of the tests for sparser networks, the same
setting as above is considered, except now with $\varepsilon_{1}=2/(5\log m)$
and with the link probability matrix $P_{\text{SBM}}$ scaled by a factor
$\rho=1/4$. The corresponding results are shown in the second row of Figure
4.1.

Figure 4.1. Simulation results for testing networks with an SBM structure for
different network orders and sample sizes. The first and the second rows are
for dense and sparse networks, respectively. (a) Null rejection rate. (b)
Power under the alternative.

We see from the first row of Figure 4.1, where the networks are dense, that
SPE-SBM and SPE-MNBS are close to the nominal level $\alpha=0.05$ under
$H_{0}$ and both achieve good power under $H_{1}$. We also observe that SPE-
AVG is the most powerful under $H_{1}$, but its rejection rates are too high
under $H_{0}$ when $m=10$. However, this issue is mitigated when we increase
the sample size to $m=50$ even though this makes $\varepsilon_{1}$ smaller,
i.e. more similar underlying SBM structures. For DFRO, it has zero rejection
rate under $H_{0}$ and increases to unit power more slowly than our proposed
tests. Similar results hold in the sparser settings on the bottom row, though
we note that SPE-AVG and DFRO perform more comparably to SPE-SBM and SPE-MNBS.

###  4.2. Graphon

In the second example, we consider a graphon structure from [5] in which

$\displaystyle
f_{0}(v_{1},v_{2})=(v_{1}^{2}+v_{2}^{2}+v_{1}^{1/2}+v_{2}^{1/2})/4\enskip.$

Then $P_{u,ij}=f_{0}(\eta_{i},\eta_{j})$, with
$\eta_{i}\overset{i.i.d}{\sim}U(0,1)$ for $i=1,2,\ldots,n$. We generate
$\\{A_{1}^{(k)}\\}_{k=1}^{m_{1}}$ from probability matrix $P_{1}$ according to
$f_{0}$. For the second group of networks, under the null hypothesis, we again
sample from $f_{0}$ to generate $\\{A_{2}^{(k)}\\}_{k=1}^{m_{2}}$. Under the
alternative hypothesis, we first randomly choose a subset
$S\subset\\{1,2,\ldots,n\\}$ with $|S|=\lfloor n/10\rfloor$, where
$\lfloor\cdot\rfloor$ is the floor operator, then generate
$\\{A_{2}^{(k)}\\}_{k=1}^{m_{2}}$ from $P_{2}$ with
$P_{2,ij}=P_{1,ij}-\varepsilon_{2}$, where

$\varepsilon_{2}=\begin{cases}1/(5\log m)&\text{if }i,j\in S\\\ 0&\text{if
}i,j\not\in S\end{cases}\enskip.$

The results are presented in the first row of Figure 4.2. As before, we set
$\varepsilon_{2}=2/(5\log m)$ for $i,j\in S$ and scale the link probability
matrix $P_{1}$ by a factor $\rho=1/4$ to yield sparser networks. The results
are shown in the second row of Figure 4.2.

Figure 4.2. Simulation results for testing networks with a graphon structure.
The first and the second rows are for dense and sparse networks, respectively.
(a) Null rejection rate. (b) Power under the alternative.

We see from Figure 4.2 that SPE-MNBS exhibits superior performance over the
other tests in terms of both null rejection rate and power. We also observe
that SPE-SBM shows a lower rejection rate than the nominal level in the dense
case, which suggests that SPE-SBM is more sensitive to network topologies that
deviate from an SBM. The behaviors of SPE-AVG and DFRO are similar to those in
the first example, as we continue to see subpar performance, especially for
small $m$.

###  4.3. Correlated Erdös-Rényi model

In the third experiment, we study the robustness of the four tests to
dependency. For this, we consider the correlated Erdős-Rényi (ER) model from
[25]. We begin by sampling two independent Erdős-Rényi networks, $A_{1}\sim
ER(n,p_{1})$ and $A_{2}\sim ER(n,p_{2})$. We generate
$\\{A_{1}^{(k)}\\}_{k=1}^{m_{1}}$ with a parameter $\varepsilon_{3}$ as
follows:

$A_{1,ij}^{(k)}\sim\begin{cases}\text{Bernoulli}(\varepsilon_{3})&\text{ if
}A_{1,ij}=1\\\ 0&\text{ if }A_{1,ij}=0\end{cases}\enskip.$

This yields $m_{1}$ networks that are marginally $ER(n,p_{1}\varepsilon_{3})$,
but whose edge sets are correlated. We similarly generate
$\\{A_{2}^{(k)}\\}_{k=1}^{m_{2}}$ conditional on $A_{2}$ with parameter
$\varepsilon_{4}$. We set $\varepsilon_{3}=\varepsilon_{4}=0.8$ and
$p_{1}=0.9$. Under the null hypothesis, we set $p_{2}=p_{1}=0.9$, whereas
$p_{2}=0.83$ under the alternative hypothesis. The results are shown in Figure
4.3.

Figure 4.3. Simulation results for testing networks with a correlated ER
structure. (a) Null rejection rate. (b) Power under the alternative.

We see that DFRO has consistently high power in the alternative setting for
the entire range of $n$, which is only matched for our tests as $n$ increases,
with SPE-AVG outperforming both SPE-SBM and SPE-MNBS. However, the rejection
rate under the null is below the nominal level for DFRO, whereas both SPE-SBM
and SPE-MNBS are very close to $\alpha=0.05$. SPE-AVG has a higher rejection
rate than expected when the sample size is $m=10$, but this improves when
$m=50$. Overall, it appears that SPE-SBM and SPE-MNBS are robust to the
independence violation when $n$ is large.

###  4.4. Beta weight distribution

In our fourth experiment, we evaluate the tests on weighted networks by
generating edge weights from Beta distributions. Since SPE-MNBS and DFRO are
for binary networks, we only consider SPE-AVG and SPE-SBM tests in this
experiment.

Each network is constructed with two communities, where the membership of the
$i$th node is given by

$M(i)=I(1\leq i\leq n/2)+2I(n/2+1\leq i\leq n)\enskip.$

When nodes $i$ and $j$ are in the same community,
$A_{1,ij}^{(k)}\sim\mathrm{Beta}(x_{1},x_{2})$ and
$A_{2,ij}^{(k)}\sim\mathrm{Beta}(x_{1}+\varepsilon_{5},x_{2}+\varepsilon_{5})$.
Otherwise, $A_{1,ij}^{(k)}\sim\mathrm{Beta}(y_{1},y_{2})$ and
$A_{2,ij}^{(k)}\sim\mathrm{Beta}(y_{1}+\varepsilon_{5},y_{2}+\varepsilon_{5})$,
where $x_{1}=2,x_{2}=8,y_{1}=4,y_{2}=1$, $1\leq i<j\leq n$. We set
$\varepsilon_{5}=0$ and $\varepsilon_{5}=1/(5\log m)$ for the null and
alternative hypotheses, respectively. Again, we scale the weighted adjacency
matrices by a factor $\rho=1/4$ to obtain sparser networks. The results are
shown in Figure 4.4.

Figure 4.4. Simulation results for testing weighted networks with a community
structure and Beta-distributed edge weights. The first and the second rows are
for dense and sparse networks, respectively. (a) Null rejection rate. (b)
Power under the alternative.

Figure 4.4 shows that SPE-SBM has good performance under both hypotheses even
though the networks are weighted. However, SPE-AVG suffers from a high type-
one error especially when the sample size $m$ is small. Both tests have
slightly higher power in the dense setting, but the difference is small.

###  4.5. Multiple sample testing

Finally, we evaluate our multiple-sample test from Section 3.2 using the test
statistic in (3.9). As in our first simulation, we consider the same SBM
structure with block matrix given in (4.1). In the null case, we generate
three groups, $\\{A_{s}^{(k)}\\}_{k=1}^{m_{s}}$, $s=1,2,3$, using
$\varepsilon_{1}=0$. In the second case, we let the first two groups be
sampled from $P_{\text{SBM}}$ with $\varepsilon_{1}=0$, while the third group
uses $\varepsilon_{1}=1/(5\log m)$. In the third case, we sample the first
group of networks from $P_{\text{SBM}}$ with $\varepsilon_{1}=0$, while the
second and third groups are with $\varepsilon_{1}=1/(5\log m)$ and $-1/(5\log
m)$, respectively. The results are given in Figure 4.5.

Figure 4.5. Simulation results for testing multiple groups of binary networks
with SBM community structure. The first row is $m=10$ and the second row is
$m=50$. (a) Null rejection rate when all three groups are the same. (b) Power
when two of the groups are the same and the third group is different. (c)
Power when all three groups are different.

We see that SPE-AVG has the highest power for small $n$ in both alternative
scenarios. However, SPE-AVG has a very inflated null rejection rate for
$m=10$, which is significantly improved, though still slightly inflated, for
$m=50$. On the other hand, SPE-SBM and SPE-MNBS have null rejection rates that
oscillate around the nominal rate of $\alpha=0.05$, indicating that the gamma
approximation for the test statistic, which is a sum of dependent $\chi^{2}$
statistics, is very accurate. Under the alternative, these methods perform
similarly: when two groups are the same and one is different, the power
quickly converges to one, and when all three groups are different, the power
almost begins at one.

##  5\. Real data examples

In this section, we apply our tests on three real datasets representing three
different settings of interest to the biological research community. The first
two are networks constructed from fMRI data that represent two distinct
streams of fMRI usage, the former being a case/control study and the latter
being task-based. The third dataset is derived from microbial measurements, an
area in which network-based representations have recently emerged as a popular
technique for studying the bacteria present within a microbiome [19].

In all three cases, the networks are weighted. Therefore, we present results
from our test for weighted networks in Section 3.1. To understand the
performance of our tests for binary networks from Section 2.1, we also present
the results as a function of thresholding the weights to binarize the networks
(as is often done in practice).

###  5.1. Data

Our first dataset comes from the Center for Cognitive Brain Imaging at
Carnegie Mellon University. As part of the StarPlus experiment, fMRI data were
collected for 6 individuals over 40 trials [22]. During each trial of the
experiment, an individual was shown an initial stimulus (either sentence or
picture) for 4 seconds, then a blank screen for 4 seconds, and finally a
second stimulus for up to 4 seconds. The trial ended when the individual
answered whether the sentence correctly described the picture. After a rest
period of 15 seconds, the trial was repeated with a new term. Images were
recorded every 0.5 seconds and thus there are a total of 54 images for each
trial which lasted for 27 seconds. In half of the 40 trials, the picture was
presented first, and in the remaining trials, the sentence was presented
first. We refer to these two datasets as PS and SP, respectively. The data are
publicly available at
http://www.cs.cmu.edu/afs/cs.cmu.edu/project/theo-81/www/.

Brain networks are constructed using the fMRI measurements from this study,
where nodes correspond to the 24 regions of interest (ROIs) given by [16] and
edges represent connectivity between the ROIs. Specifically, for each time
point $T=1,\ldots,54$, we combine the samples from the 6 individuals and
trials in the PS or SP datasets, and estimate the correlation matrix by a
shrinkage estimator proposed in [27]. As a result, we obtain 54 networks from
the PS dataset and 54 networks from the SP dataset, which we call NetPS and
NetSP, respectively. Finally, the fMRI time-series when the individuals view a
sentence consist of the first 16 networks in NetSP and the 17th-32nd networks
in NetPS. The fMRI time-series when the individuals view a picture consist of
the 17th-32nd networks in NetSP and the first 16 networks of NetPS. These
final two groups are referred to as NetS and NetP, respectively.

The StarPlus networks derive from fMRI data in a task-based setting. Another
common use of fMRI data is case/control studies. One such example is the COBRE
dataset from [1], which consists of 124 brain networks for 54 schizophrenics
and 70 controls. These brain networks are also constructed from fMRI data,
where nodes correspond to 264 ROIs defined in [26] and edges represent the
Pearson correlation coefficient between the time series of fMRI activity in
the corresponding ROIs. We use the networks constructed in [2], who provide
further detail on the pre-processing and registration steps used to construct
the networks.

The final dataset is from [8]. The authors tracked the microbiomes (MBs) of 37
women over the course of their pregnancies, 11 of which resulted in preterm
delivery. This study resulted in a time series of operational taxonomic unit
(OTU) tables, which are counts of taxa present in each sample. MB networks can
be constructed from such data, where nodes correspond to taxa and edges
represent the co-occurrence of taxa in the corresponding OTUs. We use the
networks constructed in [17] and note that network analysis is an emerging
tool in the MB literature.

###  5.2. Results for weighted tests

We begin by applying our tests for weighted networks from Section 3.1 to all
three datasets. We test if the groups defined by their respective labels –
picture/sentence, schizophrenic/control, preterm/term – are different. To do
so, we specify a null hypothesis which says the underlying random
distributions are equal against the alternative that says they are different.
We refer to this as the “alternative setting” because the two samples are
different with respect to their group label.

For $\alpha=0.05$ and $Q=1000$, we find that for the StarPlus networks, SPE-
AVG and SPE-SBM correctly reject the null with reject rates of 1 and 0.71,
separately, which are consistent with previous research on distinguishing the
cognitive states of looking at a picture and a sentence [22, 29, 21]. For both
the COBRE and MB datasets, we find a rejection rate of one for both SPE-AVG
and SPE-SBM.

Next, we perform an in silico experiment with the real data by subsampling
within one of the classes. We refer to this as the “null setting.” The
rationale for this setup is that we do not actually know if the different
groups are generated by different underlying distributions, e.g one for
schizophrenic and another for non-schizophrenic. Therefore, we want to check
if the null rejection rate is close to the nominal level in an experiment
where all of the networks are from the same group.

To do so, we test the entire NetP, non-schizophrenic, and term delivery groups
against a subsample (with half of the original sample size) of the same group
for the StarPlus, COBRE, and MB datasets, respectively.

After $1000$ random subsamples of the networks and $Q=1$ for each subsample,
for the StarPlus networks, SPE-SBM fails to reject the null hypothesis with
reject rates of 0.003, which is expected since the samples are drawn from the
same population. However, SPE-AVG rejects the null with an inflated rate of
0.096. For the COBRE networks, we obtain null rejection rates of 0.78 and 0.68
for SPE-AVG and SPE-SBM, respectively. The null rejection rates are improved
for the MB networks with 0.63 and 0.46 for SPE-AVG and SPE-SBM, respectively,
but still very inflated compared to the nominal $\alpha=0.05$.

We speculate that this is happening because even within one class, there is a
lot of variation. That is, a subsample of brain networks with schizophrenia
may look very different from another sample of brain networks with
schizophrenia because we are not controlling for potential factors such as age
and sex. We refer to this issue as having too much heterogeneity within a
class. This heterogeneity could lead to inflated null rejection rates because
the underlying distributions of the two samples are different, but the
difference is not the one we are trying to isolate.

###  5.3. Results for binary tests

As the results for the weighted tests showed inflated rejection rates in our
simulated null setting, there is reason to believe the networks are too
heterogeneous within each class. Furthermore, many of the weights could
represent spurious correlations. Therefore, this is a setting in which
binarizing the weights could improve the signal-to-noise ratio. This idea is
related to a common problem in the neuroscience literature related to the
issue of sensitivity to thresholding edges [14, 11].

To evaluate this, we apply the binary tests from Section 2.1 by binarizing the
weights, which are all correlation values in $[-1,1]$, based on thresholding
their magnitude. This threshold relates directly to the density of the
networks. With the same procedures in Section 5.2, the results are given in
Figures 5.1, 5.2, and 5.3. The dashed lines for the null rejection rate in
these figures all indicate the nominal level of 0.05.

Figure 5.1. Figures (a) and (b) show the null rejection rate and power,
respectively, for different thresholds for binarizing the StarPlus networks.
Figure 5.2. Figures (a) and (b) show the null rejection rate and power,
respectively, for different thresholds for binarizing the COBRE networks.
Figure 5.3. Figures (a) and (b) show the null rejection rate and power,
respectively, for different thresholds for binarizing the MB networks.

The plots illustrate the tradeoff between the false positive rate in our null
setting and the true positive rate in our alternative setting, which both are
functions of the threshold. As the threshold for an edge increases, the
network becomes more sparse, resulting in a higher rejection rate in our null
setting. For thresholds above 0.6, some of the networks became too sparse,
even resulting in some null graphs. On the other hand, for a low threshold,
there is less power to detect a difference in our alternative setting. Such
curves as a function of threshold could provide practitioners a way to
understand the signal-to-noise ratio of their edge weights.

For the COBRE and MB networks, we have high power for a wide range of
threshold values, which is consistent with our findings using the weighted
networks directly. However, we also see a low rejection rate in our null
setting, especially for a threshold of 0.1, which seems to provide the best
tradeoff. This suggests that the signal-to-noise ratio in the weights is too
low, which can be mitigated through thresholding. For the StarPlus networks, a
threshold between 0.1 and 0.2 seems to provide the best balance between signal
and noise for SPE-AVG, whereas 0.55 is better for SPE-SBM and SPE-MNBS.

##  6\. Discussion

In this work, we proposed a new spectral-based statistic for hypothesis
testing of populations of networks, which applies to both binary and weighted
networks under a very general framework. The test statistics are simple,
computationally friendly, and theoretically supported by our derivation of
both the limiting null distribution and asymptotic power guarantees. We have
demonstrated our method through extensive simulation study as well as real
data analysis. Future work will focus on exploring spectral-based methods for
studying inference problems for networks with additional constraints or
structures such as directed networks.

## References

  * [1] CJ Aine, H Jeremy Bockholt, Juan R Bustillo, José M Cañive, Arvind Caprihan, Charles Gasparovic, Faith M Hanlon, Jon M Houck, Rex E Jung and John Lauriello  “Multimodal neuroimaging in schizophrenia: Description and dissemination”  In _Neuroinformatics_ 15.4 Springer, 2017, pp. 343–364
  * [2] Jesús D Arroyo Relión, Daniel Kessler, Elizaveta Levina and Stephan F Taylor  “Network classification with applications to brain connectomics”  In _the Annals of Applied Statistics_ 13.3 Institute of Mathematical Statistics, 2019, pp. 1648–1677
  * [3] Zhidong Bai and Jack W. Silverstein  “Spectral Analysis of Large Dimensional Random Matrices”  New York; London: Springer, 2010
  * [4] Peter J. Bickel and Purnamrita Sarkar  “Hypothesis testing for automated community detection in networks”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 78.1, 2016, pp. 253–273
  * [5] Stanley H. Chan and Edoardo M. Airoldi  “A Consistent Histogram Estimator for Exchangeable Graph Models”  In _Proceedings of the 31st International Conference on Machine Learning_ , 2014
  * [6] Li Chen, Lizhen Lin and Jie Zhou  “A Hypothesis Testing for Large Weighted Networks With Applications to Functional Neuroimaging Data”  In _IEEE Access_ 8 IEEE, 2020, pp. 191815–191825
  * [7] Li Chen, Jie Zhou and Lizhen Lin  “Hypothesis testing for populations of networks”  In _arXiv:1911.03783_ , 2019  arXiv: https://arxiv.org/abs/1911.03783.pdf.
  * [8] Daniel B DiGiulio, Benjamin J Callahan, Paul J McMurdie, Elizabeth K Costello, Deirdre J Lyell, Anna Robaczewska, Christine L Sun, Daniela SA Goltsman, Ronald J Wong and Gary Shaw  “Temporal and spatial variation of the human microbiota during pregnancy”  In _Proceedings of the National Academy of Sciences_ 112.35 National Acad Sciences, 2015, pp. 11060–11065
  * [9] Zhishan Dong, Shuangshuang Wang and Qun Liu  “Spectral based hypothesis testing for community detection in complex networks”  In _Information Sciences_ 512 Elsevier, 2020, pp. 1360–1371
  * [10] Alberto Ferrari  “A note on sum and difference of correlated chi-squared variables”  In _arXiv:1906.09982_ , 2019  URL: https://arxiv.org/abs/1906.09982.pdf.
  * [11] Kathleen A Garrison, Dustin Scheinost, Emily S Finn, Xilin Shen and R Todd Constable  “The (in)stability of functional brain network measures across thresholds”  In _NeuroImage_ 118 Elsevier, 2015, pp. 651–661
  * [12] D. Ghoshdastidar and U. Luxburg  “Practical Methods for Graph Two-Sample Testing”  In _Proceedings Neural Information Processing Systems_ , 2018
  * [13] Debarghya Ghoshdastidar, Maurilio Gutzeit, Alexandra Carpentier and Ulrike Von Luxburg  “Two-sample hypothesis testing for inhomogeneous random graphs”  In _the Annals of Statistics_ 48.4, 2020, pp. 2208–2229
  * [14] Cedric E Ginestet, Arnaud P Fournel and Andrew Simmons  “Statistical network analysis for functional MRI: summary networks and group comparisons” Article 51  In _Frontiers in Computational Neuroscience_ 8 Frontiers, 2014
  * [15] Cedric E Ginestet, Jun Li, Prakash Balachandran, Steven Rosenberg and Eric D Kolaczyk  “Hypothesis testing for network data in functional neuroimaging”  In _The Annals of Applied Statistics_ 11.2 Institute of Mathematical Statistics, 2017, pp. 725–750
  * [16] Rebecca A Hutchinson, Radu Stefan Niculescu, Timothy A Keller, Indrayana Rustandi and Tom M Mitchell  “Modeling fMRI data generated by overlapping cognitive processes with unknown onsets using Hidden Process Models”  In _NeuroImage_ 46.1 Elsevier, 2009, pp. 87–104
  * [17] Nathaniel Josephs, Lizhen Lin, Steven Rosenberg and Eric D Kolaczyk  “Bayesian classification, anomaly detection, and survival analysis using network inputs with application to the microbiome”  In _arXiv:2004.04765_ , 2020  URL: https://arxiv.org/abs/2004.04765.pdf.
  * [18] Eric D. Kolaczyk, Lizhen Lin, Steven Rosenberg, Jackson Walters and Jie Xu  “Averages of unlabeled networks: Geometric characterization and asymptotic behavior”  In _the Annals of Statistics_ 48.1 The Institute of Mathematical Statistics, 2020, pp. 514–538
  * [19] Mehdi Layeghifard, David M Hwang and David S Guttman  “Disentangling interactions in the microbiome: A network perspective”  In _Trends in Microbiology_ 25.3 Elsevier, 2017, pp. 217–228
  * [20] László Lovász  “Large Networks and Graph Limits”  Providence, RI, USA: American Mathematical Society, 2012
  * [21] Tom M Mitchell, Rebecca Hutchinson, Marcel Adam Just, Radu S Niculescu, Francisco Pereira and Xuerui Wang  “Classifying instantaneous cognitive states from fMRI data”  In _AMIA Annual Symposium Proceedings_ 2003, 2003, pp. 465–469  American Medical Informatics Association
  * [22] Tom M Mitchell, Rebecca Hutchinson, Radu S Niculescu, Francisco Pereira, Xuerui Wang, Marcel Just and Sharlene Newman  “Learning to decode cognitive states from brain images”  In _Machine Learning_ 57.1-2 Springer, 2004, pp. 145–175
  * [23] Soumendu Sundar Mukherjee, Purnamrita Sarkar and Lizhen Lin  “On clustering network-valued data”  In _Advances in Neural Information Processing Systems_ , 2017, pp. 7071–7081
  * [24] Andrew Y. Ng, Michael I. Jordan and Yair Weiss  “On Spectral Clustering: Analysis and an algorithm”  In _Advances in Neural Information Processing Systems 14_ Cambridge, MA, USA: MIT Press, 2002, pp. 849–856
  * [25] Pedram Pedarsani and Matthias Grossglauser  “On the privacy of anonymized networks”  In _Proceedings of the 17th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2011, pp. 1235–1243
  * [26] Jonathan D Power, Alexander L Cohen, Steven M Nelson, Gagan S Wig, Kelly Anne Barnes, Jessica A Church, Alecia C Vogel, Timothy O Laumann, Fran M Miezin and Bradley L Schlaggar  “Functional network organization of the human brain”  In _Neuron_ 72.4 Elsevier, 2011, pp. 665–678
  * [27] Juliane Schäfer and Korbinian Strimmer  “A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics” Article 32  In _Statistical Applications in Genetics and Molecular Biology_ 4.1 De Gruyter, 2005
  * [28] Minh Tang, Avanti Athreya, Daniel L. Sussman, Vince Lyzinski and Carey E. Priebe  “A semiparametric two-sample hypothesis testing problem for random graphs”  In _Journal of Computational and Graphical Statistics_ 26.2 Taylor & Francis, 2017, pp. 344–354
  * [29] Xuerui Wang, Rebecca Hutchinson and Tom M Mitchell  “Training fMRI classifiers to detect cognitive states across multiple human subjects”  In _Advances in Neural Information Processing Systems 16_ MIT Press, 2004, pp. 709–716
  * [30] Yuan Zhang, Elizaveta Levina and Ji Zhu  “Estimating network edge probabilities by neighbourhood smoothing”  In _Biometrika_ 104.4 Oxford University Press, 2017, pp. 771–783
  * [31] Zifeng Zhao, Li Chen and Lizhen Lin  “Change-point detection in dynamic networks via graphon estimation”  In _arXiv:1908.01823_ , 2019  URL: https://arxiv.gg363.site/pdf/1908.01823.pdf.

##  Appendix A Background on spectral results

Here, we review some results on spectral properties of inhomogeneous networks.
We begin with the generalized Wigner matrix.

######  Definition 1 (Generalized Wigner matrix [3]).

A generalized Wigner matrix is a Hermitian random matrix $W$ whose entries on
or above the diagonal are independent.

Let $W$ be an $n\times n$ generalized Wigner matrix with eigenvalues
$\lambda_{i}$ for $i=1,\cdots,n$. The empirical spectral distribution (ESD) of
$W$ is defined as

$F_{n}(x)=\frac{1}{n}\\#\\{\lambda_{i}\leq x,1\leq i\leq n\\}\enskip,$

where $\\#$ denotes counting measure. It was proved in [3] that under some
assumptions, the ESD of a normalized $W$ approximates the semicircle law $F$:

$F(x)=\frac{1}{\sqrt{2\pi}}\sqrt{4-x^{2}},\ -2\leq x\leq 2\enskip.$

Specifically, the result is stated in the following lemma.

######  Lemma 1 (Theorem 2.9 in [3]).

Suppose that $W=\frac{1}{\sqrt{n}}X$ is an $n\times n$ Wigner matrix and the
entries above or on the diagonal of $X$ are independent but may be dependent
on $n$ and may not be necessarily identically distributed. Assume that all the
entries of $X$ are of mean zero and variance 1, and satisfy the condition
that, for any constant $\eta>0$,

(A.1)
$\lim_{n\to\infty}\frac{1}{n^{2}}\sum_{ij}\mathrm{E}(|X_{ij}|^{2})I(|X_{ij}|\geq\eta\sqrt{n})=0\enskip.$

Then the ESD of $W$ converges to the semicircular law almost surely.

Next, we give a useful theorem focusing on $F_{n}$, the ESD of a normalized
$W$ in the following.

Let $\mathcal{U}$ be an open set of the real line that contains the interval
$[-2,2]$, which is the support of the semicircular law $F$. Define
$\mathcal{F}$, the set of functions $f:\mathcal{U}\to\mathbb{R}$, where
$\mathbb{R}$ is the set of real numbers. Consider the empirical process
$G_{n}=\\{G_{n}(f)\\}$ given by:

(A.2) $G_{n}(f)=n\int_{-\infty}^{\infty}f(x)[F_{n}-F]dx,\quad
f\in\mathcal{F}\enskip.$

It is shown in [3] that $G_{n}$ converges to a Gaussian process under the
following moment conditions:

  1. (1)

For all $i$, $\mathrm{E}(|X_{ii}|^{2})=\sigma^{2}>0$, for all $i<j$,
$\mathrm{E}(|X_{ij}|^{2})=1$.

  2. (2)

(Homogeneity of fourth moments) $M=\mathrm{E}(|X_{ij}|^{4})$ for $i\neq j$.

  3. (3)

(Uniform tails) For any $\eta>0$, as $n\to\infty$,

$\lim_{n\to\infty}\frac{1}{\eta^{4}n^{2}}\sum_{ij}\mathrm{E}\big{(}|X_{ij}|^{4}I(|X_{ij}|\geq\eta\sqrt{n})\big{)}=0\enskip.$

For any $f\in\mathcal{F}$ and any integer $l\geq 0$, define

$\tau_{l}(f)=\frac{1}{2\pi}\int_{-\pi}^{\pi}f(2\cos(\theta))e^{il\theta}d\theta\enskip.$

Setting $\beta=\mathrm{E}((X_{12}-1)^{2})-2$, there is the following theorem
on the convergence of $G_{n}$.

######  Lemma 2 (Theorem 9.2 in [3]).

Under conditions $(1)$–$(3)$, the spectral empirical process
$G_{n}=(G_{n}(f))$ indexed by the set of analytic functions $\mathcal{F}$
converges weakly in finite dimension to a Gaussian process
$G=\\{G(f),f\in\mathcal{F}\\}$ with mean function $\mathrm{E}{(G(f))}$ given
by

$\frac{1}{4}(f(2)+f(-2))-\frac{1}{2}\tau_{0}(f)+(\sigma^{2}-2)\tau_{2}(f)+\beta\tau_{4}(f)\enskip,$

and the covariance function
$c(f,g)=\mathrm{E}\big{(}(G(f)-\mathrm{E}(G(f))(G(g)-\mathrm{E}(G(g))\big{)}$
given by

$\frac{1}{4\pi^{2}}\int_{-2}^{2}\int_{-2}^{2}f^{\prime}(t)g^{\prime}(t)V(t,s)dtds\enskip,$

where

$V(t,s)=(\sigma^{2}-2+\frac{1}{2}\beta
ts)\sqrt{(4-t^{2})(4-s^{2})}+2\log\left(\frac{4-ts+\sqrt{(4-t^{2})(4-s^{2})}}{4-ts-\sqrt{(4-t^{2})(4-s^{2})}}\right)\enskip.$

##  Appendix B Proofs of theorems from Section 2

###  B.1. Proof of Theorem 1

One can see that under the null hypothesis of $P_{1}=P_{2}$, $Z$ is a Wigner
matrix satisfying $\mathrm{E}(Z_{ij})=0$ and $\text{Var}(Z_{ij})=1/n$. The
remainder of the proof proceeds in three steps.

First, we show that the ESD of $Z$ converges to the semicircular law $F$
almost surely. Let $X=\sqrt{n}Z$. From Lemma 1, it is sufficient to prove that
condition (A.1) is satisfied. For any constant $\eta>0$, we have

$\displaystyle\quad\frac{1}{n^{2}}\sum_{ij}\mathrm{E}(|X_{ij}|^{2})I(|X_{ij}|\geq\eta\sqrt{n})$
$\displaystyle=\frac{1}{n^{2}}\sum_{ij}I(|X_{ij}|\geq\eta\sqrt{n})$
$\displaystyle\leq\max_{ij}\\{I(|X_{ij}|\geq\eta\sqrt{n})\\}$
$\displaystyle=\max_{i\neq j}\\{I\left(|\bar{A}_{1,ij}-\bar{A}_{2,ij}|\geq\eta
T\right),I(1\geq\eta\sqrt{n})\\}$ $\displaystyle\quad\to 0,\quad\text{as}\
n\to\infty\enskip,$

where
$T=\sqrt{n\left(\frac{1}{m_{1}}P_{1,ij}(1-P_{1,ij})+\frac{1}{m_{2}}P_{2,ij}(1-P_{2,ij})\right)}$,
the last equality is due to $X_{ii}=1$, and the last convergence comes from
the boundedness of $|\bar{A}_{1,ij}-\bar{A}_{2,ij}|$, whereas
$n/m_{u}\to\infty$ for $u=1,2$ and $\eta\sqrt{n}\to\infty$.

Next, we verify that $X$ satisfies conditions $(1)$–$(3)$ of Lemma 2. For
simplicity, we still denote the ESD of $Z$ as $F_{n}$. Then, the empirical
process $G_{n}$ defined in (A.2) converges to a Gaussian process from Lemma 2
and the asymptotic normality of $\theta$ given by (2.5) follows directly.

Condition $(1)$ holds since
$\mathrm{E}(|X_{ij}|^{2})=n\mathrm{E}(|Z_{ij}|^{2})=1$ for all $i$ and $j$.

By the central limit theorem, if $i\neq j$, then $Z_{ij}$ converges to the
same normal distribution $\mathcal{N}(0,1/n)$ when $m_{1}$ and $m_{2}$ are
large. Thus condition $(2)$ is satisfied under the large sample size
assumption.

Finally, we verify condition $(3)$ by applying a Hölder inequality. For any
$\eta>0$,

$\displaystyle\quad\frac{1}{\eta^{4}n^{2}}\sum_{ij}\mathrm{E}\big{(}|X_{ij}|^{4}I(|X_{ij}|\geq\eta\sqrt{n})\big{)}$
$\displaystyle\leq\frac{1}{\eta^{4}n^{2}}\sum_{ij}\big{(}\mathrm{E}(|X_{ij}|^{8})\big{)}^{1/2}\big{(}\mathrm{E}\big{(}I(|X_{ij}|\geq\eta\sqrt{n})\big{)}\big{)}^{1/2}$
$\displaystyle\leq\frac{C_{1}}{\eta^{4}n^{2}}\sum_{ij}\big{(}P\big{(}|X_{ij}|\geq\eta\sqrt{n}\big{)}\big{)}^{1/2}$
$\displaystyle\leq\frac{C_{1}}{\eta^{4}}\max_{i\neq
j}\left\\{\big{(}P\big{(}|X_{ij}|\geq\eta\sqrt{n}\big{)}\big{)}^{1/2},\big{(}P\big{(}1\geq\eta\sqrt{n}\big{)}\big{)}^{1/2}\right\\}$
$\displaystyle=\frac{C_{1}}{\eta^{4}}\max_{i\neq
j}\left\\{\left(P\left(|\bar{A}_{1,ij}-\bar{A}_{2,ij}|\geq\eta
T\right)\right)^{1/2},\big{(}P\big{(}1\geq\eta\sqrt{n}\big{)}\big{)}^{1/2}\right\\}$
$\displaystyle\quad\to 0,\quad\text{as}\ n\to\infty\enskip,$

where $C_{1}$ is the upper bound of
$\big{(}\mathrm{E}(|X_{ij}|^{8})\big{)}^{1/2}$ and is a positive constant. The
last convergence also follows from the boundedness of
$|\bar{A}_{1,ij}-\bar{A}_{2,ij}|$ and the infinity trend of $T$ and
$\eta\sqrt{n}$ as $n$ tends to infinity.

Taking $f(x)=x^{3}$ and combining Lemma 2 with (A.2), we have

$\displaystyle G_{n}(f)$
$\displaystyle=n\int_{-\infty}^{\infty}f(x)[F_{n}-F]dx$
$\displaystyle=n\int_{-\infty}^{\infty}f(x)F_{n}dx-n\int_{-\infty}^{\infty}f(x)Fdx$
(B.1) $\displaystyle=\sum_{i=1}^{n}f(\lambda_{i}(Z))\enskip,$

where $\lambda_{i}(Z)$ is the $i$th eigenvalue of $Z$. The third equality
holds because the first term is the expectation of $f(x)$ with respect to the
counting measure $F_{n}$ and the second term is the integral of an odd
function $f(x)F$.

By diagonalizing $Z=U^{T}\Lambda U$, where $U$ is some orthogonal matrix and
$\Lambda$ is a diagonal matrix whose diagonal entries are the eigenvalues of
$Z$. From (B.1), we obtain

$\displaystyle G_{n}(f)$
$\displaystyle=\sum_{i=1}^{n}f(\Lambda_{ii})=\text{Tr}(f(\Lambda))=\text{Tr}(U^{T}f(\Lambda)U)$
$\displaystyle=\text{Tr}(f(U^{T}\Lambda
U))=\text{Tr}(f(Z))=\text{Tr}(Z^{3})\enskip.$

$G(x^{3})$ is the same as in [9], and hence its mean and variance are known to
be $\mathrm{E}(G(x^{3}))=0$ and $\text{Var}(G(x^{3}))=15$, respectively, which
completes the proof.

###  B.2. Proof of Theorem 2

First, for any $i\neq j$, we have

(B.2)
$\hat{Z}_{ij}=\frac{\sqrt{\frac{1}{m_{1}}P_{1,ij}(1-P_{1,ij})+\frac{1}{m_{2}}P_{2,ij}(1-P_{2,ij})}}{\sqrt{\frac{1}{m_{1}}\hat{P}_{1,ij}(1-\hat{P}_{1,ij})+\frac{1}{m_{2}}\hat{P}_{2,ij}(1-\hat{P}_{2,ij})}}Z_{ij}\enskip.$

Recall that under the null hypothesis, $P_{1}=P_{2}$. By Taylor expanding the
numerator in (B.2), we obtain

$\displaystyle\quad\sqrt{\frac{1}{m_{1}}P_{1,ij}(1-P_{1,ij})+\frac{1}{m_{2}}P_{2,ij}(1-P_{2,ij})}$
$\displaystyle=\sqrt{\frac{m_{1}+m_{2}}{m_{1}m_{2}}}\sqrt{P_{1,ij}(1-P_{1,ij})}$
$\displaystyle=\sqrt{\frac{m_{1}+m_{2}}{m_{1}m_{2}}}\Big{(}\sqrt{\hat{P}_{1,ij}\big{(}1-\hat{P}_{1,ij}\big{)}}+O\big{(}P_{1,ij}-\hat{P}_{1,ij}\big{)}\Big{)}$
$\displaystyle\leq\sqrt{\frac{m_{1}+m_{2}}{m_{1}m_{2}}}\Big{(}\sqrt{\hat{P}_{1,ij}\big{(}1-\hat{P}_{1,ij}\big{)}}+o_{p}(1)\Big{)}$
$\displaystyle=\sqrt{\frac{1}{m_{1}}\hat{P}_{1,ij}\big{(}1-\hat{P}_{1,ij}\big{)}+\frac{1}{m_{2}}\hat{P}_{1,ij}\big{(}1-\hat{P}_{1,ij}\big{)}}+\sqrt{\frac{m_{1}+m_{2}}{m_{1}m_{2}}}o_{p}(1)\enskip,$

where the third equality is obtained under the assumption that
$\max_{ij}|\hat{P}_{u,ij}-P_{u,ij}|=o_{p}(1)$.

Without loss of generality, assume
$\hat{P}_{1,ij}\big{(}1-\hat{P}_{1,ij}\big{)}\leq\hat{P}_{2,ij}\big{(}1-\hat{P}_{2,ij}\big{)}$.
Then

(B.3)
$\displaystyle\quad\sqrt{\frac{1}{m_{1}}P_{1,ij}(1-P_{1,ij})+\frac{1}{m_{2}}P_{2,ij}(1-P_{2,ij})}$
$\displaystyle\leq\sqrt{\frac{1}{m_{1}}\hat{P}_{1,ij}\big{(}1-\hat{P}_{1,ij}\big{)}+\frac{1}{m_{2}}\hat{P}_{2,ij}\big{(}1-\hat{P}_{2,ij}\big{)}}+\sqrt{\frac{m_{1}+m_{2}}{m_{1}m_{2}}}o_{p}(1)\enskip.$

Similarly, we have

(B.4)
$\begin{split}&\quad\sqrt{\frac{1}{m_{1}}\hat{P}_{1,ij}\big{(}1-\hat{P}_{1,ij}\big{)}+\frac{1}{m_{2}}\hat{P}_{2,ij}\big{(}1-\hat{P}_{2,ij}\big{)}}\\\
&\leq\sqrt{\frac{1}{m_{1}}P_{1,ij}(1-P_{1,ij})+\frac{1}{m_{2}}P_{2,ij}(1-P_{2,ij})}+\sqrt{\frac{m_{1}+m_{2}}{m_{1}m_{2}}}o_{p}(1)\enskip.\end{split}$

From (B.3) and (B.4), we have

(B.5)
$\begin{split}&\quad\sqrt{\frac{1}{m_{1}}P_{1,ij}(1-P_{1,ij})+\frac{1}{m_{2}}P_{2,ij}(1-P_{2,ij})}\\\
&=\sqrt{\frac{1}{m_{1}}\hat{P}_{1,ij}\big{(}1-\hat{P}_{1,ij}\big{)}+\frac{1}{m_{2}}\hat{P}_{2,ij}\big{(}1-\hat{P}_{2,ij}\big{)}}+\sqrt{\frac{m_{1}+m_{2}}{m_{1}m_{2}}}o_{p}(1)\enskip.\end{split}$

Then, combining (B.5) with (B.2), and noticing that when $i=j$,
$\hat{Z}_{ij}=Z_{ij}=B_{ij}$, we can write

(B.6) $\hat{Z}-B=(J+H)\circ(Z-B)\enskip,$

where $\circ$ denotes the Hadamard product, $J$ is the $n\times n$ all-ones
matrix, and $H$ is an $n\times n$ matrix with entries
$H_{ij}=\max_{u=1,2}O(\hat{P}_{u,ij}-P_{ij})=o_{p}(1)$.

Moreover, $H\circ B$ is a zero matrix, thus an immediate consequence of (B.6)
is

$\hat{Z}=Z\circ(J+H)\enskip.$

We further have

$\text{Tr}(\hat{Z})^{3}-\text{Tr}(Z^{3})=3\text{Tr}\big{(}Z^{2}(Z\circ
H)\big{)}+3\text{Tr}\big{(}Z(Z\circ H)^{2}\big{)}+\text{Tr}\big{(}(Z\circ
H)^{3}\big{)}\enskip.$

Next, we prove that $\text{Tr}\big{(}Z^{2}(Z\circ
H)\big{)}=\sum_{i,k,l}Z_{ik}Z_{kl}Z_{li}H_{li}=o_{p}(1)$. The summation can be
divided into two parts in terms of $(i,k,l)$. In particular, let
$\\{1,\cdots,n\\}^{3}=S_{1}\cup S_{2},\ S_{1}\cap S_{2}=\emptyset$. In set
$S_{1}$, $Z_{ik}Z_{kl}Z_{li}\geq 0$, and in set $S_{2}$,
$Z_{ik}Z_{kl}Z_{li}<0$.

$\displaystyle\text{Tr}\big{(}Z^{2}(Z\circ H)\big{)}$
$\displaystyle=\sum_{i,k,l\in S_{1}}Z_{ik}Z_{kl}Z_{li}H_{li}+\sum_{i,k,l\in
S_{2}}Z_{ik}Z_{kl}Z_{li}H_{li}$ $\displaystyle\leq\sum_{i,k,l\in
S_{1}}Z_{ik}Z_{kl}Z_{li}\max_{li}\\{H_{li}\\}+\sum_{i,k,l\in
S_{2}}Z_{ik}Z_{kl}Z_{li}\min_{li}\\{H_{li}\\}$ (B.7)
$\displaystyle=o_{p}(1)\enskip.$

The last equality comes from as following. From Theorem 1,
$\text{Tr}(Z^{3})=O_{p}(1)$, when $S_{1}=\emptyset$ or $S_{2}=\emptyset$,
which means $Z_{ij}<0$ or $Z_{ij}\geq 0$ for any $i$ and $j$, inequality (B.7)
holds straightforwardly. When $S_{1}\neq\emptyset$ and $S_{2}\neq\emptyset$,
$\sum_{i,k,l\in S_{1}}Z_{ik}Z_{kl}Z_{li}$ can be viewed as a part of the total
sum $\sum_{i,k,l}Z_{ik}Z_{kl}Z_{li}$ when all $Z_{ij}\geq 0$, so
$\sum_{i,k,l\in S_{1}}Z_{ik}Z_{kl}Z_{li}=O_{p}(1)$. Similarly, $\sum_{i,k,l\in
S_{2}}Z_{ik}Z_{kl}Z_{li}=O_{p}(1)$. Thus, inequality (B.7) still holds.

By the similar method, we have $\text{Tr}\big{(}Z(Z\circ
H)^{2}\big{)}=o_{p}(1)$ and $\text{Tr}\big{(}(Z\circ H)^{3}\big{)}=o_{p}(1)$.
As a result,

$\text{Tr}(\hat{Z})^{3}=\text{Tr}(Z^{3})+o_{p}(1)\to\text{Tr}(Z^{3})\quad\text{as}\quad
n\to\infty\enskip,$

and our proof is complete.

###  B.3. Proof of Theorem 3

As in the proof of Theorem 2, $\hat{Z}=Z\circ(J+H)$ and

(B.8) $\text{Tr}(\hat{Z}^{3})=\text{Tr}(Z^{3})+o_{p}(1)\enskip.$

Let $\tilde{Z}$ be an $n\times n$ matrix with entries

$\tilde{Z}_{ij}=\begin{cases}\frac{(\bar{A}_{1,ij}-\bar{A}_{2,ij})-(P_{1,ij}-P_{2,ij})}{\sqrt{n\left(\frac{1}{m_{1}}P_{1,ij}(1-P_{1,ij})+\frac{1}{m_{2}}P_{2,ij}(1-P_{2,ij})\right)}}&\text{
if }i\neq j\\\ B_{ij}&\text{ if }i=j\end{cases}\enskip.$

Noting the definitions of $Z$ and $Z^{\prime\prime}$ in (2.2) and (2.9),
respectively, it is obvious that $\tilde{Z}=Z-Z^{\prime\prime}$ and
$\tilde{Z}$ is a Wigner matrix. Thus,

(B.9)
$\displaystyle\text{Tr}(\tilde{Z}^{3})=\text{Tr}(Z^{3}-(Z^{\prime\prime})^{3}+3Z(Z^{\prime\prime})^{2}-3Z^{2}Z^{\prime\prime})\enskip.$

First, similar with the proof of Theorem 4, we have

(B.10)
$\text{Tr}(\tilde{Z}^{3})/\sqrt{15}\overset{d}{\to}\mathcal{N}(0,1)\enskip.$

Next, we calculate
$\text{Tr}(3Z(Z^{\prime\prime})^{2}-3Z^{2}Z^{\prime\prime})$:

(B.11)
$\displaystyle\text{Tr}(3Z(Z^{\prime\prime})^{2}-3Z^{2}Z^{\prime\prime})=3\sum_{i,k,l}Z_{ik}Z^{\prime\prime}_{kl}Z^{\prime\prime}_{li}-3\sum_{i,k,l}Z_{ik}Z_{kl}Z^{\prime\prime}_{li}\enskip.$

When all $i,\ k,\ l$ are all different,
$\mathrm{E}(Z_{ik}Z^{\prime\prime}_{kl}Z^{\prime\prime}_{li}-Z_{ik}Z_{kl}Z^{\prime\prime}_{li})=0$
due to the independence of the elements and the fact that
$\mathrm{E}(Z_{ij})=Z^{\prime\prime}_{ij}$ for all $i$ and $j$. When at least
two subscripts are the same, we note that
$\mathrm{E}(Z_{ii})=Z^{\prime\prime}_{ii}=0$. Therefore,
$\mathrm{E}(Z_{ik}Z_{kl}Z^{\prime\prime}_{li}-Z_{ik}Z^{\prime\prime}_{kl}Z^{\prime\prime}_{li})=0$
still holds, and hence
$\mathrm{E}(\text{Tr}(3Z(Z^{\prime\prime})^{2}-3Z^{2}Z^{\prime\prime}))=0$.

Next, from (B.11),
$\text{Var}(\text{Tr}(3Z^{2}Z^{\prime\prime}-3Z(Z^{\prime\prime})^{2}))=\max_{u=1,2}O(m_{u}^{2})$.
Using Bernstein’s inequality, we have

(B.12)
$\text{Tr}(3Z(Z^{\prime\prime})^{2}-3Z^{2}Z^{\prime\prime})=\max_{u=1,2}O_{p}(m_{u}^{2})\enskip.$

Finally, we have

(B.13) $\displaystyle|\text{Tr}((Z^{\prime\prime})^{3})|$
$\displaystyle=\left|\sum_{i,k,l}Z^{\prime\prime}_{ik}Z^{\prime\prime}_{kl}Z^{\prime\prime}_{li}\right|\enskip.$

The sum in (B.13) can be divided into two parts in terms of the sets $S_{a}$
and $S_{b}$ as

$\displaystyle|\text{Tr}((Z^{\prime\prime})^{3})|$
$\displaystyle\geq\left|\sum_{(i,k,l)\in
S_{a}}Z^{\prime\prime}_{ik}Z^{\prime\prime}_{kl}Z^{\prime\prime}_{li}\right|-\left|\sum_{(i,k,l)\in
S_{b}}Z^{\prime\prime}_{ik}Z^{\prime\prime}_{kl}Z^{\prime\prime}_{li}\right|$
$\displaystyle\geq an^{3}\min_{(i,k,l)\in
S_{a}}|{Z^{\prime\prime}}_{ik}|^{3}+bn^{3}\min_{(i,k,l)\in
S_{b}}|{Z^{\prime\prime}}_{ik}|^{3}>0,\enskip$

where with a slight abuse of notation, $\min\\{Z_{ik}\\}$ is taken over all
the pairs of indices among $(i,k,l).$ It is obvious to see
$|\text{Tr}((Z^{\prime\prime})^{3})|\geq\max_{u=1,2}O((m_{u}n)^{3/2})$.

Similarly, we can write

$\displaystyle|\text{Tr}((Z^{\prime\prime})^{3})|\geq-bn^{3}\max_{(i,k,l)\in
S_{b}}|{Z^{\prime\prime}}_{ik}|^{3}-an^{3}\max_{(i,k,l)\in
S_{a}}|{Z^{\prime\prime}}_{ik}|^{3}>0\enskip.$

Again have
$|\text{Tr}((Z^{\prime\prime})^{3})|\geq\max_{u=1,2}O((m_{u}n)^{3/2})$.
Therefore, in either case, we see that

(B.14)
$|\text{Tr}((Z^{\prime\prime})^{3})|\geq\max_{u=1,2}O((m_{u}n)^{3/2})\enskip.$

Combining (B.10), (B.12) and (B.14), we are ready to prove the asymptotic
power of our test.

$\displaystyle\quad P(|\theta|>\mu_{\alpha/2})$
$\displaystyle=P\left(\frac{1}{\sqrt{15}}|\text{Tr}(Z^{3})|>\mu_{\alpha/2}\right)$
$\displaystyle=P\left(\frac{1}{\sqrt{15}}|\text{Tr}\big{(}\tilde{Z}^{3}+(Z^{\prime\prime})^{3}-(3Z(Z^{\prime\prime})^{2}-3Z^{2}Z^{\prime\prime})\big{)}|>\mu_{\alpha/2}\right)$
$\displaystyle=P\left(\frac{1}{\sqrt{15}}\text{Tr}(\tilde{Z}^{3})>\mu_{\alpha/2}-\frac{1}{\sqrt{15}}\text{Tr}((Z^{\prime\prime})^{3})+\frac{1}{\sqrt{15}}\text{Tr}\big{(}3Z(Z^{\prime\prime})^{2}-3Z^{2}Z^{\prime\prime}\big{)}\right)$
$\displaystyle\quad+P\left(\frac{1}{\sqrt{15}}\text{Tr}(\tilde{Z}^{3})<-\mu_{\alpha/2}-\frac{1}{\sqrt{15}}\text{Tr}((Z^{\prime\prime})^{3})+\frac{1}{\sqrt{15}}\text{Tr}\big{(}3Z(Z^{\prime\prime})^{2}-3Z^{2}Z^{\prime\prime}\big{)}\right)$
$\displaystyle\quad\to 1\enskip.$

Combing this with (B.8), we have

$\displaystyle P(|\hat{\theta}|>\mu_{\alpha/2})$
$\displaystyle=P(|\theta+o_{p}(1)|>\mu_{\alpha/2})\to 1\enskip,$

which completes our proof.

###  B.4. Proof of Theorem 5

It is not difficult to observe that for any $i\neq j$,

(B.15)
$\hat{Z}_{ij}=\frac{\sqrt{\frac{1}{m_{1}}\Sigma_{1,ij}+\frac{1}{m_{2}}\Sigma_{2,ij}}}{\sqrt{\frac{1}{m_{1}}\hat{\Sigma}_{1,ij}+\frac{1}{m_{2}}\hat{\Sigma}_{2,ij}}}Z_{ij}\enskip.$

Under the null hypothesis $H_{0}$, we have $\Sigma_{1}=\Sigma_{2}$. Without
loss of generality, assume $\hat{\Sigma}_{1,ij}\leq\hat{\Sigma}_{2,ij}$. For
the numerator in (B.15), utilizing a Taylor expansion, we have

$\displaystyle\sqrt{\frac{1}{m_{1}}\Sigma_{1,ij}+\frac{1}{m_{2}}\Sigma_{2,ij}}$
$\displaystyle=\sqrt{\frac{m_{1}+m_{2}}{m_{1}m_{2}}\Sigma_{1,ij}}$
$\displaystyle=\sqrt{\frac{m_{1}+m_{2}}{m_{1}m_{2}}}\Big{(}\sqrt{\hat{\Sigma}_{1,ij}}+O\big{(}\Sigma_{1,ij}-\hat{\Sigma}_{1,ij}\big{)}\Big{)}$
$\displaystyle\leq\sqrt{\frac{\hat{\Sigma}_{1,ij}}{m_{1}}+\frac{\hat{\Sigma}_{2,ij}}{m_{1}}}+\sqrt{\frac{m_{1}+m_{2}}{m_{1}m_{2}}}O\big{(}\Sigma_{1,ij}-\hat{\Sigma}_{1,ij}\big{)}$
(B.16)
$\displaystyle\leq\sqrt{\frac{\hat{\Sigma}_{1,ij}}{m_{1}}+\frac{\hat{\Sigma}_{2,ij}}{m_{1}}}+\sqrt{\frac{m_{1}+m_{2}}{m_{1}m_{2}}}o_{p}(1)\enskip,$

where the last inequality comes from the condition
$\max_{ij}|\hat{\Sigma}_{u,ij}-\Sigma_{u,ij}|=o_{p}(1)$.

Similarly, we have

(B.17)
$\quad\sqrt{\frac{1}{m_{1}}\hat{\Sigma}_{1,ij}+\frac{1}{m_{2}}\hat{\Sigma}_{2,ij}}\leq\sqrt{\frac{1}{m_{1}}\Sigma_{1,ij}+\frac{1}{m_{2}}\Sigma_{2,ij}}+\sqrt{\frac{m_{1}+m_{2}}{m_{1}m_{2}}}o_{p}(1)\enskip.$

From (B.16) and (B.17),

(B.18)
$\quad\sqrt{\frac{1}{m_{1}}\Sigma_{1,ij}+\frac{1}{m_{2}}\Sigma_{2,ij}}=\sqrt{\frac{1}{m_{1}}\hat{\Sigma}_{1,ij}+\frac{1}{m_{2}}\hat{\Sigma}_{2,ij}}+\sqrt{\frac{m_{1}+m_{2}}{m_{1}m_{2}}}o_{p}(1)\enskip.$

Combining (B.18) with (B.15) and noting that $Z_{ii}=\hat{Z}_{ii}=B_{ii}$, we
have

(B.19) $\hat{Z}-B=(J+H)\circ(Z-B)\enskip,$

where $H$ is an $n\times n$ matrix with entries $H_{ij}=o_{p}(1)$.

The remainder of the proof follows just the same as in Theorem 2.

###  B.5. Proof of (3.6)

Denote $G_{1}$ and $G_{2}$ as two $n\times n$ matrices, with elements
$G_{u,ij}=\mathrm{E}\big{(}\big{(}{A_{u,ij}^{(k_{u})}}\big{)}^{2}\big{)}$ for
$u=1,2$, then

$\displaystyle\quad|\hat{\Sigma}_{u,ij}^{2}-\Sigma_{u,ij}^{2}|$
$\displaystyle=\left|\frac{\sum_{k_{u}=1}^{m_{u}}\big{(}{A_{u,ij}^{(k_{u})}}\big{)}^{2}-m_{u}\big{(}\bar{A}_{u,ij}\big{)}^{2}}{m_{u}-1}-G_{u,ij}+P_{u,ij}^{2}\right|$
$\displaystyle=\left|\frac{\sum_{k_{u}=1}^{m_{u}}\big{(}{A_{u,ij}^{(k_{u})}}\big{)}^{2}}{m_{u}}-G_{u,ij}+\frac{\sum_{k_{u}=1}^{m_{u}}\big{(}{A_{u,ij}^{(k_{u})}}\big{)}^{2}}{m_{u}(m_{u}-1)}-\Bigg{(}\big{(}\bar{A}_{u,ij}\big{)}^{2}-P_{u,ij}^{2}+\frac{\big{(}\bar{A}_{u,ij}\big{)}^{2}}{m_{u}-1}\Bigg{)}\right|$
$\displaystyle\leq\left|\frac{\sum_{k_{u}=1}^{m_{u}}\big{(}{A_{u,ij}^{(k_{u})}}\big{)}^{2}}{m_{u}}-G_{u,ij}\right|+\frac{\max\limits_{k_{u}=1,2,\dots,m_{u}}\big{(}{A_{u,ij}^{(k_{u})}}\big{)}^{2}}{m_{u}-1}$
$\displaystyle\quad+|\bar{A}_{u,ij}-P_{u,ij}||\bar{A}_{u,ij}+P_{u,ij}|+\frac{\big{(}\bar{A}_{u,ij}\big{)}^{2}}{m_{u}-1}$
$\displaystyle=O_{p}(m_{u}^{-1/2})+\frac{\max\limits_{k_{u}=1,2,\dots,m_{u}}\big{(}{A_{u,ij}^{(k_{u})}}\big{)}^{2}}{m_{u}-1}+O_{p}(m_{u}^{-1/2})|\bar{A}_{u,ij}+P_{u,ij}|+\frac{\big{(}\bar{A}_{u,ij}\big{)}^{2}}{m_{u}-1}$
(B.20) $\displaystyle=O_{p}(m_{u}^{-1/2}),$

where the second to last equality comes from Bernstein’s inequality and the
last equality is due to the boundedness of distributions $\\{F_{u,ij}\\}$ in
our framework. It is noted that $\hat{\Sigma}_{u,ij}$ and $\Sigma_{u,ij}$ are
bounded, so $|\hat{\Sigma}_{u,ij}+\Sigma_{u,ij}|$ is also bounded. As a
result,
$\max\limits_{i,j}|\hat{\Sigma}_{u,ij}-\Sigma_{u,ij}|=O_{p}(m_{u}^{-1/2})$.