Controlling the False Discovery Rate in Complex Multi-Way Classified Hypotheses

In modern inferential problems arising from large datasets, hypotheses often occur in complex structures of groups. Compelling examples can be found in neuro-imaging studies like functional magnetic resonance imaging (fMRI) or electroencephalography (EEG). Large sets of hypotheses obtained from such studies can be classified according to (a) locations of interest (voxels in an fMRI study and electrodes in an EEG study) and (b) their broader classifications (clusters of voxels known as Regions of Interest (ROIs) and regions of the cerebral cortex of the brain). The groups formed by the locations of interest can be nested within the groups formed by their broader classifications. In a different scenario, a set of hypotheses can be classified by more than one independent criteria that unlike the previous case, cannot be ordered among themselves. For example, Stein et al. (2010) conducted a study to detect voxelwise genomewide association in the human brain. The study involved the analysis of the relationship between $448293$ single-nucleotide polymorphisms(SNPs) and volume change in $31622$ voxels. The hypotheses obtained in this study can be classified according to groups formed by their SNPs and/or groups formed by a family for each voxel. Furthermore, in such studies, when the grouping criteria are defined by the nature of the scientific experiment and statisticians have little or no control on their demarcation, it is quite likely the groups overlap with one or more hypotheses becoming members of more than one group. For example, in an EEG experiment, some electrodes are placed at the boundaries of two adjacent brain regions, thereupon, hypotheses that are associated to these electrodes can be classified as members of either brain regions.

Though multiple hypotheses testing is a highly active research area, besides the p-filter algorithm proposed by Ramdas et al. (2019), not much development has been done to adequately justify overlapping groups, simultaneous or nested group structures in newly designed testing procedures. All the same, the grouping information provide valuable insight into the nature of the hypotheses which can be profitably used to create sharper distinction of signals from noise. The hypotheses classified into a group share common characteristics with other member hypotheses of the group. Consequently, a grouping structure levies a generic effect on its members, and a multiple testing method that utilizes all such valuable information in its algorithm would be naturally poised to make sharper distinction between the significant and non-significant state of an individual hypothesis. Some researchers have addressed the task of developing new multiple testing procedures adjusted to simple classification structure where hypotheses are classified into distinct groups by one norm of classification (Hu et al. (2010), Benjamini and Bogomolov (2014), Ignatiadis et al. (2016), Li and Barber (2019), etc.). However, to the best of our knowledge, there does not exist a well-defined multiple testing procedure that controls the overall False Discovery Rate, by utilizing structural information gleaned from multiple classifications, that are either simultaneous, or hierarchical or both, and make provisions for overlapping groups. In the absence of such a procedure, an investigator has to prioritize a classification over the rest and choose a multiple testing procedure that accommodates non-overlapping groups formed by a single classification norm. In a situation where each of the multiple classifications contain equally important information, the investigator is obligated to compromise the accuracy of the analysis.

In this article we propose the theory of a generalized testing procedure, and its data-adaptive version, in an effort to alleviate this problem. We systematically develop a weighted Benjamini-Hochberg (BH)-type procedure, the Generalized Grouped BH (Gen-GBH) procedure from combining two weighted BH-type procedures, the Heir-GBH and the $S$-way GBH, designed for hierarchical classification and simultaneous classification structures respectively. Our proposed methods control the overall False Discovery Rate (FDR) in their respective classification settings. The theoretical formulations of the methods control FDR when the p-values satisfy some positive dependence criteria, while their data adaptive analogs control FDR under assumptions of independence. Due to the detailed classification information used in our approach, our proposed methods are poised to control FDR more accurately and possess more power than existing procedures that utilize partial or no classification information.

