Hypothesis Testing for Network Data with Power Enhancement

TL;DR

This paper introduces a new network comparison method that improves power and flexibility, allowing for global and individual link testing in non-normal network data, with applications in brain connectivity and genomics.

Contribution

It develops a power-enhanced testing procedure for network data that does not require normality and can handle multiple network comparisons simultaneously.

Findings

Power enhancement improves test sensitivity.

Method controls false discovery rate.

Effective in brain connectivity analysis.

Abstract

Comparing two population means of network data is of paramount importance in a wide range of scientific applications. Many existing network inference solutions focus on global testing of entire networks, without comparing individual network links. Besides, the observed data often take the form of vectors or matrices, and the problem is formulated as comparing two covariance or precision matrices under a normal or matrix normal distribution. Moreover, many tests suffer from a limited power under a small sample size. In this article, we tackle the problem of network comparison, both global and simultaneous inferences, when the data come in a different format, i.e., in the form of a collection of symmetric matrices, each of which encodes the network structure of an individual subject. Such data format commonly arises in applications such as brain connectivity analysis and clinical genomics. We no longer require the underlying data to follow a normal distribution, but instead impose some moment conditions that are easily satisfied for numerous types of network data. Furthermore, we propose a power enhancement procedure, and show that it can control the false discovery, while it has the potential to substantially enhance the power of the test. We investigate the efficacy of our testing procedure through both an asymptotic analysis and a simulation study under a finite sample size. We further illustrate our method with an example of brain structural connectivity analysis.

Tables

Table 1. Table 1: The empirical FDR and empirical power, in percentage, for Algorithms 1 and 2 based on 100 data replications. The significance level is α = 5 % 𝛼 percent 5 \alpha=5\% .

	$n_1=n_2=100$			$n_1=n_2=25$
	$0.2q$	$0.15q$	$0.1q$	$0.2q$	$0.15q$	$0.1q$
Bernoulli	Empirical FDR
Algorithm 1	4.0	4.5	4.9	5.8	6.4	7.3
Algorithm 2	2.9	2.8	3.5	3.9	5.0	5.9
	Empirical power
Algorithm 1	88.1	86.7	84.6	44.6	42.2	40.1
Algorithm 2	91.6	91.5	90.7	55.3	54.7	53.2
Bernoulli mixture	Empirical FDR
Algorithm 1	4.1	4.4	4.6	5.9	5.8	7.5
Algorithm 2	1.5	1.6	2.4	3.5	4.1	5.4
	Empirical power
Algorithm 1	89.8	88.6	87.6	42.3	41.4	41.1
Algorithm 2	95.7	95.8	95.6	54.6	54.7	54.5
Transformed Wishart	Empirical FDR
Algorithm 1	4.4	4.9	5.4	4.9	5.5	6.2
Algorithm 2	2.3	2.5	2.7	3.0	3.8	4.2
	Empirical power
Algorithm 1	52.2	55.1	59.3	37.3	39.4	42.0
Algorithm 2	60.1	64.4	69.4	41.6	44.7	49.8

Table 2. Table 2: Significantly different links found by the proposed multiple testing procedure for the KKI-42 dataset. Reported are the significant links identified by the power-enhanced Algorithm 2 but missed by Algorithm 1 .

Differentiating links
r.rostralanteriorcingulate $\leftrightarrow$ r.superiorparietal	l.caudalmiddlefrontal $\leftrightarrow$ l.fusiform
r.fusiform $\leftrightarrow$ r.postcentral	r.lateralorbitofrontal $\leftrightarrow$ l.isthmuscingulate
l.corpuscallosum $\leftrightarrow$ l.lingual	r.paracentral $\leftrightarrow$ l.corpuscallosum
r.precentral $\leftrightarrow$ l.isthmuscingulate	l.lingual $\leftrightarrow$ l.parsopercularis
r.corpuscallosum $\leftrightarrow$ r.inferiortemporal	r.temporalpole $\leftrightarrow$ l.medialorbitofrontal
r.lingual $\leftrightarrow$ l.precuneus	r.parsorbitalis $\leftrightarrow$ r.superiortemporal
r.corpuscallosum $\leftrightarrow$ r.fusiform	r.fusiform $\leftrightarrow$ r.rostralmiddlefrontal
r.lingual $\leftrightarrow$ l.parahippocampal	r.temporalpole $\leftrightarrow$ l.entorhinal
r.frontalpole $\leftrightarrow$ l.rostralmiddlefrontal	r.cuneus $\leftrightarrow$ l.pericalcarine
l.parahippocampal $\leftrightarrow$ l.pericalcarine	l.isthmuscingulate $\leftrightarrow$ l.frontalpole
r.medialorbitofrontal $\leftrightarrow$ l.isthmuscingulate	r.middletemporal $\leftrightarrow$ r.rostralmiddlefrontal
r.parsorbitalis $\leftrightarrow$ r.precuneus	l.superiorfrontal $\leftrightarrow$ l.transversetemporal

Equations

$$H_0:{\mathbf{s}}_1={\mathbf{s}}_2\text{mbox}versus{H}_1:{\mathbf{s}}_1\ne {\mathbf{s}}_2.$$

$$H_{0,i,j}:s_{1,i,j}=s_{2,i,j}\text{mbox}versus{H}_{1,i,j}:s_{1,i,j}\ne s_{2,i,j},\text{mbox}for1\le i<j\le p.$$

$$\text{E}\left[\exp \left\{\eta (S_{d,l,i,j}-s_{d,i,j}{)}^2/\text{Var}(S_{d,l,i,j})\right\}\right]\le K,\text{ for }1\le i<j\le p,l=1,\dots ,n_d.$$

$$\text{E}\left\{{\left|(S_{d,l,i,j}-s_{d,i,j})/\text{Var}(S_{d,l,i,j}{)}^{1/2}\right|}^{4{\gamma }_0+2+ϵ}\right\}\le K,\text{ for }1\le i<j\le p,l=1,\dots ,n_d.$$

$$W_{i,j}={\bar{S}}_{1,i,j}-{\bar{S}}_{2,i,j}.$$

$$V_{1,i,j}=n_1^{-1}\sum _{l=1}^{n_1}(S_{1,l,i,j}-{\bar{S}}_{1,i,j}{)}^2,\text{ and }{V}_{2,i,j}=n_2^{-1}\sum _{l=1}^{n_2}(S_{2,l,i,j}-{\bar{S}}_{2,i,j}{)}^2,$$

$$T_{i,j}=\frac{W_{i,j}}{(V_{1,i,j}/n_1+V_{2,i,j}/n_2{)}^{1/2}},1\le i<j\le p.$$

$$M_n=\underset{1\le i<j\le p}{\max }{T}_{i,j}^2.$$

$${\text{P}}_{H_0}\left(M_n-2\log q+\log \log q\le x\right)\to \exp \left\{-{\pi }^{-1/2}\exp (-x/2)\right\},\text{as}{n}_1,n_2,q\to \infty .$$

$${\Psi }_{\alpha }=I(M_n\ge 2\log q-\log \log q+q_{\alpha }),$$

$$\mathcal{U}(c)=\left\{({\mathbf{s}}_1,{\mathbf{s}}_2):\underset{1\le i<j\le p}{\max }\frac{|s_{1,i,j}-s_{2,i,j}|}{{\left\{\text{Var}(S_{1,l,i,j})/n_1+\text{Var}(S_{2,l,i,j})/n_2\right\}}^{1/2}}\ge c(\log q{)}^{1/2}\right\}.$$

$$\underset{(s_1,s_2)\in \mathcal{U}(2\sqrt{2})}{\inf }\text{P}\left({\Psi }_{\alpha }=1\right)\to 1,\text{as}{n}_1,n_2,q\to \infty .$$

$$\underset{(s_1,s_2)\in U(c_0)}{\inf }\underset{T_{\alpha }\in {\mathcal{T}}_{\alpha }}{\sup }\text{P}\left(T_{\alpha }=1\right)\le 1-\beta ,$$

$$\text{FDP}(h)=\frac{R_0(h)}{R(h)\vee 1},\text{FDR}(h)=\text{E}\{\text{FDP}(h)\}.$$

$$\hat{h}=\inf \left\{h:0\le h\le (2\log q{)}^{1/2},\widehat{\text{FDP}}(h)\le \alpha \right\}.$$

