Asymptotic Uncertainty of False Discovery Proportion

Multiple hypothesis testing has been a popular topic in statistical research; it has wide applications in many scientific areas, such as biology, medicine, genetics, neuroscience, and finance. Consider a multiple testing problem where we need to test $p$ hypotheses simultaneously; denote by $p_{0}$ and $p_{1}$ the number of true null and false null hypotheses respectively. The testing outcomes can be summarized in Table 1,

where $V$ and $R$ denote the number of hypotheses that are wrongly rejected (i.e., the number of type I errors made) and rejected (i.e., the number of rejections made) out of the $p$ hypotheses, respectively.

Conventional methods propose to control the familywise error rate (FWER): $\text{FWER}=P(V>0)$ (Bonferroni, 1936; Šidák, 1967; Holm, 1979; Simes, 1986; Holland and Copenhaver, 1987; Hochberg, 1988; Rom, 1990), or the generalized familywise error rate (gFWER): $\text{gFWER}=P(V>k),\text{ for }k\geq 1$ (Dudoit et al., 2004; Pollard and van der Laan, 2004; Lehmann and Romano, 2012). These methods control the probability of zero or a limited number of false positives, therefore they are usually conservative, especially when the number of hypotheses $p$ is large. We observe that large $p$ is common in practice. For example, in a genome-wide association study, hundreds of thousands or even millions of single nucleotide polymorphisms (SNPs) of each subject and the associated disease status are measured. It is of interest to find which SNPs are associated with the disease status. A common approach is to conduct a hypothesis test for each SNP; this results in a huge number of hypothesis tests.

When $p$ is large, FDP and its statistical characteristics have been introduced. False discovery rate (FDR) has been the most popularly studied; it was introduced by Benjamini and Hochberg (1995) and was proven to be a reasonable criterion; in particular,

$$\mathrm{FDP}=\frac{V}{R},\qquad\mathrm{FDR}=E(\mathrm{FDP}),$$

where $\mathrm{FDP}=\mathrm{FDR}=0$ if $R=0$. A list of FDR research can be found in Benjamini and Hochberg (1995), Benjamini et al. (2001), Storey (2002), Storey and Tibshirani (2003), Storey et al. (2004), Ferreira and Zwinderman (2006), Clarke et al. (2009), and the references therein. We observe that most of these methods are founded on the assumption that the test statistics are independent or satisfy some restrictive dependence structure.

Because of the complicated nature of the practical data, the resultant test statistics may follow any arbitrary dependence structure. The effects of dependence on FDR have been studied extensively, e.g., by Benjamini et al. (2001), Finner and Roters (2002), Owen (2005), Sarkar (2006), Efron (2007), among others. Recently, new methods for the control and/or estimation of FDR have been proposed to allow complex and general dependence structure or even incorporate it in the testing procedure. Examples include the local index of significance (LIS) testing procedure for dependence modeled by a hidden Markov model (Sun and Cai, 2009), the pointwise and clusterwise analysis for spatial dependence (Sun et al., 2015), and the principal factor approximation (PFA) method for arbitrary dependence (Fan et al., 2012; Fan and Han, 2017; Fan et al., 2019). In particular, Fan et al. (2012) studied FDP and FDR under arbitrary dependence with a two step procedure. They first proposed a PFA method to reduce arbitrary dependence to weak dependence; and then established theoretically the strong consistency of the FDP and FDR. Therefore, their method is valid for test statistics under arbitrary dependence; under weak dependence, both FDP and FDR converge to the same asymptotic limit.

In addition to FDR, another popular multiple testing criterion is called false discovery exceedance (FDX, also referred to as TPPFP in the literature), which is defined as the tail probability of the event that FDP exceeds a threshold $q\in(0,1)$:

$$\text{FDX}=P(\mathrm{FDP}>q).$$

