On the testing of multiple hypothesis in sliced inverse regression

In the general framework of the regression analysis, the goal is to infer the relation between a univariate response variable $y$ and a $p\times 1$ vector $\bm{x}$. One would like to know $y|\bm{x}$, namely, how the distribution of $y$ depends on the value of $\bm{x}$. Among the literature of sufficient dimension reduction [25, 8, 9, 28, 29], the fundamental idea is to replace the predictor by its projection to a subspace without loss of information. In other words, we seek for a subspace $\mathcal{S}_{y|\bm{x}}$ of the predictor space such that

$$y\rotatebox[origin={c}]{90.0}{$\models$}\bm{x}|\mathbf{P}_{\mathcal{S}_{y|\bm{x}}}\bm{x}.$$

(1)

Here $\models$ indicates independence, and $P_{(\cdot)}$ stands for a projection operator. The subspace $\mathcal{S}_{y|\bm{x}}$ is called the central subspace. Let $d$ be the dimension of this central subspace. Let $\mathbf{B}$, a $p\times d$ matrix, be a basis of the central subspace $\mathcal{S}_{y|\bm{x}}$. Then the equation (1) is equivalent to

$$y\rotatebox[origin={c}]{90.0}{$\models$}\bm{x}|\mathbf{B}\bm{x}.$$

(2)

To further reduce the dimensionality especially when the number of predictors $p$ diverges with respect to $n$, it is commonly assumed that $y$ depends on $\bm{x}$ through a subset of $\bm{x}$, known as the Markov blanket and denoted as $\mathcal{MB}(y,\bm{x})$ [33, 36, 5], such that

$$y\rotatebox[origin={c}]{90.0}{$\models$}\bm{x}|\mathcal{MB}(y,\bm{x}).$$

For each predictor, one would like to know whether $x_{j}\in\mathcal{MB}(y,\bm{x})$, which could be formulated as a multiple testing problem. The null hypothesis stating that $x_{j}\notin\mathcal{MB}(y,\bm{x})$ is equivalent to

$$\mathcal{H}_{j}:\mathbf{P}_{span(x_{j})}(\mathcal{S}_{y|\bm{x}})=\mathcal{O}_{p},$$

(3)

where $\mathcal{O}_{p}$ is the origin point in the $p$-dimensional space [10]. In other words, it is equivalent to saying that the $j$-th row of the matrix $\mathbf{B}$ consists of all zeros.

There are many attempts targeting estimating the central subspace in the existing literature on sufficient dimension reduction. The most widely used method is the sliced inverse regression (SIR) which was first introduced in [25]. Later, there are many attempts to extend SIR, including, but not limited to, sliced average variance estimation [14, 8], directional regression [24], and constructive estimation [41]. Nevertheless, most of the existing methods and theories in sufficient dimension reduction focus on the estimation of the central subspace $\mathcal{S}_{y|\bm{x}}$. The result of statistical inference is very limited when $p$ diverges, not to mention the procedure of controlling the false discovery rate (FDR) when testing these hypotheses simultaneously.

The challenge arises from two perspectives. First, among the literature on sufficient dimension reduction, the result on the limiting distribution of the estimator of the central subspace is very limited when $p$ diverges. When $p$ is fixed, [20, 46, 30, 23] have derived the asymptotic distribution of the sliced inverse regression. To the best of the authors’ knowledge, there are no results on the limiting distribution when $p$ diverges unless assuming the signal is strong and the total number of false hypotheses is fixed [40].

Second, after the test statistic is determined for each hypothesis, how to combine these test statistics to derive a method that controls the false discovery rate is challenging. Many existing procedures, such as [2, 3, 42, 35], work on the (generalized) linear regression models. In [5], the authors considered an arbitrary joint distribution of $y$ and $\bm{x}$ and proposed the model-X knockoff to control FDR. However, this method requires that the distribution of the design matrix is known, which is not feasible in many practice.

The study conducted by [21] explored variable selection for the linear regression model, assuming the condition of weak and rare signals. It is noted in this paper that the selection consistency is not possible and allowing for false positives is necessary. While several existing penalization-based methods on sufficient dimension reduction require imposing uniform signal strength condition to achieve consistent results ([27, 29, 38, 44]), this work tackles a more challenging scenario. Specifically, we develop the central limit theorem of SIR, utilizing recent theories on Gaussian approximation [6, 7] without relying on the uniform signal strength conditions. This theoretical result is the first of its kind in the literature on sufficient dimension reduction and is a necessity for simultaneous inference when the effects of some relevant predictors are either moderate or weak.