$$\underset{(n_1,n_2,q)\to \infty }{\lim }\frac{\text{FDR}}{\alpha q_0/q}=1,\text{and}\frac{\text{FDP}}{\alpha q_0/q}\to 1\text{ in probability},\text{as}{n}_1,n_2,q\to \infty .$$

$$(i,j)\notin \mathcal{I}\text{ implies that }{s}_{1,i,j}-s_{2,i,j}=0,1\le i<j\le p.$$

$$A_{i,j}=\frac{{\bar{S}}_{1,i,j}+{\hat{\kappa }}_{i,j}{\bar{S}}_{2,i,j}}{(V_{1,i,j}/n_1+{\hat{\kappa }}_{i,j}^2{V}_{2,i,j}/n_2{)}^{1/2}},1\le i<j\le p,$$

$$a_{i,j}=\frac{s_{1,i,j}+{\kappa }_{i,j}{s}_{2,i,j}}{{\left\{\text{Var}(S_{1,l,i,j})+{\kappa }_{i,j}^2\text{Var}(S_{1,l,i,j})\right\}}^{1/2}},\text{ where }{\kappa }_{i,j}=\frac{n_2\text{Var}(S_{1,l,i,j})}{n_1\text{Var}(S_{2,l,i,j})}.$$

$${\text{P}}_{H_{0,i,j}}\left(|T_{i,j}|\ge h,|A_{i,j}|\ge \lambda \right)=\{1+o(1)\}G(h)\text{P}\left(|N(0,1)+a_{i,j}|\ge \lambda \right)+O(q^{-M}),$$

$${\text{P}}_{H_{0,i,j}}\left(|T_{i,j}|\ge h,|A_{i,j}|<{\lambda }_k\right)=\{1+o(1)\}G(h)\text{P}\left(|N(0,1)+a_{i,j}|<{\lambda }_k\right)+O(q^{-M}),$$

$$w_i={\left(\sum _{k=1}^K\frac{q_k{\hat{\pi }}_k}{1-{\hat{\pi }}_k}\right)}^{-1}\frac{q{\hat{\pi }}_k}{(1-{\hat{\pi }}_k)},i\in {\mathcal{G}}_k,1\le k\le K.$$

$$\mathcal{J}=\left\{(C_{1}N-1)\sqrt{\log q}/N,C_1\sqrt{\log q},\dots ,(C_{2}N-1)\sqrt{\log q}/N,C_2\sqrt{\log q}\right\}.$$

$${\text{FDP}}_{\text{adj}}=\frac{{\sum }_{i\in {\mathcal{H}}_0}I\left\{p_i^w\le p_{(\hat{\tau })}^w\right\}}{{\sum }_{i=1}^{q}I\left\{p_i^w\le p_{(\hat{\tau })}^w\right\}\vee 1},$$

$$\underset{(n_1,n_2,q)\to \infty }{\lim }\frac{{\text{FDR}}_{\text{adj}}}{\alpha q_0/q}=1,\text{and}\frac{{\text{FDP}}_{\text{adj}}}{\alpha q_0/q}\to 1\text{ in probability, as }{n}_1,n_2,q\to \infty .$$

$$\Psi =\text{E}\left\{\frac{{\sum }_{(i,j)\in {\mathcal{H}}_1}I(|T_{i,j}|\ge \hat{h})}{|{\mathcal{H}}_1|}\right\},{\Psi }_{\text{adj}}=\text{E}\left[\frac{{\sum }_{i\in {\mathcal{H}}_1}I\left\{p_i^w\le p_{(\hat{\tau })}^w\right\}}{|{\mathcal{H}}_1|}\right].$$

$${\Psi }_{adj}\ge \Psi +o(1),\text{ as }q\to \infty .$$

Full text

*Keywords: * Auxiliary information; False discovery rate; Multiple testing; Network data; Power enhancement.

1 Introduction

With prevalence of network data in recent years, the problem of comparing two populations of networks is gaining increasing attention. Our motivation is brain connectivity analysis, which studies functional and structural brain architectures through neurophysiological measures of brain activities and synchronizations (Varoquaux and Craddock, 2013). Accumulated evidences have suggested that, compared to a healthy brain, the brain connectivity network alters in the presence of numerous neurological disorders, for example, Alzheimer’s disease, autism spectrum disorder, among many others. Such alternations are believed to hold crucial insights of disease pathologies (Fox and Greicius, 2010). A typical brain connectivity study collects imaging scans, such as functional magnetic resonance imaging, or diffusion tensor imaging, from groups of subjects with and without disorder. Based on the imaging scan, a network is constructed for each individual subject, with the nodes corresponding to a common set of brain regions, and the edges encoding the functional or structural associations between the regions. A fundamental scientific question of interest is to compare the brain networks and to identify local connectivity patterns that alter between the two populations. Network comparison is equally interesting in many other scientific areas as well, for instance, clinical genomics, where of crucial interest is to understand and compare gene regulatory networks of patients with and without cancer (Luscombe et al., 2004).

In the context of brain connectivity analysis, there has been a rich literature on network estimation methods (Ahn et al., 2015; Qiu et al., 2016; Han et al., 2016; Wang et al., 2016; Zhu and Li, 2018, among many others). There is, however, a relative paucity of inference methods, especially simultaneous inference for individual links. Even though both can produce, in effect, a concise representation of the network structure, network inference is a fundamentally different problem than network estimation. Among the few existing network inference solutions, Kim et al. (2014) first summarized the network through a set of network metrics then employed a standard two-sample test. This strategy is commonly employed in the neuroscience literature and is easy to implement. However, it remains unclear to what extent each network metric provides a meaningful representation of brain function and structure (Fornito et al., 2013). Ginestet et al. (2017) characterized the geometry of the space of labeled, undirected networks with edge weights, established a central limit theorem for an appropriate notion of a network empirical mean, then developed an analog of the classical two-sample test. Kolaczyk et al. (2019) further extended this idea to unlabeled networks. However, these methods focused on the comparison of two entire networks and developed some global tests. None addressed the inference of individual links of two networks, nor considered simultaneous tests with false discovery control. Chen et al. (2015) developed a method to detect differentially expressed connectivity subnetworks under different clinical conditions by searching clusters of the graph. They resorted to a permutation test to obtain the $p$-value of the selected subnetwork, but did not provide the statistical significance of individual links and their differences. Besides, they only controlled the family-wise error rate, instead of the false discovery rate. Xia et al. (2015) first encoded the connectivity network by a partial correlation matrix computed from vector-valued data under a normal distribution. They then proposed a multiple testing procedure to compare the partial correlation matrices from the two populations, along with a proper false discovery control. Xia and Li (2019) further extended the test to matrix-valued data under a matrix normal distribution. In both cases, the test statistics were constructed based on the vector or matrix-valued data, which, as we explain next, may not be directly observable. Moreover, the underlying data distribution may not always be normal or matrix normal. Durante and Dunson (2018) developed a fully Bayesian solution for network comparison, which is very flexible and can handle the data format of our problem, but it requires specification of a series of prior distributions and can be computationally intensive.

Applications such as brain connectivity analysis actually raise new challenges for network inference. First, the observed data come in the form of $p\times p$ matrices, each of which encodes the network structure for one individual subject, and $p$ is the number of network nodes. For instance, in brain structural connectivity, what one observes are the numbers of white matter fibers between pairs of brain anatomical regions, and this matrix of counts forms a network, with brain regions constituting the nodes and the fiber counts the links. This is ultimately different from the data format studied in most network methods, where a network structure is inferred from some vector-valued or matrix-valued data and usually takes the form a Pearson correlation or partial correlation matrix. This fundamental difference in terms of the available data format would thus require a completely new problem formulation and inferential procedure. Second, in a multitude of applications including brain connectivity analysis, the sample size is usually very small, e.g., in tens. This calls for a testing procedure that is powerful enough to detect differentially expressed links under a limited sample size. In this article, we address the problem of comparing two populations of network data, more precisely, the two population means of networks. We aim to consider both global and simultaneous inferences, tackle the new data format, and explicitly enhance the power of the test.