Compared to FDR that is only a single-value summary, FDX has the potential to fully describe the uncertainty and/or distribution of FDP as the threshold value $q$ varies. However, it has not been fully characterized how the dependence among the test statistics would impact FDX as opposed to FDR in the literature. In particular, some works on FDX assume independence among tests (Genovese and Wasserman, 2004; Guo and Romano, 2007; Ge and Li, 2012); more importantly, most works focus on developing procedures that control FDX as a single-value characteristic of FDP other than investigating how dependence affects the uncertainty and/or distribution of FDP (Korn et al., 2004; van der Laan et al., 2004; Lehmann and Romano, 2005; Genovese and Wasserman, 2006; Delattre and Roquain, 2015; Hemerik et al., 2019; Döhler and Roquain, 2020; Basu et al., 2021). Thus, these works still lack the valuable information about how FDP varies from study to study.

As far as we are aware, most of the above works have focused on proposing the estimates of FDP and/or its related statistical characteristics. On the one hand, studying these estimates is sufficient to ensure that the corresponding FDP and/or its statistical characteristics are under control. On the other hand, we observe that the asymptotic uncertainty of the FDP may rely heavily on the dependence structure of the corresponding test statistics; and it is of great practical value to quantify this variability, as it can serve as an indicator of the quality of the FDP estimate and benefit the research of the FDR and FDX methods. More specifically, with the FDP variance estimate available, one can conclude more confidently how well the FDP is controlled when controlling its corresponding mean at a given level. It may also benefit the development on the methods for estimating and controlling FDX. The research on the uncertainty of FDP has high potential impact, but it is still limited in the community. For example, Ge and Li (2012) and Delattre and Roquain (2011) derived the asymptotic distribution of FDP under independence and Gaussian equi-correlation dependence, respectively; Delattre and Roquain (2016) established the asymptotic distribution of $\mathrm{FDP}$ under some assumptions on the dependence structure of the test statistics and that the expected values of the test statistics in the alternatives are equivalent.

In this paper, our focus is to provide a rigorous theoretical investigation of the asymptotic uncertainty of FDP in the framework that the test statistics are weakly dependent. In particular, we derive the theoretical expansion of the FDP: it is a linear combination of the tests, based on which we obtain its asymptotic variance. We also examine how this variance varies under different dependence structures among the test statistics both theoretically and numerically. One interesting observation is that the FDP based on weakly dependent test statistics converges to the same limit as that based on independent test statistics (Fan et al., 2012); however, the asymptotic variance of the former is almost always significantly greater than the latter both theoretically and numerically. With these observations, we recommend that in a multiple testing performed by an FDP procedure, we may report both the mean and the variance estimates of FDP to enrich the study outcome. Finally, with a real GWAS dataset, we demonstrate how our findings may impact real studies.

Consider a genome-wide association study with $n$ subjects; for each subject, its disease status and $p$ SNPs are measured. The data can be represented by $\mathbf{Y}=(Y_{1},\ldots,Y_{n})^{T}$ and a $n\times p$ matrix $\mathbf{X}=(X_{ij})_{i=1,\ldots,n;j=1,\ldots,p}$, where $Y_{i}$ denotes the disease status of subject $i$, and the $i$th row of $\mathbf{X}$ collects its associated SNPs. When $Y_{i}$’s are measured in a continuous scale, the association between the $j$th SNP and the disease status can be modelled by the marginal linear regression model:

$$Y_{i}=\alpha_{j}+\beta_{j}X_{ij}+\epsilon_{ij},\quad i=1,\ldots,n,$$

(1)

where $\alpha_{j}$ and $\beta_{j}$ are regression parameters and $\epsilon_{ij}$ are random errors. The estimates $\hat{\beta}_{j},\ j=1,\ldots,p$ and the corresponding test statistics $Z_{j},\ j=1,\ldots,p$ for

$$H_{0j}:\beta_{j}=0\quad\text{versus}\quad H_{1j}:\beta_{j}\neq 0,\quad j=1,\dots,p,$$

(2)