We proceed by constructing a statistic for each hypothesis based on sliced inverse regression. Applying Gaussian approximation theory, we demonstrate that this statistic is asymptotically symmetric about zero when the null hypothesis holds. We refer to this as a Angular Balanced Statistic (ABS). We then develop a single-step procedure that rejects a hypothesis when its ABS exceeds a certain threshold. Additionally, we provide an estimator for the false discovery proportion. For a designated FDR level $q$, we adaptively select a threshold such that the estimated false discovery proportion is no greater than $q$. This method is referred to as the Model-free multiple Testing procedure using Angular balanced statistics (MTA). Theoretical analysis confirms that MTA asymptotically controls the FDR at the $q$-level under regularity conditions. Simulation results and data analysis demonstrate that MTA significantly outperforms its competitors in terms of power while controlling the FDR.

The paper is organized as follows. In Section 2, we derive the central limit theorem of SIR when both the dimension $p$ and the number of important predictors $s$ diverge with $n$ using the recently developed Gaussian approximation theory. In Section 3, we construct ABS based on the estimator of SIR and propose the MTA method. It is shown that the FDR of the MTA method is less than or equal to a designated level asymptotically. In Sections 4 and 5, we provide numerical evidence including extensive simulations and a real data analysis to demonstrate advantages of MTA. We conclude the paper in Section 6, and include all technical details in the appendix.

Recall that $y$ is the response and $\bm{x}$ is a p-dimensional vector. In the literature of sufficient dimension reduction, people aim to find the central subspace $\mathcal{S}_{y|\bm{x}}$ defined in (1). The sliced inverse regression (SIR) introduced in [25] is the first and the most popular method among many existing ones.

Assume that the covariance matrix of $\bm{x}$ is $\bm{\Sigma}$. Let $\bm{\Omega}=\bm{\Sigma}^{-1}$ be the precision matrix of $\bm{x}$. Let $s_{j}=\#\{k:\bm{\Omega}_{jk}\neq 0\}$ and $s=\max_{j}s_{j}$. When the distribution of $\bm{x}$ is elliptically symmetric, it is shown in [25] that

$\displaystyle\bm{\Sigma}\mathcal{S}_{y|\bm{x}}=col(\bm{\Lambda}),$

(4)

where $\bm{\Lambda}=var(\mathbb{E}[\bm{x}|y])$ and $col(\bm{\Lambda})$ is the column space spanned by $\bm{\Lambda}$.

Given $n$ $i.i.d.$ samples $(y_{i},\bm{x}_{i})$, $i=1,\cdots,n$. To estimate $\bm{\Lambda}$, divide the data into $H$ slices according to the order statistics $y_{(i)}$, $i=1,\ldots,n$. Let $\bm{x}_{(i)}$ be the concomitant associated with $y_{(i)}$. Note that slicing the data naturally forms a partition of the support of the response variable, denoted as $\mathcal{H}$. Let $\mathcal{P}_{h}$ be the $h$-th slice in the partition $\mathcal{H}$. Here we let $\mathcal{P}_{1}=(-\infty,y_{(\lceil n/H\rceil)}]$ and $\mathcal{P}_{H}=(y_{(\lceil n/H\rceil*(H-1)+1)},+\infty)$. Let $\bar{\bm{x}}$ be the mean of all the $\bm{x}$’s and $\bar{\bm{x}}_{h,\cdot}$ be the sample mean of the vectors $\bm{x}_{(j)}$’s such that its concomitant $y_{(j)}\in\mathcal{P}_{h}$ and estimate $\bm{\Lambda}\triangleq var(\mathbb{E}[\bm{x}|y])$ by

$$\widehat{\bm{\Lambda}}_{\mathcal{H}}=\frac{1}{H}\sum_{h=1}^{H}(\bar{\bm{x}}_{h,\cdot}-\bar{\bm{x}})(\bar{\bm{x}}_{h,\cdot}-\bar{\bm{x}})^{\tau}.$$