Specifically, suppose we observe two groups of $p\times p$ network data, $\{\boldsymbol{S}_{1,1},\ldots,\boldsymbol{S}_{1,n_{1}}\}$ and $\{\boldsymbol{S}_{2,1},\ldots,\boldsymbol{S}_{2,n_{2}}\}$, where $n_{1}$, $n_{2}$ are the number of network samples for the two groups, respectively. Suppose $\boldsymbol{S}_{d,l}=(S_{d,l,i,j})_{p\times p}\sim\mathcal{F}_{d}(\boldsymbol{s}_{d})$, $l=1,\ldots,n_{d},d=1,2$, where $\mathcal{F}_{d}$ is some distribution with a symmetric mean matrix $\boldsymbol{s}_{d}=(s_{d,i,j})_{p\times p}$. Our goal is to test whether the two population means are the same; i.e.,

[TABLE]

If the global null in (1) is rejected, we further aim to identify at which locations the two mean matrices are different. That is, we wish to simultaneously test,

[TABLE]

In Xia et al. (2015), the observed data, $\boldsymbol{X}_{d,l}\in\mathbb{R}^{p}$, $l=1,\ldots,n_{d}$, $d=1,2$, is of the vector form, and is assumed to follow a normal distribution with the covariance matrix $\boldsymbol{\Sigma}_{d}$. Let $\boldsymbol{R}_{d}$ denote the corresponding partial correlation matrix, i.e., the standardized version of $\boldsymbol{\Sigma}_{d}^{-1}$, $d=1,2$. The network structure is then encoded by $\boldsymbol{R}_{d}$, and the problem becomes testing if $\boldsymbol{R}_{1}=\boldsymbol{R}_{2}$. Xia and Li (2019) followed a similar setup, except that the observed data $\boldsymbol{X}_{d,l}\in\mathbb{R}^{p\times t}$ becomes a matrix, and is assumed to follow a matrix normal distribution with the covariance $\boldsymbol{\Sigma}_{d}\otimes\boldsymbol{\Lambda}_{d}$, and the network is still encoded by the standardized version of $\boldsymbol{\Sigma}_{d}^{-1}$. The key difference for our setting is that, we do not always observe $\boldsymbol{X}_{d,l}$ directly, but instead $\boldsymbol{S}_{d,l}$ only. This difference in data format completely distinguishes our method from nearly all existing solutions such as Xia et al. (2015) and Xia and Li (2019). Moreover, we do not impose that the underlying data follows a normal or matrix normal distribution. Instead, we consider a general class of distributions for $\mathcal{F}_{d}$ satisfying some moment condition. Our method works for many different types of network links, for instance, binary links when $\mathcal{F}_{d}$ follows a light tailed distribution, or count links when $\mathcal{F}_{d}$ follows a heavy-tailed distribution.

For the global test (1), we develop a global test statistic taken as the maximum of a set of individual test statistics. We then derive its limiting null distribution, and show the resulting global test is power minimax optimal asymptotically. For the simultaneous test (2), we first develop a multiple testing procedure, and show that it can asymptotically control the false discovery at the pre-specified level. Next we propose a method to substantially enhance the power of the simultaneous inference procedure for (2). Specifically, we extend the grouping-adjusting-pooling idea of Xia et al. (2019a), and modify it for our inference of network data.

Our proposal differs from the existing solutions and makes several useful contributions. First, to the best of our knowledge, there has been no solution directly targeting simultaneous hypothesis testing of individual links for the network data in the format of $\boldsymbol{S}_{d,l}$. Our method bridges this gap, and offers a timely solution to a range of scientific applications where this form of problem and data is commonly encountered. Second, our global test statistic is constructed as the maximum of the individual test statistics for all links. This type of maximum statistic enjoys various advantages and has been commonly employed in the hypothesis testing literature (e.g., Cai et al., 2013; Cai and Zhang, 2016; Xia et al., 2019b). However, the derivation of its asymptotics, as well as the properties of the subsequent multiple testing procedure, are far from trivial in our new context of network comparison. Moreover, we remark that, in some network data applications, the individual test statistics may be correlated, and a global test statistic that utilizes such correlations may result in a more powerful test. However, this may not always be the case. For instance, in our brain connectivity application, the nodes are usually the brain anatomical regions, which can scatter at distant locations of the brain. As a result, there is no obvious correlation structure for the individual test statistics built on the pairs of brain regions. Therefore, we do not explicitly impose or employ any correlation structure when constructing the global test statistic. On the other hand, in our power enhancement procedure, we implicitly utilize the fact that some individual test statistics may be correlated and clustered. We then use a data driven approach to find such clusters and incorporate this information in our test. Finally, the power enhancement approach we develop is particularly useful in many applications, e.g., brain connectivity analysis, where the sample size is limited. Although motivated by Xia et al. (2019a), our enhancement method differs from Xia et al. (2019a) considerably in several ways. We explicitly compare the two power enhancement procedures in Section 4.5. Overall, we feel our method provides a useful addition to the general toolbox of network inference.

We adopt the following notation throughout this article. For a symmetric matrix $\boldsymbol{A}_{d}$, let $\lambda_{\max}(\boldsymbol{A}_{d})$ and $\lambda_{\min}(\boldsymbol{A}_{d})$ denote the largest and smallest eigenvalues of $\boldsymbol{A}_{d}$, respectively. For a set $\mathcal{H}$, let $|\mathcal{H}|$ denote its cardinality. For two sequences of real numbers $\{a_{n}\}$ and $\{b_{n}\}$, write $a_{n}=O(b_{n})$ if there exists a constant $C$ such that $|a_{n}|\leq C|b_{n}|$ holds for all $n$, write $a_{n}=o(b_{n})$ if $\lim_{n\rightarrow\infty}a_{n}/b_{n}=0$, and write $a_{n}\asymp b_{n}$ if there are positive constants $c$ and $C$ such that $c\leq a_{n}/b_{n}\leq C$ for all $n$. Write $n={n_{1}n_{2}}/{(n_{1}+n_{2})}$ and assume that $n_{1}\asymp n_{2}$.

The rest of the article is organized as follows. Section 2 presents the moment conditions for the distribution of $\mathcal{F}_{d}$ and show they are easily satisfied in numerous types of network data. Section 3 develops the global testing and the simultaneous testing procedures for the two-sample network comparison, and Section 4 studies power enhancement, both of which are key to our proposal. Section 5 presents the simulations, and a brain connectivity analysis as an illustration. The Supplementary Material collects additional lemmas and the proofs.

2 Moment Conditions and Examples

We begin with some moment conditions imposed on $\mathcal{F}_{d}$. We then give a number of examples and show that those conditions are easily satisfied in numerous types of network data.

2.1 Moment conditions

We assume that the distribution $\mathcal{F}_{d}$ of the network data $\boldsymbol{S}_{d,l}$ satisfies one of the following two conditions: a sub-Gaussian-type tail, or a polynomial-type tail, as stated below.

1. (C1)

(Sub-Gaussian-tail). Suppose that $\log p=o(n^{1/5})$, and that there exist some constants $\eta>0$ and $K>0$, such that, for $d=1,2$,

[TABLE] 2. (C2)

(Polynomial-tail). Suppose that $p\leq cn^{\gamma_{0}}$ for some constants $\gamma_{0}$ and $c>0$, and that there exist some constants $\epsilon>0$ and $K>0$, such that, for $d=1,2$,

[TABLE]

We first comment that, both conditions are common, and similar conditions have been often assumed in the high-dimensional setting (Cai et al., 2014; Van de Geer et al., 2014). These moment conditions are much weaker than the Gaussian assumption as usually required in the testing literature (Schott, 2007; Srivastava and Yanagihara, 2010). Next we discuss a number of network data examples that satisfy the above moment conditions, including Bernoulli and mixture Bernoulli data, Poisson data, as well as correlation and partial correlation data. Furthermore, we discuss some examples where the distributions are heavy-tailed, but after some data transformation, they still satisfy the moment conditions. Examples include transformed normal count data and transformed Wishart count data.

2.2 Binary network data

Binary network is arguably the most commonly seen network data type, where each link is a binary indicator. The Bernoulli distribution is often assumed; i.e., for $\boldsymbol{S}_{d,l}=(S_{d,l,i,j})_{p\times p}$, $S_{d,l,i,j}$ follows a Bernoulli distribution with mean $s_{d,i,j}$, $u<s_{d,i,j}<1-u$ for a constant $0<u<1$, $l=1,\ldots,n_{d},d=1,2$, and $1\leq i<j\leq p$. In such case, $\boldsymbol{S}_{d,l}$ satisfies the sub-Gaussian-tail condition in (C1), e.g., with $\eta=1$ and $K=(1-u)\exp\{u(1-u)^{-1}\}+u\exp\{(1-u)u^{-1}\}$.