can be derived with a standard least squares procedure. Proposition 1 in Fan et al. (2012) verified that $\hat{\beta}_{j},\ j=1,\ldots,p$ are correlated and derived their joint distribution, under appropriate regularity conditions. As a consequence, the test statistics $Z_{j},\ j=1,\ldots,p$ are also correlated. In particular, they derived that conditioning on $\mathbf{X}$,

$$(Z_{1},\ldots,Z_{p})^{T}\sim N((\mu_{1},\ldots,\mu_{p})^{T},\boldsymbol{\Sigma}),$$

(3)

where $(\mu_{1},\ldots,\mu_{p})^{T}$ and $\boldsymbol{\Sigma}=(\sigma_{ij})_{i,j=1,\ldots,p}$ are appropriate population mean and variance matrix for $(Z_{1},\ldots,Z_{p})^{T}$. Furthermore, the testing problem (2) becomes equivalent to

$$H_{0j}:\mu_{j}=0\quad\text{versus}\quad H_{1j}:\mu_{j}\neq 0,\quad j=1,\ldots,p.$$

(4)

Therefore, we shall work on the multiple testing problem under the setup (3) and (4) hereafter. Let $\mathcal{H}_{0}=\{j:H_{0j}\text{ is true}\}$ and $\mathcal{H}_{1}=\{j:H_{0j}\text{ is false}\}$ respectively be the sets of indices for the true nulls and the true alternatives, and let $p_{0}$ and $p_{1}$ denote their cardinalities. For a given threshold $t\in(0,1)$, the $j$th test is given by $t_{j}=\mathrm{1}(|Z_{j}|>|z_{t/2}|)$, where $t_{j}=1$ indicates that the null hypothesis is rejected, and $t_{j}=0$ otherwise; here $z_{t/2}$ denotes the $(t/2)$th quantile of a standard normal distribution. Define $V(t)=\sum_{j\in\mathcal{H}_{0}}t_{j}$, $S(t)=\sum_{j\in\mathcal{H}_{1}}t_{j}$, and $R(t)=\sum_{j=1}^{p}t_{j}$; they are respectively the number of false discoveries, correct discoveries, and the total number of discoveries. As a consequence, we have $\mathrm{FDP}(t)=V(t)/R(t)$.

Fan et al. (2012) studied the estimation of $\mathrm{FDP}(t)$ under an arbitrary dependence structure of $\boldsymbol{\Sigma}$. One milestone finding is that if $\boldsymbol{\Sigma}$ follows the weak dependence structure:

$$\sum_{i,j}|\sigma_{ij}|=O\left(p^{2-\delta}\right)\text{ for some }\delta>0,$$

(5)

then,

$$\lim_{p\to\infty}\left[\mathrm{FDP}(t)-\frac{p_{0}t}{\sum_{i=1}^{p}\left\{\Phi(z_{t/2}+\mu_{i})+\Phi(z_{t/2}-\mu_{i})\right\}}\right]=0,$$

(6)

almost surely, where $\Phi(\cdot)$ denotes the cumulative distribution function of a standard normal distribution. This indicates that under the aforementioned weak dependence assumption, $\mathrm{FDP}(t)$ converges to a limit not depending on the covariance structure of the test statistics. It is noteworthy that this limit holds not only for $\mathrm{FDP}(t)$, but also for $\mathrm{FDR}(t)$.

We observe that this conclusion does not hold for the asymptotic variance of $\mathrm{FDP}(t)$; the asymptotic variance of $\mathrm{FDP}(t)$ varies under different covariance structures among the test statistics even when they are weakly dependent. In this article, we shall derive the asymptotic variance of $\mathrm{FDP}(t)$ theoretically and study how it varies under different dependence structures of the test statistics numerically.

In this section, we study the asymptotic variance of $\mathrm{FDP}(t)$ in the framework of Section 2. We first establish the asymptotic expansion of $\mathrm{FDP}(t)$, and its asymptotic variance under a weak dependence structure. Then, we show that with mild conditions, this asymptotic variance under dependence is higher than that under independence. This implies that under weak dependence, although the asymptotic mean of $\mathrm{FDP}(t)$ remains the same, its asymptotic variance can be inflated. For space limitations and presentational continuity, we have relegated the technical details in the supplementary materials.