(5)

The $\widehat{\bm{\Lambda}}_{\mathcal{H}}$ was shown to be a consistent estimator of $\bm{\Lambda}$ under some technical conditions [18, 20, 46, 25, 28].

Alternatively, we could view SIR through a sequence of ordinary least squares regressions. Let $f_{h}(y),h=1,2,\cdots,H$ be a sequence of transformations of $Y$. Following the proof of [43, 40], one knows that under the linearity condition [25],

$$Ef(y)\phi(y)\in\mathcal{S}_{y|\bm{x}},$$

where $\phi(y)=\bm{\Sigma}^{-1}E(\bm{x}|y)$. Let $\bm{\beta}_{h}\in R^{p\times 1}$ be defined as

$$\bm{\beta}_{h}=argmin_{\bm{\beta}_{h}}E(f_{h}(y)-\bm{x}\bm{\beta}_{h})^{2}.$$

Assuming the coverage condition [32, 13, 40], then

$$Span(\mathbf{B})=\mathcal{S}_{y|\bm{x}},$$

where $\mathbf{B}=(\bm{\beta}_{1},\cdots,\bm{\beta}_{H})\in R^{p\times H}$.

Note that different choices of $f_{h}(y)$ lead to different methods [40, 17]. To name a few here, [43] suggested $f_{h}(y)=y^{h}$ where $h\leq H$. After slicing the data into $H$ slices according to the value of the response variable $y$, [13] suggested $f_{h}(y)=y$ if $y$ is in the $h$-th slice and 0 otherwise. If we choose $f_{h}(y)=\mathbbm{1}(y\in\mathcal{P}_{h})$, this will lead to SIR, which is the main focus in this paper [43, 40].

After obtaining data $(\bm{x}_{i},y_{i})$ based on a sample of $n$ subjects, let

$$f_{h}(\bm{y})=\mathbbm{1}(\bm{y}\in\mathcal{P}_{h})=\left(\mathbbm{1}(y_{1}\in\mathcal{P}_{h}),\mathbbm{1}(y_{2}\in\mathcal{P}_{h}),\cdots,\mathbbm{1}(y_{n}\in\mathcal{P}_{h})\right)^{T}.$$

Let $\hat{\bm{\beta}}_{h}$ be defined as

$$\hat{\bm{\beta}}_{h}=argmin_{\bm{\beta}_{h}}||f_{h}(\bm{y})-\bm{x}^{T}\bm{\beta}_{h}||^{2}=(\bm{x}\bm{x}^{T})^{-1}\bm{x}f_{h}(\bm{y}),$$

or a general form

$$\hat{\bm{\beta}}_{h}=argmin_{\bm{\beta}_{h}}||f_{h}(\bm{y})-\bm{x}^{T}\bm{\beta}_{h}||^{2}=\frac{1}{n}\hat{\bm{\Omega}}\bm{x}f_{h}(\bm{y}),$$

(6)

where $\hat{\bm{\Omega}}$ is a suitable approximation of the inverse of the Gram matrix $\hat{\bm{\Sigma}}=\bm{x}^{T}\bm{x}/n$. Let

$$\hat{\bm{B}}=(\hat{\bm{\beta}}_{1},\cdots,\hat{\bm{\beta}}_{H}).$$

(7)

There are many methods to estimate the precision matrix. As an example, we could consider the one given by the Lasso for the node-wide regression on the design matrix $\bm{x}$ [31]. Next, we will derive the central limit theorem of SIR when $p$ converges to the infinity.

Our derivation is built upon the Gaussian approximation (GAR) theory recently developed in [6, 7]. Let $\mathcal{P}_{H}=\cup_{h}\mathcal{P}_{h}$ be a partition of the sample space of $y$ and $p_{h}=P(y\in\mathcal{P}_{h})$. Define

$$\tilde{\bm{\beta}}_{h}=\frac{1}{n}\bm{\Omega}\bm{x}f_{h}(\bm{y}).$$

(8)

Let $\tilde{\bm{B}}=(\tilde{\bm{\beta}}_{1},\cdots,\tilde{\bm{\beta}}_{H})$. For $i=1,2,\cdots,n,j=1,2,\cdots,p,h=1,2,\cdots,H$, let $z_{ijh}$’s be normal random variables such that