The paper is arranged as follows. In Section 2, we present a brief review of multiple testing research closely related to ours and some preliminary results at the base of our proposed methods. In Section 3, we present the p-value weighting schemes for hypotheses structured according to three different types of classification, hierarchical, multi-way simultaneous, and generalized classifications. The corresponding weighted BH methods in their oracle forms with proven FDR control under positive dependence and our proposals of their data-adaptive versions are also given in this section. The generalized classification integrates hierarchical classification into a simultaneous setting. Simulation studies in Section 4 investigate the performances of the oracle and data-adaptive weighted BH for hierarchically grouped hypotheses relative to their BH counterparts that ignore the underlying structure. Though the p-filter algorithm (Ramdas et al. (2019)) considers a general classification setup such as ours, it serves FDR controlling objectives that are quite different from the objectives in this paper. Nevertheless, we include it in our comparative studies, noting the caveats caused by the differences in our goals. The results of these studies indicate that the p-value weighting scheme in our proposed weighted BH under hierarchical classification setting can capture the underlying structure quite effectively. This carries over to our proposed data-adaptive method under the generalized classification, as seen from its application to an EEG dataset, illustrated in Section 5. The paper concludes with Section 6 where some remarks are made on the effectiveness of multiple testing as an important statistical tool for analyzing complexly structured datasets. Proofs of some results associated with the FDR control of our proposed methods are given in the Appendix.

An example demonstrating the advantage of the Heir-GBH over the BH and respectively, their data-adaptive versions, in a hierarchical classification setup. Through a short example, we illustrate the superior distinguishing power of the Heir-GBH and the DAHeir-GBH over existing comparable testing procedures, the oracle BH and its data-adaptive version. We consider the setup similar to that in the first row of Figure 1, consisting of $N=25$ hypotheses that are hierarchically classified in two levels, with the p-values corresponding to significant hypotheses highlighted. The Adaptive BH disregards classification structure and assigns to each p-value $P_{i},i=1,\ldots,N$, the constant estimate of weight (see Storey et al. (2004)):

where $R_{N}(\lambda)=\sum_{i=1}^{N}\mathds{1}(P_{i}\leq\lambda)$, with $\lambda$ being set at $0.5$, for the Adaptive BH as well as for the DAHeir-GBH. In Figure 2, the layout of the p-values is shown, alongwith a table showing the weights assigned by the four methods in three stages. At level $0$, since there is no grouping, the hierarchical methods (with choice of weights as in (3.2), (3.3), and (3.7), (3.10)) are equivalent to the BH in its oracle and data-adaptive forms, respectively. As seen from the table of weights, the Adaptive BH, equivalent to the DAHeir-GBH at $L=0$, assigns constant weight $0.56$, while the oracle BH, equivalent to the Heir-GBH at $L=0$ assigns the weight $0.6$ to all the p-values. At level $L=1$, the weights assigned by the Heir-GBH and DAHeir-GBH are more sensitive to the significant state of the hypotheses than the oracle BH or the Adaptive BH. The five hypotheses at the intersection of the two groups at level one are assigned the smallest weights, appropriately accounting for the effects of both overlapping parent groups. At level $L=2$, the hierarchical procedures are more sensitive to the significant states of the p-values; the highest weights are assigned to the members of group 4 at level two, that does not contain any non-null p-values, and the smallest weights are assigned to the members of group 3 at level two, that are in the intersection of the overlapping groups at level one and contain the largest number of non-null p-values. The weights are sensitive to the effect of grouping on the significance of the individual p-value, and accuracy of the weights increase with the details of the hierarchical structure.

We further investigate the performance of the Heir-GBH in comparison with the oracle BH and the p-filter algorithm. The DAHeir-GBH is compared to the Adaptive BH and the pfilter algorithm with adaptive weights. The comparisons are made, respectively, in terms of (i) control on FDR and (ii) average power at varying densities of signals. We consider $N=5000$ hypotheses in a hierarchical setup of two levels. The layout resembles a section of the hierarchical setup of the EEG dataset we analyzed in the following secion. At the first level, there are two groups, i.e., $m_{1}=2$, the first group $G(g_{1}=1)$ comprising of $n_{g_{1}=1}=3000$ hypotheses and the second group $G(g_{1}=2)$ consisting of $n_{g_{1}=2}=2500$ member hypotheses. The groups overlap and there are $500$ hypotheses at the intersection of $G(g_{1}=1)$ and $G(g_{1}=2)$. At level 2, each of these two groups are further classified into smaller groups each containing $100$ members, i.e., $m_{2}(G(g_{1}=1))=30$, $n_{g_{1}g_{2}}=n_{1g_{2}}=100,\forall g_{2}=1,\ldots,m_{2}(G(g_{1}=1))$ and $m_{2}(G(g_{1}=2))=25$, $n_{g_{1}g_{2}}=n_{2g_{2}}=100,\forall g_{2}=1,\ldots,m_{2}(G(g_{1}=2))$. For the sake of simplification, the hypotheses at the intersection of the two groups at level one are classified into $5$ groups, each of size $100$ at level two. Thus at level $L=2$, there are $50$ groups, each of size $100$.