In Section 3.1, we have established the asymptotic variance of FDP in the framework of weak dependence. We observe that this variance depends on the parameters $(\mu_{1},\ldots,\mu_{p})^{T}$ and $\boldsymbol{\Sigma}$. We need to estimate them so that we can quantify this variance to make it practically applicable. Throughout, we assume that $\boldsymbol{\Sigma}$ is known or can be estimated from other resources (see Section 2) and estimate only $(\mu_{1},\ldots,\mu_{p})^{T}$.

One direct approach is to estimate $(\mu_{1},\ldots,\mu_{p})^{T}$ by their corresponding test statistics $(Z_{1},\ldots,Z_{p})^{T}$. However, in practical data, we often have $p_{0}\gg p_{1}$; this implies that most $\mu_{j}$’s are zeros. Estimating these “zeros” using the corresponding observed $Z_{j}$’s tends to diminish the efficiency of the method, as they can be viewed as noise. In our numerical studies, we consider the following estimation approach.

• (1)

  Find an estimate $\hat{\pi}_{0}$ for $\pi_{0}=p_{0}/p$, which is the proportion of the number of tests that belong to $\mathcal{H}_{0}$ over the total number of tests.

  Estimation methods for $\pi_{0}$ are available in the literature. For example, Storey and Tibshirani (2003) proposed $\hat{\pi}_{0}(\lambda)=\sum_{i=1}^{p}1(P_{i}>\lambda)/\{m(1-\lambda)\}$ with an $\lambda\in(0,1)$; the basic idea of this estimator is to pretend that $p$-values greater than $\lambda$ are from $\mathcal{H}_{0}$. In this framework, smoothed and bootstrap estimators of $\pi_{0}$ are proposed by Storey and Tibshirani (2003) and Storey et al. (2004), respectively. Additionally, Langaas et al. (2005) proposed a nonparametric maximum likelihood estimator of $\pi_{0}$. More recently, Wang et al. (2010) proposed a “SLIM” estimator of $\pi_{0}$; their estimator accounted for the dependence among the test statistics. We shall adopt these estimators in our numerical studies and compare their performance in the estimation of the FDP variance.

• (2)

  Estimate $p_{1}$ by $\hat{p}_{1}=p(1-\hat{\pi}_{0})$. We then collect the hypotheses with the top $\hat{p}_{1}$ statistic values to form $\widehat{\mathcal{H}}_{1}$, and the rest are assigned to $\widehat{\mathcal{H}}_{0}$.

• (3)

  Set $\hat{\mu}_{j}=0$ for $j\in\widehat{\mathcal{H}}_{0}$; and set $\hat{\mu}_{j}=Z_{j}$ for $j\in\widehat{\mathcal{H}}_{1}$; and use them to estimate the variance of $\mathrm{FDP}(t)$ based on Corollary 1.

We conduct simulation studies to validate our developments in Sections 3 and 4. In Section 5.1, we use various numerical examples to illustrate how the asymptotic variances of FDP are different for varied dependence structures among the test statistics. In Section 5.2, we illustrate how the choice of the methods for estimating $\hat{\pi}_{0}$ reviewed in Section 4 affects the estimation for the variance of FDP.

We apply our method to a genome-wide association study (GWAS); more specifically, an eQTL mapping study. The data are formed of 210 subjects from the International Hapmap Project, including three groups: (1) 90 subjects from Asian formed of 45 Japanese from Tokyo, 45 Han Chinese from Beijing, named “Grp1”; (2) 60 Utah residents with ancestry from northern and western Europe, named “Grp2”; and (3) 60 Nigerian from Ibadan, named “Grp3”. More details of the data can be found from the web link: ftp://ftp.sanger.ac.uk/pub/genevar/; see also Bradic et al. (2011). We adopt the existing methods to construct test statistics for testing the association between the expression level of the gene CCT8 and the SNP genotypes for each group of the subjects (Bradic et al., 2011; Fan et al., 2012).