• •

  $Ez_{ijh}=p_{h}\bm{\Omega}_{\cdot j}^{T}E(\bm{x}|y\in\mathcal{P}_{h})$;

• •

  $V(z_{ijh})=p_{h}\bm{\Omega}_{\cdot j}^{T}E(\bm{x}\bm{x}^{T}|y\in\mathcal{P}_{h})\bm{\Omega}_{\cdot j}-p_{h}^{2}(\bm{\Omega}_{\cdot j}^{T}E(\bm{x}|y\in\mathcal{P}_{h}))^{2}$;

• •

  $Cov(z_{ijh},z_{ikh})=p_{h}\bm{\Omega}_{\cdot j}^{T}E(\bm{x}\bm{x}^{T}|y\in\mathcal{P}_{h})\bm{\Omega}_{\cdot k}-p_{h}^{2}\bm{\Omega}_{\cdot j}^{T}E(\bm{x}|y\in\mathcal{P}_{h})\bm{\Omega}_{\cdot k}^{T}\\ E(\bm{x}|y\in\mathcal{P}_{h})$;

• •

  $Cov(z_{ijh_{1}},z_{ikh_{2}})=-p_{h_{1}}p_{h_{2}}\bm{\Omega}_{\cdot j}^{T}E(\bm{x}|y\in\mathcal{P}_{h_{1}})\bm{\Omega}_{\cdot k}^{T}E(\bm{x}|y\in\mathcal{P}_{h_{2}})$.

Note that when $\mathcal{H}_{j}$ is true and the coverage condition holds, then $Ez_{ijh}=0$ for any $i$ and $h$ (see Remark LABEL:remark:z in the appendix).

Let $\bm{z}_{i\cdot\cdot}$ be a $p$ by $H$ matrix consisting of $z_{ijh}$ where $j=1,2,\cdots,p$ and $h=1,2,\cdots,H$. Let $\mathcal{A}^{re}$ consist of all sets $A$ of the form

$$A=\{\omega\in R^{pH}:a_{j}\leq\omega_{j}\leq b_{j},j=1,2,\cdots,pH\}$$

(9)

for some $-\infty\leq a_{j}\leq b_{j}\leq\infty$. We will develop a bound of the quantity

$$\rho_{n}(\mathcal{A}^{re})=\sup_{A\in\mathcal{A}^{re}}\left|\mathbb{P}(\sqrt{n}\hat{\mathbf{B}}\in A)-\mathbb{P}(\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\bm{z}_{i\cdot\cdot}\in A)\right|.$$

(10)

The theoretical investigation requires the following assumptions

Assumption 1 is assumed in [45], describing the tail property of the predictors. The bounds of the eigenvalues of $\bm{\Sigma}$ are commonly assumed when considering the estimation of the covariance matrix and the precision matrix. In the sliced inverse regression, the data are sliced according to the response variable $y$. This creates a partition on the support of the response variable. We assume the following condition on the partition.

The condition (D1) on the partition requires that the probability data falls in each slice is inversely proportional to the number of slices. This condition is assumed in [28] when establishing the phase transition of SIR. The condition (D2) requires a sufficient variation of the $\bm{x}$, which falls in each slice. Otherwise, if the smallest eigenvalue of that covariance matrix in a given slice converges to zero, the corresponding estimator becomes unstable. In Remark 2, we give examples that (D2) holds.

The proof is put in the appendix. From this theorem, it is seen that the estimator $\hat{\mathbf{B}}$ could be approximated by Gaussian random vectors which preserve the mean and covariance structure.

In this section, we consider testing the hypothesis defined in (3) simultaneously such that the false discovery rate is controlled at a designated level $q$. We start with the construction of the test statistic.

First, we randomly split the sample as two parts, denoted as $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ [15]. Let $\hat{\mathbf{B}}^{1}$ and $\hat{\mathbf{B}}^{2}$ be the estimators of the sliced inverse regression based on these two subsets respectively. We define the following mirror statistics

$$M_{j}=\sum_{h=1}^{H}\hat{\mathbf{B}}^{1}_{jh}\hat{\mathbf{B}}^{2}_{jh}.$$

(12)