The simulation setting was designed using the following steps:

1. 1.

   Generate the following arrays to simulate the impact of classification at each level of grouping.

   1. (a)

      $\bm{\Theta}_{L=0}$ as an $m\times n$ random matrix of i.i.d. Ber$(1-\pi_{0})$, where $m=50,n=100$. The elements of the matrix $\Theta_{L=0}$ represent the state of the individual hypotheses at level $L=0$.

   2. (b)

      $\bm{\Theta}_{L=1}$ as defined below demonstrates the influence of second level of grouping on the set of hypotheses.

      $\displaystyle\bm{\Theta}_{L=1}=\begin{pmatrix}\theta_{G(g_{1}=1)}\\ \\ \theta_{G(g_{1}=1)\cap G(g_{1}=2)}\\ \\ \theta_{G(g_{1}=2)}\end{pmatrix}\star\begin{pmatrix}\bm{1}_{25}\bm{1}^{T}_{n}\\ \\ \bm{1}_{5}\bm{1}^{T}_{n}\\ \\ \bm{1}_{20}\bm{1}^{T}_{n}\\ \end{pmatrix}.$

      $\theta_{G(g_{1}=1)}$ and $\theta_{G(g_{1}=2)}$ are random observations obtained from Ber$(1-\pi_{1})$ to simulate the effect of the groups at level $L=1$. In order to reflect the overlapping significance of $G(g_{1}=1)$ and $G(g_{1}=2)$ we draw a sample $\theta_{G(g_{1}=1)\cap G(g_{1}=2)}\sim$ Ber$(1-\pi^{*}_{L=1})$, where $\pi^{*}_{1}\leq\pi_{1}$. Each of these effects are imposed on the respective hypotheses by $\bm{\Theta}_{L=1}$.

   3. (c)

      $\bm{\theta}_{L=2}$ as a random vector of $m=50$ i.i.d. Ber$(1-\pi_{2})$, reflects the grouping effect at level $L=2$.

2. 2.

   Obtain

   $\displaystyle\bm{\Theta}=\bm{\Theta}_{L=0}\star\bm{\Theta}_{L=1}\star\bm{\theta}_{L=2}\bm{1}^{T}_{n},$

   (4.2)

   where $A\star B$ denotes the Hadamard product between matrices $A$ and $B$, and $\bm{1}_{a}$ representing the $a$-dimensional vector of $1$’s.

3. 3.

   Given $\bm{\Theta}$, generate a random $m\times n$ matrix $\mathbf{X}=((X_{gh}))$ as follows:

   	$\displaystyle\mathbf{X}=$	$\displaystyle\quad\mu\bm{\Theta}+\sqrt{(1-\rho_{L1})(1-\rho_{L2})}\mathbf{Z}_{mn}+\sqrt{(1-\rho_{L1})\rho_{L2}}\mathbf{Z}_{m}\bm{1}_{n}^{T}+$
   		$\displaystyle\quad\sqrt{\rho_{L1}(1-\rho_{L2})}\bm{1}_{m}\mathbf{Z}_{n}^{T}+\sqrt{\rho_{L1}\rho_{L2}}Z_{0}\bm{1}_{m}\bm{1}_{n}^{T},$		(4.3)

   having generated $\mathbf{Z}_{mn}$ as $m\times n$ random matrix, $\mathbf{Z}_{m}$ as $m$-dimensional random vector, and $\mathbf{Z}_{n}$ as $n$-dimensional random vector, where $m=50,n=100$, each comprising i.i.d. $N(0,1)$ samples, and $Z_{0}$ as an additional $N(0,1)$ random variable. The quantities $\rho_{L1}$ and $\rho_{L2}$ measure the strength of correlation among the p-values in each group at levels $L=1$ and $L=2$ respectively. For simplicity, we assume same $\rho_{L1}$ and $\rho_{L2}$ in all groups at respective levels of grouping.