The same holds true for the mixture Bernoulli distribution as discussed in Durante and Dunson (2018). That is, for some integer $H>0$ and randomly selected $\{\phi_{1},\ldots,\phi_{H}\}$ subject to $\sum_{h=1}^{H}\phi_{h}=1$ and $\phi_{h}>0$, $\textsf{P}(S_{d,l,i,j}=x)=\sum_{h=1}^{H}\phi_{h}\left\{s^{(h)}_{d,i,j}\right\}^{x}\left\{1-s^{(h)}_{d,i,j}\right\}^{1-x}$, with $u<s^{(h)}_{d,i,j}<1-u$ for some constant $0<u<1$, $x=0,1$, $h=1,\ldots,H$, $l=1,\ldots,n_{d},d=1,2$ and $1\leq i<j\leq p$. For this example, $\boldsymbol{S}_{d,l}$ again satisfies the sub-Gaussian-tail condition in (C1), with $\eta=1$ and $K=(1-u)\exp\{u(1-u)^{-1}\}+u\exp\{(1-u)u^{-1}\}$.

2.3 Correlation network data

Correlation network is another equally common network data type. In brain functional connectivity analysis and many other applications, the network is often encoded by a Pearson correlation or a partial correlation matrix. Take the Pearson correlation network as an example. The functional imaging data is usually summarized as a spatial-temporal matrix. That is, for the $l$th subject in the $d$th group, the observed data is of the form $\boldsymbol{X}_{d,l}\in\mathbb{R}^{p\times t_{d}}$, $l=1,\ldots,n_{d},d=1,2$, where $p$ is the number of brain regions, and $t_{d}$ is the number of repeated measures. Then the brain functional connectivity network is encoded by the sample correlation matrix $\boldsymbol{S}_{d,l}=t_{d}^{-1}\sum_{j=1}^{t_{d}}\{\boldsymbol{X}_{d,l,(\cdot,j)}-\bar{\boldsymbol{X}}_{d,l}\}\{\boldsymbol{X}_{d,l,(\cdot,j)}-\bar{\boldsymbol{X}}_{d,l}\}^{{\mathsf{\scriptscriptstyle T}}}$, where $\boldsymbol{X}_{d,l,(\cdot,j)}$ denotes the $j$th column of the matrix $\boldsymbol{X}_{d,l}$ and $\bar{\boldsymbol{X}}_{d,l}=t_{d}^{-1}\sum_{j=1}^{t_{d}}\boldsymbol{X}_{d,l,(\cdot,j)}$ denotes the sample mean vector (Fornito et al., 2013). Next we show that, as long as $\boldsymbol{X}_{d,l}$ satisfies one of the conditions in Lemma 1, then $\boldsymbol{S}_{d,l}$ satisfies the sub-Gaussian-tail condition (C1).

Lemma 1**.**

Suppose $\boldsymbol{X}_{d,l}$ satisfies one of the following conditions: (i) $\log p=o(t^{1/5})$, and there exist constants $\eta^{\prime}>0,K^{\prime}>0$ such that $\textsf{E}\left(\exp\left[\eta^{\prime}\{X_{d,l,i,j}-\textsf{E}(X_{d,l,i,j})\}^{2}/\textsf{Var}(X_{d,l,i,j})\right]\right)$ $\leq K^{\prime}$, where $t=\max\{t_{1},t_{2}\}$ and $t_{1}\asymp t_{2}$; (ii) $p\leq c^{\prime}t^{\gamma_{0}^{\prime}}$, for some $\gamma_{0}^{\prime},c^{\prime}>0$, and there exist constants $\epsilon^{\prime}>0,K^{\prime}>0$ such that $\textsf{E}\left[\left|\left\{X_{d,l,i,j}-\textsf{E}(X_{d,l,i,j})\right\}/\textsf{Var}(X_{d,l,i,j})^{1/2}\right|^{4\gamma_{0}^{\prime}+4+\epsilon^{\prime}}\right]\leq K^{\prime}$, for $i=1,\ldots,p,j=1,\ldots,t_{d}$. Then $\boldsymbol{S}_{d,l}$ satisfies the sub-Gaussian-tail condition in (C1), with $\eta=1/4$ and $K=2$, as $t\rightarrow\infty$.

We remark that a similar result as Lemma 1 can be obtained for the partial correlation network, by using the inverse regression techniques as in Liu (2013). Xia and Li (2019) tackled the network comparison problem assuming $\boldsymbol{X}_{d,l}$ is directly observable and follows a matrix normal distribution. Lemma 1 suggests that, the test we develop later is still applicable when $\boldsymbol{X}_{d,l}$ is available, even though it may not be as powerful as the test of Xia and Li (2019) in this case. On the other hand, the main focus of this article is to develop a test of comparing two networks even when $\boldsymbol{X}_{d,l}$ is not observed, but only $\boldsymbol{S}_{d,l}$ is. As such, our test is more general than that of Xia and Li (2019).

2.4 Poisson network data

Count network is another common network data type, where each link is a count. For instance, in brain structural connectivity analysis, the link is the number of white matter fibers between anatomical brain regions. The Poisson distribution is often imposed; i.e., $S_{d,l,i,j}$ follows a Poisson distribution with mean $s_{d,i,j}$, $0<u_{1}<s_{d,i,j}<u_{2}$, $l=1,\ldots,n_{d},d=1,2$, $1\leq i<j\leq p$. For any constant $\epsilon>0$, let $M$ be the smallest integer that is no smaller than $4\gamma_{0}+2+\epsilon$, where $\gamma_{0}$ is as defined in (C2). Then $\boldsymbol{S}_{d,l}$ satisfies the polynomial-tail condition (C2), with $K$ upper bounded by $u_{1}^{-(M-1)/2}\left[\sum_{i=0}^{M}u_{2}^{i}\genfrac{\{}{\}}{0.0pt}{}{M}{i}+u_{2}^{M}(u_{2}/2+1)\right]$, and $\genfrac{\{}{\}}{0.0pt}{}{M}{i}$ is the number of ways to partition a set of $M$ objects into $i$ non-empty subsets.

2.5 Transformed network data

We next consider some examples where the original network data $\boldsymbol{G}_{d,l}=(G_{d,l,i,j})_{p\times p}\sim\tilde{\mathcal{F}}_{d}(\tilde{\boldsymbol{s}}_{d})$, $l=1,\ldots,n_{d},d=1,2$, and $\tilde{\mathcal{F}}_{d}$ is some heavy-tailed distribution that only differs in the mean matrix $\tilde{\boldsymbol{s}}_{d}=(\tilde{s}_{d,i,j})\in\mathbb{R}^{p\times p}$ between the two groups. In such cases, testing the means of the original samples are equivalent to testing the means of the transformed data, $S_{d,l,i,j}=f(G_{d,l,i,j})$, where $f$ is some one-to-one transformation function. We next give two examples, where after transformation, the transformed data satisfies (C1) or (C2) or both.

One example is the log-normal count network. After the logarithmic transformation of $\boldsymbol{G}_{d,l}$, the transformed data $\boldsymbol{S}_{d,l}$ follows a normal distribution, and thus both (C1) and (C2) are satisfied. This can be further extended to the transformed normal mixture network.

Another example is the transformed Wishart count network, where the transformed data $\boldsymbol{S}_{d,l}$ follows the Wishart distribution with the scale matrix $m^{-1}\boldsymbol{\Sigma}_{d}$ and the degrees of freedom $m$. For this case, $\boldsymbol{S}_{d,l}$ satisfies the sub-Gaussian-tail condition (C1). Moreover, in this case, the testing problems (1) and (2) are closely related to the covariance matrix testing problems studied in Li and Chen (2012) and Cai et al. (2013). The key difference between our method and the existing ones is that, we only observe $\boldsymbol{S}_{d,l}$, but not the original vector samples. This example can also be further extended to the case of the product of Gaussian mixtures network, or the Wishart mixtures network.

3 Two-sample Test on Network Data