The SNP genotypes have three possibilities: $\text{SNP}=0$ means no polymorphism, $\text{SNP}=1$ indicates one nucleotide has polymorphism, and $\text{SNP}=2$ implies both nucleotides have polymorphisms. Adopting the strategy in Bradic et al. (2011), we transform the SNP data into two dummy variables, namely $(d_{1},d_{2})$:

• •

  $(d_{1},d_{2})=(0,0)$ for $\text{SNP}=0$;

• •

  $(d_{1},d_{2})=(1,0)$ for $\text{SNP}=1$;

• •

  $(d_{1},d_{2})=(0,1)$ for $\text{SNP}=2$.

For subject $i$, denote by $Y_{i}$ and $(d_{1,i,j},d_{2,i,j})$ the logarithm-transformed value of the expression level of gene CCT8 and the aforementioned dummy variables of the $j$th SNP. We consider two sets of marginal linear regression models:

	$\displaystyle\mbox{Set 1:}\qquad Y_{i}$	$\displaystyle=$	$\displaystyle\alpha_{1,j}+\beta_{1,j}d_{1,i,j}+\epsilon_{1,i,j},\quad\mbox{for }j=1,\ldots,p;$
	$\displaystyle\mbox{Set 2:}\qquad Y_{i}$	$\displaystyle=$	$\displaystyle\alpha_{2,j}+\beta_{2,j}d_{2,i,j}+\epsilon_{2,i,j},\quad\mbox{for }j=1,\ldots,p.$

Here we impute the missing SNP data by 0, and exclude the SNPs that are identical to each other for all the subjects.

For each group of the data, we combine aformentioned two sets of models, and follow the procedure in Fan et al. (2012) to generate the test statistics; this leads to $2p$ test statistics. In detail, we use the procedure described in Section 2 to establish an initial set of test statistics and their corresponding covariance matrix. We then use the dependence-adjusted procedure in Fan et al. (2012) to remove the top ten principal factors from the covariance matrix, and take a standardization step to generate a new set of test statistics of unit marginal variances; this helps reduce the dependence among the test statistics.

With the resultant test statistics for each group, we estimate the asymptotic limits and standard deviations of $\mathrm{FDP}(t)$ based on the following setups:

• •

  For each group of test statistics, we consider two thresholds of $p$-value, which are chosen such that the estimates of the $\mathrm{FDP}(t)$ limits are reasonable values.

• •

  The asymptotic limits are estimated from (6), denoted by $\widehat{\text{FDP}}(t)$.

• •

  The asymptotic standard deviations are computed based on the methods in Section 4. As we have observed in Section 5.2 that the standard deviation estimates based on the $\pi_{0}$ estimates from the “Langaas” and “SLIM” methods outperform the others; therefore, we present only the results from these methods. Furthermore, to assess how dependence may affect the estimation of these standard deviation estimates, we report these estimates based on the assumptions that the test statistics are independent and dependent, respectively. Therefore, for each group of subjects, we have four estimates, respectively denoted by $S_{\widehat{\text{FDP}}(t)}^{\text{Langaas}}$, $S_{\widehat{\text{FDP}}(t),\text{Ind}}^{\text{Langaas}}$, $S_{\widehat{\text{FDP}}(t)}^{\text{SLIM}}$, and $S_{\widehat{\text{FDP}}(t),\text{Ind}}^{\text{SLIM}}$.