4. 4.

   Apply each of the DAHeir-GBH and BH procedures at FDR level $\alpha=0.05$, and the p-filter process with $\alpha=0.05$ at each level of grouping, for testing $H_{i}:\text{E}(X_{i})=0$ against $K_{i}:\text{E}(X_{i})>0$, simultaneously for all $i=1,\ldots,N$, in terms of the corresponding weighted p-values. Note the proportions of false rejections among all rejections and correct rejections among all false nulls.

5. 5.

   Repeat Steps 1-4 500 times to simulate the values of FDR and power for each procedure by averaging out the corresponding proportions obtained in Step 4.

The formulation in equation (4.2) allows us to determine the significant state of each hypothesis, subject to the state of significance of its parent groups at all levels. This helps to demonstrate the true effect of the underlying hierarchical structure in the simulation studies. Since the groups at level $L=1$ overlap, we can regulate the density of signals in intersecting region using the following

and in the non-intersecting regions by the following

Here $\pi$ represents the proportion of true nulls in the entire set of $mn$ hypotheses in terms of the proportions of groups at level 2 i.e., ($1-\pi_{2}$), and at level $1$, i.e, ($1-\pi_{1}$) (or $(1-\pi^{*}_{1})$ in the overlapping region) containing signals, and the proportion of individual hypotheses that are significant $(1-\pi_{0})$, when no classification is imposed, i.e., at level $L=0$.

We set $\mu=0$ for true null hypotheses and $=3$ for true signals. We also set $\pi_{2}=0.5$, $\pi_{L=1}=0.5$, $\pi^{*}_{L=1}=0.25$ and increase $1-\pi_{0}$ from $0$ through $1$ in order to simulate different densities of signals using (4.4) and (4.5).

When the p-values satisfy the PRDS condition, the DAHeir-GBH still controls FDR under certain conditions. We set $\pi_{1}=\pi^{*}_{1}=\pi_{2}=0$ so that the signals are uniformly distributed in all groups at levels one and two, and the variability comes from $\pi_{0}=0.3$. We choose $\lambda<\alpha=0.05$ for the DAHeir-GBH and the adaptive BH procedures. For the p-filter algorithm, we set $\lambda^{(m)}=\lambda$, for $m=0,1,2$. The within group dependence at levels one and two are characterized by $\rho_{L1}=0.3$ and $\rho_{L2}=0.4$ respectively. The comparative performances of the three procedures are shown in figure 6. All procedures conservatively control the FDR at all levels of $\lambda$, and the proposed DAHeir-GBH has comparable power with the adaptive BH.

Interested readers are referred to Nandi et al. (2021) for numerical results on the simultaneous classification setup.

The simulation results here and in Nandi et al. (2021) assert that our proposed methods in the hierarchically as well as simultaneously grouped setting are able to control FDR, at least under independence, and are more powerful than existing procedures. Since the Gen-GBH and DAGen-GBH combine the hierarchical and simultaneous procedures, their performances would evidently be better than existing competing methods.

We apply our DAGen-GBH procedure to analyze the data in http://kdd.ics.uci.edu/databases/eeg/eeg.data.html, which arises from a large study examining EEG correlates of genetic predisposition to alcoholism. The experiment was conducted on a control group and a treatment group comprising of alcoholic subjects. For our analysis, we consider $10$ subjects each in the control group and alcoholic group. On each subject, ten trials were performed. During each trial, a picture from Snodgrass and Vanderwart (1980) was presented to the subject, while EEG activity in the form of voltage fluctuations (in microvolts) were recorded at $256$ time points (each time point being $1/256$th of a second) from $61$ electrodes placed on the subject’s scalp. Figure 7 shows an example of the EEG pattern, averaged over measurements obtained from ten trials, for two subjects, one each from the control and alcoholic groups. The $x$ and $y$ axes correspond to time and scalp electrodes, respectively, and the $z$-axis corresponds to the EEG measurements. It illustrates a clear difference in the voltage fluctuation patterns for the two subjects over time and electrodes, which leads us to the following driving question underlying this study:

At which pairs of electrodes, are the EEG measurements significantly different between the two groups of subjects?

Answering this question can shed further insight into brain dysfunction and regions impacted by alcoholism. In the literature, case-control studies on EEG measurements, such as this, has often been used to diagnose brain dysfunction and evaluate the performance of relevant treatment procedures (for examples, see Centorrino et al. (2002), Oscar-Berman and Marinkovic (2007), etc.). Specifically, we consider simultaneous testing of the following null hypotheses.