We begin with the construction of test statistic for the two testing problems (1) and (2). We then develop a global testing procedure for (1), and a simultaneous testing procedure for (2). For each test, we derive its corresponding asymptotic properties.

3.1 Test statistics

We first observe that the testing problem (1) is equivalent to the test, $H_{0}^{\prime}:\max_{1\leq i<j\leq p}|s_{1,i,j}-s_{2,i,j}|=0$. This motivates us to construct the test statistic based on

[TABLE]

where $\bar{S}_{d,i,j}=n_{d}^{-1}\sum_{l=1}^{n_{d}}S_{d,l,i,j}$. We standardize $W_{i,j}$, and estimate the variance of $S_{d,l,i,j}$ by

[TABLE]

respectively. This leads to our test statistic,

[TABLE]

3.2 Global test

In brain connectivity analysis and many other applications, it is generally postulated that the differences between two network structures concentrate on a small number of brain regions. This translates to a sparse alternative in our global test. Correspondingly, we construct the global test statistic as,

[TABLE]

Let $\boldsymbol{\Gamma}_{d}\in\mathbb{R}^{q\times q}$ denote the covariance matrix of vech($\boldsymbol{S}_{d,l}$), where $q=p(p-1)/2$, and vech$(\cdot)$ is the operator that turns the upper triangular part of $\boldsymbol{S}_{d,l}$ into a vector. Let $\boldsymbol{R}_{d}=(r_{d,i,j})\in\mathbb{R}^{q\times q}$ denote the corresponding correlation matrix. We introduce two conditions.

1. (A1)

$C_{0}^{-1}\leq\lambda_{\min}(\boldsymbol{\Gamma}_{d})\leq\lambda_{\max}(\boldsymbol{\Gamma}_{d})\leq C_{0}$ for some constant $C_{0}>0$, $d=1,2$. 2. (A2)

$\max_{d=1,2}\max_{1\leq i<j\leq q}|r_{d,i,j}|<r<1$ for some constant $0<r<1$.

Both conditions are mild. Particularly, Condition (A1) on the eigenvalues of the covariance matrix is common in the high dimensional setting (Bickel et al., 2008; Rothman et al., 2008; Yuan, 2010; Cai et al., 2014). Condition (A2) is also mild, because if $\max_{1\leq i<j\leq q}|r_{d,i,j}|=1$, then $\boldsymbol{\Gamma}_{d}$ is singular. We next obtain the limiting distribution of our test statistic $M_{n}$.

Theorem 1**.**

Suppose that (A1)-(A2), and one of (C1) and (C2) hold. Then for any $x\in\mathbb{R}$,

[TABLE]

Based on this limiting null distribution, we define the asymptotic $\alpha$-level test as,

[TABLE]

where $q_{\alpha}=-\log\pi-2\log\log(1-\alpha)^{-1}$.

We next study the power and the asymptotic optimality of the test $\Psi_{\alpha}$. Toward that end, define the sparsity of $\boldsymbol{s}_{1}-\boldsymbol{s}_{2}$ as $k_{q}=|\{(i,j):s_{1,i,j}-s_{2,i,j}\neq 0,1\leq i<j\leq p\}|$. We also introduce a class of $(\boldsymbol{s}_{1},\boldsymbol{s}_{2})$,

[TABLE]

Theorem 2**.**

Suppose that one of (C1) and (C2) holds. Then,

[TABLE]

Furthermore, suppose that $k_{q}=o(q^{r})$ for some $r<1/2$. Let $\alpha,\beta>0$ and $\alpha+\beta=1$. Then there exists a constant $c_{0}>0$ such that, for all sufficiently large $n_{d}$ and $q$,

[TABLE]

where $\mathcal{T}_{\alpha}$ is the set of all $\alpha$-level tests, i.e., $\textsf{P}_{H_{0}}({T}_{\alpha}=1)\leq\alpha$ for all $T_{\alpha}\in\mathcal{T}_{\alpha}$.

This theorem shows that the null hypothesis in (1) can be rejected by $\Psi_{\alpha}$ with a high probability if the pair of the network means belong to the class $\mathcal{U}(2\sqrt{2})$. In addition, with the mild sparsity condition $k_{q}=o(q^{r})$, the lower bound rate of $(\log q)^{1/2}$ cannot be further improved, because for a sufficiently small $c_{0}$, any $\alpha$-level test is unable to reject the null correctly uniformly over $\mathcal{U}(c_{0})$ with probability tending to 1. Henceforth, the global test $\Psi_{\alpha}$ reaches the power minimax optimality asymptotically.

3.3 Simultaneous test

We next develop a multiple testing procedure for (2) based on the test statistic $T_{i,j}$ in (4). Let $h$ be the threshold level such that $H_{0,i,j}$ is rejected if $|T_{i,j}|\geq h$. Let $\mathcal{H}_{0}=\{(i,j):s_{1,i,j}=s_{2,i,j},1\leq i<j\leq p\}$ be the set of true nulls, and $\mathcal{H}_{1}=\mathcal{H}\setminus\mathcal{H}_{0}$ the set of true alternatives, where $\mathcal{H}=\{(i,j):1\leq i<j\leq p\}$. Denote by $R_{0}(h)=\sum_{(i,j)\in\mathcal{H}_{0}}I(|T_{i,j}|\geq h)$ and $R(h)=\sum_{1\leq i<j\leq p}I(|T_{i,j}|\geq h)$ the total number of false positives and rejections, respectively. Then we define the false discovery proportion and false discovery rate by

[TABLE]

An ideal choice of $h$ would reject as many true positives as possible while controlling the FDP at the pre-specified level $\alpha$. That is, we select $h_{0}=\inf\left\{h\,:\,0\leq h\leq(2\log q)^{1/2},\,\text{FDP}(h)\leq\alpha\right\}$. Since $R_{0}(h)$ is unknown, we estimate it conservatively by $2q\{1-\Phi(h)\}$, where $\Phi(h)$ is the standard normal cumulative distribution function. This leads to our multiple testing procedure as summarized in Algorithm 1.

We next show that this testing procedure controls the FDR and FDP asymptotically at the pre-specified level. For notation simplicity, we write FDP=FDP($\hat{h}$) and FDR=FDR($\hat{h}$), where $\hat{h}$ is obtained in Algorithm 1. Define $\mathcal{A}_{i}(\xi)=\{j:\max(|r_{1,i,j}|,|r_{2,i,j}|)\geq(\log q)^{-2-\xi}\}$, and $\mathcal{S}_{\rho}=\{(i,j):1\leq i<j\leq p,|s_{1,i,j}-s_{2,i,j}|/\{\textsf{Var}(S_{1,l,i,j})/n_{1}+\textsf{Var}(S_{2,l,i,j})/n_{2}\}^{1/2}\geq(\log q)^{1/2+\rho}\}$. We further introduce some conditions.

1. (B1)

$|\mathcal{S}_{\rho}|\geq[{1}/\{\pi^{1/2}\alpha\}+\delta]({\log q})^{1/2}$, for some constant $\delta>0$ and any sufficiently small constant $\rho>0$. 2. (B2)

$\max_{1\leq i\leq q}|\mathcal{A}_{i}(\xi)|=o(q^{\nu})$ for some constants $\xi>0$ and $0<\nu<(1-r)/(1+r)$. 3. (B3)

$q_{0}=|\mathcal{H}_{0}|\geq c_{1}q$ for some constant $c_{1}>0$.

Condition (B1) on $\mathcal{S}_{\rho}$ is mild, as it only requires a small number of $\boldsymbol{s}_{1}$ and $\boldsymbol{s}_{2}$ having standardized difference with the order of $(\log q)^{1/2+\rho}$ for any sufficiently small constant $\rho>0$. Condition (B2) is mild, as it requires that not too many $S_{d,l,i,j}$ are highly correlated, but still allows the number of highly correlated pairs to grow in the order of $o(q^{1+\nu})$. Condition (B3) is also a natural and mild assumption, because if it does not hold, i.e., $q_{0}=o(q)$, then we can simply reject all the hypotheses. As a result, we would have $|R_{0}|=q_{0}$, $|R|=q$, and the FDR would tend to zero. Under these conditions, we obtain the asymptotic properties of our multiple testing procedure in terms of false discovery control.

Theorem 3**.**

Suppose that (A2), (B1)-(B3), and one of (C1) and (C2) hold. Then,