We summarize our estimation results in Table 5. From this table, we observe that under different assumptions of whether the dependence is present among the test statistics significantly affect the asymptotic standard deviation of $\mathrm{FDP}(t)$. The ratios $S_{\widehat{\text{FDP}}(t)}^{\text{Langaas}}/S_{\widehat{\text{FDP}}(t),\text{Ind}}^{\text{Langaas}}$ and $S_{\widehat{\text{FDP}}(t)}^{\text{SLIM}}/S_{\widehat{\text{FDP}}(t),\text{Ind}}^{\text{SLIM}}$ range from 4 to 7 folds. Moreover, the ratios of the FDP standard deviation estimates to its limit values $\widehat{\text{FDP}}(t)$ are between 0.25 and 0.82 and thus the uncertainty of FDP is nonnegligible. These observations reinforce the statement that in an FDP related inference, it may not be sufficient to consider only the limit value of FDP (or equivalently, FDR); we suggest that it is important to correctly account for the dependence; otherwise, it might be over optimistic in the conclusion of how well the FDP is controlled. For example, when $t=0.005$ in Grp1, the 95 percentile (say, two standard deviations above the mean) of the $\mathrm{FDP}(t)$ is estimated as $0.315+0.153\times 2=0.621$ under dependence; in contrast, if we assume that the test statistics are independent, the corresponding estimate is $0.315+0.031\times 2=0.377$. Therefore, ignoring the dependence leads to an over optimistic belief of how well the upper bound of FDP is controlled; such a nonignorable difference may significantly change the statistical interpretation of the results in practice.

In this article, we have established the theoretical results for the asymptotic variance of $\mathrm{FDP}$ in the framework that the test statistics are weakly dependent. We observe that with mild conditions, the asymptotic variance of $\mathrm{FDP}$ under dependence is usually greater than that under independence; this is also demonstrated by our numerical studies: even under weak dependence, the variances of the FDPs can be numerically much larger than those under independence. It is thus crucial for us to take into account the dependence in the estimation of the variance of the FDP and report it together with its mean. Furthermore, we have observed in the real data example that the standard deviation estimate of $\mathrm{FDP}$ under the dependence assumption is significantly different from (usually greater than) that under the independence assumption. In practice, we suggest that we presumably include the dependence in the estimation, unless there is a strong evidence of otherwise. This will be helpful for us to better understand how well the FDP is controlled, and avoid over optimistic conclusions.

In this article, we have studied the asymptotic uncertainty of FDP in the same framework as Fan et al. (2012). Under the assumption that the test statistics are weakly dependent, we obtain the asymptotic expansion of $\mathrm{FDP}(t)$ and establish the explicit formula for computing the asymptotic variance of $\mathrm{FDP}(t)$. Within this framework, we observe many interesting open problems that can be studied in the future. We list several as examples; there are many others.

First, our proposed method for estimating the variance of $\mathrm{FDP}(t)$ in Section 4 relies on the estimation of $\pi_{0}$; there, we adopted existing methods in the literature. From simulation studies in Section 5.2, we observe that the performance of the variance estimate for $\mathrm{FDP}(t)$ relies heavily on that of $\pi_{0}$; although the Langaas and SLIM methods perform reasonably well, we believe that they are not yet optimal, as none of these $\pi_{0}$ estimation methods have sufficiently accounted for the dependence among the test statistics. Therefore, we expect that establishing improved $\pi_{0}$ estimation methods that can better accommodate the dependence among the test statistics would be of particular interest, and can be incorporated into our estimation method in Section 4.

Second, we have imposed some weak dependence assumptions in our development; this may limit the method in some applications. Practically, one possible solution is not to use the original test statistics, but the dependence-adjusted ones, as they are more powerful and often weakly dependent (Friguet and Causeur, 2011; Fan et al., 2012); yet there is no theoretical guarantees. Therefore, it is challenging, but would be of interest to study the asymptotic variance of $\mathrm{FDP}(t)$ under arbitrary dependence.

Third, we assume that the test statistics are normally distributed in this article. It would be of great both practical and theoretical value to generalize our study to accommodate multiple testing problems where the test statistics follow other popular distributions, e.g., $t$, $F$, and chi-square distributions. We refer to Zhuo et al. (2020) as an example of the recent works on dependent $t$-tests and leave the relevant investigation to our future work.

SUPPLEMENTARY MATERIALS

The supplementary materials contain technical details for the theorems in Section 3.