Note that when the j-th hypothesis is true, namely, $\mathbf{P}_{span(x_{j})}(\mathcal{S}_{y|\bm{x}})=\mathcal{O}_{p}$, one would expect that both vectors $\hat{\mathbf{B}}_{j\cdot}^{1}$ and $\hat{\mathbf{B}}_{j\cdot}^{2}$ would center around zero. On the other hand, when $H_{j}$ is false, both vectors center around a common non-zero vector $\bm{B}_{j\cdot}$. The angles between $\hat{\bm{B}}_{j\cdot}^{k},k=1,2,$ and $\bm{B}_{j\cdot}$ are approximately the same. Consequently, we call this statistic $M_{j}$ a Angular Balanced Statistic (ABS). Intuitively, the ABS statistic would be symmetric with respect to zero under the null and it tends to have a positive value under the alternative. This is depicted in Figure 1. In the left panel, when considering the predictor $x_{1}\in\mathcal{MB}(y,\bm{x})$, the estimated coefficient vectors $\hat{\mathbf{B}}_{1\cdot}^{1}$ and $\hat{\mathbf{B}}_{1\cdot}^{2}$ center around a non-zero vector. On the other hand when considering the estimated coefficients of $x_{2}$ which does not belong to $\mathcal{MB}(y,\bm{x})$, the estimated coefficients $\hat{\mathbf{B}}_{2\cdot}^{1}$ and $\hat{\mathbf{B}}_{2\cdot}^{2}$ center around the origin. In the right panel, we plot the histogram of all the mirror statistics when generating the data according to Setting 1 in Section 4, it is seen that $M_{j}$’s are roughly symmetric with respect to zero when the corresponding null hypotheses are true. When the corresponding null hypothesis is false, this statistic tends to have a large positive value.

The following theorem provides a rigorous argument for the above intuition. Namely, we show that the ABS is symmetric with respect to zero under the null hypothesis asymptotically.

We now turn to the false discovery rate control. To simplify our notation, let $\theta_{j}$ be the indicator of whether the $j$-th hypothesis defined in (3) is false. Let $\mathcal{S}_{0}=\{j:\theta_{j}=0\}$ be the set consisting of all the true hypotheses. Consider the single-step procedure where the $j$-th hypothesis is rejected if the ABS $M_{j}$ is greater than some threshold $t$. Namely, we choose a decision $\bm{\delta}=(\delta_{1},\cdots,\delta_{p})$ as $\delta_{j}=1(M_{j}>t)$. Then the false discovery proportion is

$$FDP(t)=\frac{\sum_{j:\theta_{j}\in\mathcal{S}_{0}}\delta_{j}}{(\sum_{j}\delta_{j})\vee 1}.$$

According to Theorem 2, the probabilities $\mathbb{P}(M_{j}>t)$ and $\mathbb{P}(M_{j}<-t)$ are approximately the same when the hypothesis $\mathcal{H}_{j}$ is true. Namely,

	$\displaystyle\mathbb{E}\sum_{j:\theta_{j}\in\mathcal{S}_{0}}\delta_{j}=\sum_{j:\theta_{j}\in\mathcal{S}_{0}}\mathbb{P}(M_{j}>t)$
	$\displaystyle\approx\sum_{j:\theta_{j}\in\mathcal{S}_{0}}\mathbb{P}(M_{j}<-t)\leq\sum_{j}\mathbb{P}(M_{j}<-t).$

Consequently, we estimate the false discovery proportion by

$$\widehat{FDP(t)}=\frac{\sum_{j}1(M_{j}<-t)}{\sum_{j}1(M_{j}>t)}.$$

(14)

Note that $\widehat{FDP(t)}$ tends to overestimate the true false discovery proportion. However, this is not a serious issue due to two reasons: (i) the over-estimation protects us from large false positives and (ii) the probability that $\mathbb{P}(M_{j}<-t)$ is negligible when $\mathcal{H}_{j}$ is false.

For any designated level, say $q$, we could choose a data-driven threshold $\tau_{q}$ as

$$\tau_{q}=\min\left\{t:\widehat{FDP(t)}\leq q\right\}.$$

We call this method Model free multiple Testing procedure based on ABS (MTA). The steps of MTA are summarized in Algorithm 1.

We now turn our attention to the asymptotic property of MTA. We assume the following assumption:

For linear model, a similar assumption stating that the sum of partial correlation is bounded by $C_{1}p_{0}^{\alpha}$, where $\alpha\in(0,2)$ is imposed in [42]. When dealing with nonlinear model, Assumption 3 requires that the sum of pair-wise partial correlation in each slice being bounded by $Cp_{0}^{\alpha}$ for a constant $C>0$ any $\alpha\in(0,2)$. However, we consider a more general non-linear model where the contribution of the $j$-th variable is related to vectors with length $H$. Assumption 3(a) is equivalent to bound the pairwise canonical correlations which is a general procedure for investigating the relationships between two sets of variables. In Assumption 3(b), we consider the absolute sum of the correlation coefficients of corresponding variables in two sets. Since the correlation measure is not unique, we consider two different types to enlarge the suitable cases for our theory. The following two lemmas characterize the sum of pairwise covariances of mirror statistics, which is the key to establishing the asymptotic property of MTA.

In Lemma 1, we build a relationship between the correlations among mirror statistics and the correlations among covariates. Next, we derive asymptotic properties of the MTA method. We define a few quantities:

	$$\displaystyle V(t)=\frac{\#\{j:j\in\mathcal{S}_{0},M_{j}\leq-t\}}{p_{0}},\,V^{\prime}(t)=\frac{\#\{j:j\in\mathcal{S}_{0},M_{j}\geq t\}}{p_{0}},\,$$
	$$\displaystyle V^{1}(t)=\frac{\#\{j:j\in\mathcal{S}_{1},M_{j}\geq t\}}{p_{1}},$$
	$$\displaystyle G_{0}(t)=\lim_{p}\frac{\sum_{j\in\mathcal{S}_{0}}\mathbb{E}1(M_{j}\leq-t)}{p_{0}},\,G_{1}(t)=\lim_{p}\frac{\sum_{j\in\mathcal{S}_{1}}\mathbb{E}1(M_{j}\geq t)}{p_{1}},$$

and

	$\displaystyle FDP(t):=\frac{V^{\prime}(t)}{(V^{\prime}(t)+r_{p}V^{1}(t))\vee 1/p},$
	$\displaystyle FDP^{\infty}(t):=\lim_{p}\frac{V^{\prime}(t)}{(V^{\prime}(t)+r_{p}V^{1}(t))\vee 1/p},$

where $r_{p}=p_{0}/p_{1}$ and $FDP^{\infty}(t)$ is the pointwise limit of $FDP(t)$. We have the following results.

Lemma 2 is based on the weak dependence assumption 3(a) or 3(b) under a model-free assumption. In literature, similar week dependence conditions on p-values [37] and linear models [15, 42, 16] are commonly used to study the asymptotical property of FDR. With the aid of this lemma, we show that the FDR of the proposed MTA method is less than or equal to $q$ asymptotically.

The proof of Theorem 3 is included in the appendix.

In this section, we use simulations to evaluate the finite sample performance of MTA with state-of-the-art competitors. First, we consider the model-X knockoff, which constructs knockoffs $\tilde{\bm{x}}$ by assuming the distribution of $\bm{x}$ is Gaussian with a known covariance matrix. We construct the feature statistic as

$$W_{j}=\sum_{h=1}^{H}\hat{\mathbf{B}}_{jh}\hat{\mathbf{B}}_{\tilde{j}h},$$

where $\hat{\mathbf{B}}$ is the solution of (7) when replacing $\bm{x}$ by $(\bm{x},\tilde{\bm{x}})$. Second, we consider the marginal non-parametric test to test the dependence between the response variable and each predictor and calculate the Hilbert-Schmidt independence criterion (HSIC) test statistic proposed in [19]. After obtaining the p-values, we apply the Benjamini-Hochberg (BH) procedure (HSIC+BH). Third, we apply the dimension reduction coordinate test (DR-coor) in Cook [10] to test the contribution of each predictor and calculate the p-value followed by the BH method (DR-coor+BH). In the following sections, we consider three different multi-index models with a wide range of nonlinear dependencies between the response and predictors. For predictors, we consider two different designs including Gaussian design and a real-data design.

Neural imaging studies hold great promise for predicting and reducing psychiatric disease burden and advancing our understanding of the cognitive abilities that underlie humanity’s intellectual feats. Functional magnetic resonance imaging (fMRI) data are increasingly being used for the ambitious task of relating individual differences in studying brain function to typical variations in complex psychological phenotypes. In functional brain imaging research, e.g., via fMRI, functional network analysis, which focuses on understanding the connections between different brain areas using some quantitative methods, has become a booming area popular approach. The human connectome project (HCP) dataset aims at characterizing human brain connectivity in about 1,200 healthy adults to enable detailed studies of the association between brain functions and behaviors at the level of individual subjects. Customized scanners were used to produce high-quality and consistent data to measure brain connectivity. The data containing various traits and fMRI data can be easily accessed through (db.humanconnectome.org)

To obtain functional connectomes, we used a state-of-the-art fMRI data preprocessing framework – population-based functional connectome (PSC) mapping. The dataset used in this work comes from the HCP Q1 release [39] consisting of fMRI scans from 68 subjects during resting state (rfMRI). The acquisition parameters are as follows: 90×104 matrix, 220mm FOV, 72 slices, TR=0.72s, TE=33.1ms, flip angle=52°. Preprocessing steps of the rfMRI data have been conducted within the Q1 release, including motion correction, spatial smoothing, temporal pre-whitening, slice time correction, and global drift removal. The fMRI data from each individual in each task was then linearly registered to the MNI152 template.

Since we are particularly interested in cognition, we extract four cognition-related measures as y from HCP, including:

• •

  Oral reading recognition test: participants on this test are asked to read and pronounce letters

• •

  Crystallized composite score: crystallized cognition composite can be interpreted as a global assessment of verbal reasoning. We use the age-corrected standard score.

Our predictors are brain functional connections between each pair of ROIs. Based on the brain parcellation, we have $68*68/2$ connections. We first applied a screening procedure to select top $1/5$ of the connections as our candidates. We applied the proposed MTA method to select important connections associated with the cognitive traits with controlled false discovery rate at $10\%$.

As shown in the left panel of Figure 4, we observe dense connections within both the left and right frontal lobes. In particular, we see that brain regions such as r23 (rostral middle frontal), and l27 (caudal middle frontal) are densely involved (nodes with high degrees) in the selected connections. We also observe that the four nodes in the temporal lobe (regions 8, 14, 29), all appear in the selected connections to the frontal lobe. The frontal lobe processes how you interpret language. The temporal lobes are also believed to play an important role in processing auditory information and with the encoding of memory [34]. Hence, our results show that richer connections between the visual and auditory are strongly associated with higher learning ability in vocabulary. As shown in the right panel of Figure 4, we observe more connections between the front lobe and the temporal lobe. Temporal lobe regions are linked to memory, language, and visual perception [1]. Additionally, we observe a high level of communication between the left and right brain hemispheres, which contributes to higher cognition [4].

In this paper, we consider the multiple hypothesis testing procedure which aims at controlling FDR without assuming a model of the conditional distribution of the response. We combine the idea of data-splitting, Gaussian mirror, and sliced inverse regression (SIR) to construct mirror statistics. With the aid of the developed central limit theorem of SIR, it is shown that the mirror statistic is approximately symmetric with respect to zero when the null hypothesis is true. We then provide an estimator of the false discovery proportion for the single-step testing method and suggest a data-driven threshold based on such an estimator. It is theoretically shown that the FDR of the proposed method is less than or equal to a designated level asymptotically.

The theoretical investigation of the SIR using the newly developed Gaussian approximation theory paves the road for the statistical inference for sufficient dimension reduction with a diverging dimension. It has the potential to be used for other sufficient dimension reduction methods such as sliced average variance estimator (SAVE, [8]), principle Hessian directions (PHD, [26]), and many others. We will leave them in future discussions.

The proposed approach is flexible due to the lean assumptions on the link function. The popular deep neural network could be viewed as the multiple index model. In an ongoing project, it would be interesting to see how the proposed method could be used to reduce the dimension under this framework.

Data used in the preparation of this article were obtained from the Human Connectome Project (www.humanconnectome.org). The HCP WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) were funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. We would like to thank Bing Li, Wen Zhou, Asaf Weinstein for their comments and suggestions. Xing’s research is supported by NSF (DMS-2124535).