[TABLE]

4 Power Enhancement

In brain connectivity analysis and many other applications, the sample size $n_{d}$ is often small, whereas the number of nodes $p$ can be moderate to large. This results in a limited power for the proposed testing procedure. We explore in this section an explicit power enhancement method that has potential to substantially improve the power of the simultaneous inference developed in Section 3.3. We borrow the idea of grouping, adjusting and pooling (GAP) that was first proposed in Xia et al. (2019a). However, our method differs from Xia et al. (2019a) in many ways, including a different, and actually less restrictive, assumption, a different set of primary and auxiliary statistics, and a different modification of the multiple testing procedure. We show that the modified procedure is asymptotically more powerful, meanwhile it can still control FDR and FDP asymptotically. We obtain these properties assuming the sub-Gaussian-tail condition (C1). Parallel results can be obtained under the polynomial-tail condition (C2) as well, but are technically more involved. We begin by describing the intuition behind our power enhancement solution, then derive the proper auxiliary statistic for our inference problem. We then develop the modified simultaneous testing procedure, and study its asymptotic properties in terms of power improvement and false discovery control. We also compare in detail our method with the GAP method of Xia et al. (2019a).

4.1 Intuition

We recognize that there exists additional information in the data that is potentially useful to improve the simultaneous testing procedure of Algorithm 1. We first discuss our intuition, then use some simple example to illustrate where the auxiliary information is and how it can facilitate our multiple testing procedure.

In a multitude of applications including brain connectivity analysis, it is often believed that the difference between the two networks under different biological conditions is small. This means $\boldsymbol{s}_{1}-\boldsymbol{s}_{2}$ is sparse. Accordingly, one can find a baseline matrix $\boldsymbol{s}_{0}$, such that $\boldsymbol{s}_{1}^{\prime}=\boldsymbol{s}_{1}-\boldsymbol{s}_{0}$ and $\boldsymbol{s}_{2}^{\prime}=\boldsymbol{s}_{2}-\boldsymbol{s}_{0}$ are individually sparse. Let $\mathcal{I}_{d}=\{(i,j):s_{d,i,j}^{\prime}\neq 0,1\leq i<j\leq p\}$ denote the support of $\boldsymbol{s}_{d}^{\prime}$, $d=1,2$, and $\mathcal{I}=\mathcal{I}_{1}\cup\mathcal{I}_{2}$ denote the union support. Note that the set of alternative hypotheses $\mathcal{H}_{1}$ defined in Section 3.3 is the same as $\mathcal{I}$, if $s_{1,i,j}\neq s_{2,i,j}$ for every $(i,j)\in\mathcal{I}_{1}\cap\mathcal{I}_{2}$. In general, $\mathcal{H}_{1}$ is a proper subset of $\mathcal{I}$. Since $\boldsymbol{s}_{1}^{\prime}$ and $\boldsymbol{s}_{2}^{\prime}$ are both sparse, we realize that the cardinality of $\mathcal{I}$ is small. Moreover, the following relationship holds true:

[TABLE]

Therefore, the knowledge about $\mathcal{I}$ is useful to help narrow down the search in multiple testing. In other words, if one can find a way to identify possible entries $(i,j)$ in $\mathcal{I}$, it would provide useful information about the set of true alternatives $\mathcal{H}_{1}$, or equivalently, the set of true nulls $\mathcal{H}_{0}$. As a consequence, it can potentially increase the power of the testing procedure.

A key observation is then, while the test statistic is built on the difference between $\bar{S}_{1,i,j}$ and $\bar{S}_{2,i,j}$ as defined in Section 3.1, the sum of $\bar{S}_{1,i,j}$ and $\bar{S}_{2,i,j}$ can provide crucial information about $\mathcal{I}$. Consider a toy example where the network data is binary, and $S_{d,l,i,j}$ follows a Bernoulli distribution with mean $s_{d,l,i,j}$, $l=1,\ldots,n_{d},d=1,2,1\leq i<j\leq p$. Assume that $s_{1,i,j}=s_{2,i,j}=s_{0,i,j}=0.1$ for $80\%$ of the $(i,j)$ pairs, $s_{1,i,j}=s_{2,i,j}=s_{0,i,j}=0.9$ for $10\%$ of the $(i,j)$ pairs, and for the rest of the $(i,j)$ pairs, $s_{1,i,j},s_{2,i,j}\sim\text{Uniform}(0.1,0.9)$ and $s_{0,i,j}=0.1$. In this example, for the pairs $(i,j)\notin\mathcal{I}$, the sum of $s_{1,i,j}$ and $s_{2,i,j}$ is either very small, which is 0.2, or very large, i.e., 1.8. Meanwhile, for the pairs $(i,j)\in\mathcal{I}$, the sum is in between. Henceforth, this sum contains useful information about $\mathcal{I}$, and can potentially enhance the power of the multiple testing procedure.

Based on the above discussion, we can see that, the more sparsity structure information the auxiliary statistics can capture, the more information they can provide about the union support $\mathcal{I}$, and the more substantial power gain the test can achieve. In general, the sparser the true difference $\boldsymbol{s}_{1}-\boldsymbol{s}_{2}$ is, the more information the auxiliary statistics can offer.

4.2 Auxiliary statistics

We next formally construct the auxiliary statistic that provides useful information about the union support $\mathcal{I}$. It is important to note that, the auxiliary statistic should be constructed so that they are asymptotically independent of the test statistic $T_{i,j}$ in (4). This way the null distribution of $T_{i,j}$ would not be distorted by the incorporation of the auxiliary statistic.

Recall that $V_{d,i,j}$ is the sample variance of $S_{d,l,i,j}$ as defined in (3). We construct the auxiliary statistic $A_{i,j}$ as,

[TABLE]

where $\hat{\kappa}_{i,j}=(n_{2}V_{1,i,j})/(n_{1}V_{2,i,j})$. The next proposition shows that the test statistic $T_{i,j}$ and the auxiliary statistic $A_{i,j}$ are asymptotically independent under the null hypothesis. Define

[TABLE]

Proposition 1**.**

Suppose (C1) holds. For any constants $M>0$ and $C>0$, we have

[TABLE]

uniformly for $0\leq h\leq C\sqrt{\log q}$, $0\leq\lambda\leq C\sqrt{\log q}$, and $1\leq i<j\leq p$, with $G(h)=2\{1-\Phi(h)\}$. Furthermore, for all $0\leq k\leq CN$ with an integer constant $N$,

[TABLE]

uniformly for $0\leq h\leq C\sqrt{\log q}$ and $1\leq i<j\leq p$, where $\lambda_{k}=(k/N)\sqrt{\log q}$.

4.3 Power enhanced simultaneous test

Based on $(T_{i,j},A_{i,j})$, we now modify the simultaneous testing procedure of Algorithm 1. We first describe the main idea. We next summarize the modified testing procedure in Algorithm 2. Finally, we discuss some specific choices of the key parameters of the algorithm.

Since there are totally $q=p(p-1)/2$ tests to carry out simultaneously, we rearrange the pairs of $\{(T_{i,j},A_{i,j}),1\leq i<j\leq p\}$ into $\{(T_{i},A_{i}),i=1,\ldots,q\}$. After obtaining all the $p$-values, $p_{i}=2\{1-\Phi(|T_{i}|)\}$, from Algorithm 1, our basic idea is to adjust those $p$-values by $p_{i}^{w}=\min\{p_{i}/w_{i},1\}$, with $w_{i}$ being the adjusting weights, $i=1,\ldots,q$. We utilize the auxiliary statistics $A_{i}$ to help compute the adjusting weights $w_{i}$, by groups. Specifically, we consider a set of grid values, $\mathcal{J}=\{(C_{1}N-1)\sqrt{\log q}/N,C_{1}\sqrt{\log q},\ldots,(C_{2}N-1)\sqrt{\log q}/N,C_{2}\sqrt{\log q}\}$, where $C_{1}$, $C_{2}$ and $N$ are some pre-specified constants. We divide the index set $\{1,\ldots,q\}$ into $K$ groups according to the auxiliary statistics $(A_{1},\ldots,A_{q})$. As an example, we take $K=3$. That is, we choose two grid points $\mathcal{J}_{K}=\{\lambda_{1},\lambda_{2}\}$ in $\mathcal{J}$, and obtain $K=3$ groups of indices, $\mathcal{G}_{1}=\{i\,:\,1\leq i\leq q,\,-\infty<A_{i}\leq\lambda_{1}\}$, $\mathcal{G}_{2}=\{i\,:\,1\leq i\leq q,\,\lambda_{1}<A_{i}\leq\lambda_{2}\}$, and $\mathcal{G}_{3}=\{i\,:\,1\leq i\leq q,\,\lambda_{2}<A_{i}\leq\infty\}$. For each group $\mathcal{G}_{k}$, we compute its cardinality, $q_{k}=|\mathcal{G}_{k}|$. We also estimate the proportion, $\pi_{k}$, of alternatives in $\mathcal{G}_{k}$, $k=1,\ldots,K$. To do so, we employ the method of Schweder and Spjøtvoll (1982) and Storey (2002) to obtain an estimate $\tilde{\pi}_{k}$ first, then stablize it by $\hat{\pi}_{k}=(\epsilon\vee\tilde{\pi}_{k})\wedge(1-\epsilon)$, where $\epsilon$ is a small positive number; we set $\epsilon=10^{-5}$. Then for all the indices in $\mathcal{G}_{k}$, we compute the group-wise adjusting weight:

[TABLE]

This idea of adjusting the weights $w_{i}$ by groups is motivated by our intuition discussed in Section 4.1. After obtaining the weights, we adjust the $p$-values and apply the Benjamini-Hochberg procedure (Benjamini and Hochberg, 1995, BH) to the adjusted $p$-values $p_{i}^{w}$. Finally, we search all possible choices of $\mathcal{J}_{K}$ among $\mathcal{J}$, and find the one that yields the largest number of rejections. We apply the BH procedure again to the adjusted $p$-values under this choice of $\mathcal{J}_{K}$ to obtain the final adjusted rejection region. We summarize this modified simultaneous testing procedure in Algorithm 2.

We discuss some specific choices of the parameters in Algorithm 2. First, the number of groups $K$ is usually set at $K=3$. As shown in Xia et al. (2019a), when $K\geq 4$, there is little additional power gain, but a more expensive computation. Second, the constants $C_{1}$ and $C_{2}$ can be chosen so that $C_{1}\sqrt{\log q}$ is equal to the smallest value of the auxiliary statistics and $C_{2}\sqrt{\log q}$ is equal to the largest value of the auxiliary statistics. If the absolute values of the smallest and largest auxiliary statistics exceed $16\sqrt{\log q}$, we truncate at $C_{1}=-16$ and $C_{2}=16$ to stabilize and expedite the computation. We note here that, if the network data are non-negative, such as the binary and poisson network data, then both $C_{1}$ and $C_{2}$ are non-negative. By contrast, in Xia et al. (2019a), $C_{1}$ and $C_{2}$ were fixed at $-4$ and $4$. Finally, $N$ can be any integer for the theoretical validity. Numerically, a larger value of $N$ implies a more precise grid search, but at the cost of a heavier computational burden. We choose $N$ such that the gap between two adjacent grid points, $(\log p)^{1/2}/N$, equals $0.1$ approximately.

4.4 FDR control and power enhancement

We next show that the modified inference of Algorithm 2 is asymptotically more powerful than Algorithm 1, meanwhile it can still asymptotically control the false discovery.

Denote $\{p_{i}^{w}:1\leq i\leq q\}$ the adjusted $p$-values from Algorithm 2, and $\{p_{(i)}^{w}:1\leq i\leq q\}$ the ordered adjusted $p$-values. The corresponding adjusted FDP is:

[TABLE]

where $\hat{\tau}$ is the cutoff obtained from Step 3.4 of Algorithm 2, and $I(\cdot)$ is the indicator function. Accordingly, $\textrm{FDR}_{\textrm{adj}}=\textsf{E}(\textrm{FDP}_{\textrm{adj}})$. The next theorem shows that the modified procedure can still control FDR and FDP asymptotically.

Theorem 4**.**

Suppose (A2), (B1)-(B3), and (C1) hold. Suppose $\log q=o(n^{1/C})$ for some constant $C>5$. Then,

[TABLE]

Next, denote the power of the testing procedures of Algorithms 1 and 2 by $\Psi$ and $\Psi_{\textrm{adj}}$, respectively. That is,

[TABLE]

Then the next theorem shows that, by incorporating the auxiliary statistics $A_{i,j}$, the modified simultaneous testing procedure of Algorithm 2 is asymptotically more powerful than Algorithm 1, which is solely based on the test statistics $T_{i,j}$.

Theorem 5**.**

Suppose the same conditions in Theorem 4 hold. Then,

[TABLE]

4.5 Comparison to GAP

Although motivated by the GAP method of Xia et al. (2019a), our power enhancement procedure is also considerably different from GAP. While Xia et al. (2019a) tackled the problem of mean comparison of vector-valued samples, we target the problem of network mean comparison. This leads to a different set of test and auxiliary statistics, but a number of additional intrinsic differences as well.

First, the two methods impose different assumptions. A key requirement for GAP to enhance the power is that the parameters of interest from each group are individually sparse. In our setup, however, the parameters may all be non-negative. For instance, in a binary network or a count network, all the entries of both means $\boldsymbol{s}_{1}$ and $\boldsymbol{s}_{2}$ are usually non-negative. As such, the means may not be individually sparse. Our procedure instead only requires the difference of the two means $\boldsymbol{s}_{1}-\boldsymbol{s}_{2}$ is sparse, which reasonably holds and is often imposed in numerous applications including brain connectivity analysis (Zhu and Li, 2018).

Second, the two methods differ in terms of the range of the auxiliary statistics that contribute most to the power enhancement. Consider the case when $K=3$. In Xia et al. (2019a), since both means are assumed to be individually sparse, the tests that are more likely to be adjusted and rejected are those with the corresponding auxiliary statistics either being negative and small, or positive and large. That is, the power enhancement hinges more on those tests in $\mathcal{G}_{1}$ and $\mathcal{G}_{3}$ with small or large auxiliary statistics. However, in our setup, the individual means $\boldsymbol{s}_{1}$ and $\boldsymbol{s}_{2}$ can both be dense and their entries are all positive. Instead we only assume that $\boldsymbol{s}_{1}-\boldsymbol{s}_{2}$ is sparse. Take a binary brain connectivity network as an example. The observed networks are often sparse, in that most links are zero, since it is known that brain connections are energy consuming and biological units tend to minimize energy-consuming activities (Raichle and Gusnard, 2002; Bullmore and Sporns, 2009). This translates to small connection probabilities for most entries of $\boldsymbol{s}_{1}$ and $\boldsymbol{s}_{2}$, while all these probabilities are positive. Moreover, the difference of the means between the two populations is often sparse, which translates to equal connection probabilities for most entries of $\boldsymbol{s}_{1}$ and $\boldsymbol{s}_{2}$, or equivalently zero difference for most entries of of $\boldsymbol{s}_{1}-\boldsymbol{s}_{2}$. This is similar to the toy example we discuss in Section 4.1. For such cases, as a consequence of Algorithm 2, the tests whose corresponding auxiliary statistics are too small or too large would be adjusted so that they are less likely to be rejected. Instead, those tests whose auxiliary statistics are in between would be adjusted so that they are more likely to be rejected. In other words, the power enhancement in our setup may hinge more on $\mathcal{G}_{2}$, rather than $\mathcal{G}_{1}$ and $\mathcal{G}_{3}$.

Third, due to the above difference, the grid construction in Step 1.4 of Algorithm 2 is noticeably different from that of GAP in Xia et al. (2019a). Specifically, in Xia et al. (2019a), to ensure the inclusion of important locations in $\mathcal{G}_{1}$ and $\mathcal{G}_{3}$, the constants $C_{1}$ and $C_{2}$ can be simply fixed at $-4$ and $4$, respectively, so that the upper bound of those small negative auxiliary statistics and lower bound of those large positive auxiliary statistics can be attained in the grid $\mathcal{J}$. By contrast, for our problem, the upper bound of the auxiliary statistics in the union support $\mathcal{I}$ can go beyond the bound in Xia et al. (2019a), i.e., $4\sqrt{\log q}$, and the lower bound of the auxiliary statistics in the union support can be non-negative. Since the tests in $\mathcal{G}_{2}$ are more likely to be adjusted and rejected, we need to do a more thorough grid construction and choose the constants $C_{1}$ and $C_{2}$ based on the smallest and largest values of the auxiliary statistics as described in Section 4.3.

5 Numerical Analysis

We first study and compare the finite-sample performance of the two simultaneous inference procedures, Algorithms 1 and 2, through simulations. We then illustrate our method with a brain structural connectivity analysis.

5.1 Simulations

We consider a $68\times 68$ network, which is of the same dimension as the real data in Section 5.2. This results in $q=68(68-1)/2=2278$ hypotheses to test simultaneously. We consider three network data distributions, three sparsity levels, and two different sample sizes. More specifically, we consider Bernoulli, Bernoulli mixture, and transformed Wishart distributions.

• •

**Bernoulli: ** Select the sets $\mathcal{M}_{d,1}$ and $\mathcal{M}_{0}$ from $q$ hypotheses, uniformly and randomly, with $|\mathcal{M}_{d,1}|=|\mathcal{M}_{0}|=k_{q}/2$, $d=1,2$. Here $k_{q}$ is a parameter that controls the sparsity level, and is defined later. Let $\mathcal{M}_{d}=\mathcal{M}_{d,1}\cup\mathcal{M}_{0}$. For $(i,j)\notin\mathcal{M}_{d}$, generate $S_{d,l,i,j}\sim\text{Bernoulli}(1,0.3)$. For $(i,j)\in\mathcal{M}_{d}$, generate $S_{d,l,i,j}\sim\text{Bernoulli}(1,r_{d,i,j})$, where $r_{1,i,j}$ is equal to 0.5 with probability 0.1, and is equal to 0.8 otherwise, whereas $r_{2,i,j}$ is equal to 0.8 with probability 0.1, and is equal to 0.5 otherwise.

• •

**Bernoulli mixture: ** Generate $\mathcal{M}_{d}$ in the same way as before. Generate $S_{d,l,i,j}\sim\text{Bernoulli}(1,r_{d,i,j})$, where $r_{d,i,j}=\pi_{i,j}*r_{d,1,i,j}+(1-\pi_{i,j})*r_{d,2,i,j}$, with $\pi_{i,j}\sim\text{Uniform}(0,1)$. For $(i,j)\notin\mathcal{M}_{d}$, $r_{d,1,i,j}=r_{d,2,i,j}=0.3$, $d=1,2$. For $(i,j)\in\mathcal{M}_{d}$, $r_{1,1,i,j}$ is equal to 0.5 with probability 0.1, and is equal to 0.7 otherwise, whereas $r_{2,1,i,j}$ is equal to 0.7 with probability 0.1, and is equal to 0.5 otherwise, and $r_{d,2,i,j}=r_{d,1,i,j}+0.2$.

• •

**Wishart with logarithm transformation: ** Select the sets $\mathcal{M}_{d,1}$ and $\mathcal{M}_{0}$ from $q$ hypotheses, uniformly and randomly, with $|\mathcal{M}_{d,1}|=k_{q}/4$, and $|\mathcal{M}_{0}|=3k_{q}/4$, $d=1,2$. Let $\mathcal{M}_{d}=\mathcal{M}_{d,1}\cup\mathcal{M}_{0}$. Generate $\Sigma_{d}$ such that $\Sigma^{\prime}_{d,i,j}=\text{Uniform}(3,5)$ if $(i,j)\in\mathcal{M}_{d}$ and $\Sigma^{\prime}_{d,i,j}=0$ otherwise. Let $\Sigma^{\prime}_{d,j,i}=\Sigma^{\prime}_{d,i,j}$ and $\Sigma_{d}=\Sigma^{\prime}_{d}+\{|\lambda_{\min}(\Sigma^{\prime}_{d})|+0.5\}I$, where $I$ is the identify matrix. Generate $S_{d,l}^{\prime}\sim\text{Wishart}(m^{-1}\Sigma_{d},m)$, with $m=100$, and $S_{d,l}=\log[\text{round}\{\exp(S_{d,l}^{\prime})\}]$, where round$(\cdot)$ rounds a number to the nearest integer.

The sparsity level between the two network population means is controlled by the parameter $k_{q}$, and we consider three sparsity levels, $k_{q}=0.2q,0.15q$ and $0.1q$. We also examine two sample sizes, $n_{1}=n_{2}=100$ and $n_{1}=n_{2}=25$, and the latter mimics the real data setting where the sample size is very limited.

We apply both Algorithms 1 and 2, and set the nominal level at $\alpha=0.05$. We report the empirical FDR and power, both in percentage, based on 100 replications in Table 1. It is seen that, in all cases, the empirical FDRs are generally controlled under the nominal level by both algorithms. Algorithm 2 is slightly more conservative than Algorithm 1, which is mainly due to the normalization step of the weight calculation as shown in (5). A similar phenomenon has also been observed in Xia et al. (2019a). For the empirical power, it is seen that Algorithm 2 has a clear power improvement over Algorithm 1, thanks to its utilization of the auxiliary information. Furthermore, the performance under the varying sample size confirms the power enhancement of Algorithm 2 as theoretically revealed in Section 4.4. We also observe that the power gain becomes more substantial when the true difference $\boldsymbol{s}_{1}-\boldsymbol{s}_{2}$ becomes more sparse, which agrees with our intuition explained in Section 4.1.

5.2 Structural connectivity analysis

We illustrate our method with a brain structural connectivity analysis of the KKI-42 dataset, which is available at http://openconnecto.me/data/public/MR/archive/, and its detailed description can be found in Landman et al. (2011). The data consists of 21 subjects with no history of neurological conditions, aging from 22 to 61 years old. Each subject received two resting-state diffusion tensor images (DTI) under a scan-rescan imaging session. For simplicity, we treat the data as if those images were from independent samples, which is common for the analysis of this dataset (Wang et al., 2017). It results in a total sample size of $42$ for this study. Brain regions are constructed following the Desikan Atlas (Desikan et al., 2006), leading to $p=68$ regions equally divided in the left and right hemispheres. DTI is a magnetic resonance imaging technique that measures the diffusion of water molecules to map white matter tractography in the brain. Each DTI image was preprocessed, and was summarized in the form of a $68\times 68$ network, where the edges record the total number of white matter fibers between the pair of nodes. It is also equally common to focus on the form of a binary network, where the edges become the binary indicators of presence of absence of white matter fibers (Wang et al., 2017). See Zhang et al. (2018) for more details on the construction of a brain structural network from DTI images. We partition the subjects into two age groups, the ones whose are younger than 30 years, and the ones who are 30 or older. Age 30 is a transition period, usually known as the “age 30 transition”, when the first phase of early adulthood comes to a close, and the basis for the next life structure is formed. Moreover, this partition yields about the same number of subjects for each group, with $n_{1}=22$ for the younger-than-30 age group, and $n_{2}=20$ for the older-than-30 age group. We study the age-related difference in structural connectivity patterns, which is of universal interest, as aging is the main risk factor for progressive loss of either structure or function of brain neurons (Morrison and Hof, 1997).

We applied the power-enhanced test to this dataset, first the binary network, then the count network with a logarithm transformation. We set the significance level at 0.05. For the binary network, out of the total of 2278 links, Algorithm 1 identified 38 significantly different links, whereas the power-enhanced Algorithm 2 identified 62 links, including all 38 linked found by Algorithm 1 plus 24 additional links. This agrees with both our theory and the simulations, in that Algorithm 2 is usually able to recognize more significant links than Algorithm 1. Figure 1 plots all the identified links by Algorithm 2 and the associated brain regions, visualized using the BrainNet Viewer (Xia et al., 2013). Table 2 records the 24 additional links found by Algorithm 2. For the count network, out of the total of 2278 links, Algorithm 1 identified 140 significantly different links, whereas Algorithm 2 identified 177 links. Among the 140 links found by Algorithm 1, 125 links were also found by Algorithm 2. Moreover, Algorithm 2 identified 52 additional links. Figure 2 shows all the identified links by Algorithm 2. Those additional links found by Algorithm 2 capture the smaller changes that were missed by Algorithm 1. The biological roles of those links are not immediately clear, but they point to potentially interesting connectivity patterns that require further scientific validation.

6 Supplementary Material

The additional lemmas and theorem proofs are available in the online supplementary material.