A Unified and Optimal Multiple Testing Framework based on $\rho$-values

Some of the most popular FDR rules use p-value as significance index for ranking (e.g., Benjamini and Hochberg, , 1995; Cai et al., , 2022; Lei and Fithian, , 2018; Genovese et al., , 2006). The standard way of obtaining p-values is to apply a probability integral transform to some well-known test statistics. For example, Li and Barber, (2019) uses a permutation test, Roquain and Van De Wiel, (2009) employs a Mann–Whitney U test and Cai et al., (2022) adopts a t-test. However, the p-value based methods can be inefficient because the conventional p-values do not incorporate information from the alternative distributions (e.g., Sun and Cai, , 2007; Leung and Sun, , 2021). The celebrated Neyman–Pearson lemma states that the optimal statistic for testing a single hypothesis is the likelihood ratio, while an analog of the likelihood ratio statistic for multiple testing problems is the Lfdr (Efron et al., , 2001; Efron, , 2003; Aubert et al., , 2004). It has been shown that a ranking and thresholding procedure based on Lfdr is asymptotically optimal for FDR control (Sun and Cai, , 2007; Cai et al., , 2019). Nevertheless, the validity of many Lfdr based methods relies crucially on the consistent estimation of the Lfdr statistics (Basu et al., , 2018; Gang et al., , 2023; Xie et al., , 2011), which can be difficult in practice (e.g., Cai et al., , 2022). To overcome this, the weighted p-value based approaches are proposed to emulate the Lfdr method; see, for example, Lei and Fithian, (2018), Li and Barber, (2019) and Liang et al., (2023). However, all of those Lfdr-mimicking procedures either require strong model assumptions or are suboptimal.

To address the aforementioned issues arising from p-value and Lfdr frameworks, this article proposes a new concept called $\rho$-value, which takes the form of likelihood ratios and allows wide flexibility in the choice of density functions. Based on such $\rho$-values (or weighted $\rho$-values, similarly as weighted p-values), we aim to develop a new flexible multiple testing framework that unifies p-value based and Lfdr based approaches. Specifically, the proposed $\rho$-value based methods also operate in two steps: first rank all hypotheses according to (weighted) $\rho$-values (with the goal of mimicking Lfdr ranking), then reject those hypotheses with (weighted) $\rho$-values less than or equal to some threshold, where the threshold is determined similarly as in p-value based procedures.

Compared to existing frameworks, the new framework has several advantages. First, if the $\rho$-values are carefully constructed, then their ranking coincides with that produced by Lfdr statistics. Thus, methods based on $\rho$-values can be asymptotically optimal. Second, the strategies for choosing the threshold for $\rho$-value based methods are similar to those for p-value based methods. Therefore, the proposed approaches enjoy theoretical properties similar to methods based on p-values. Moreover, the FDR guarantee of the $\rho$-value approaches does not require consistent estimations of Lfdr statistics, and hence the proposed approaches are much more flexible than the Lfdr framework. Third, compared to Lfdr methods, side information can be easily incorporated into $\rho$-value based procedures through the simple weighting scheme to improve the ranking of the hypotheses. Finally, we provide a unified view of p-value based and Lfdr based frameworks. In particular, we show that these two frameworks are not as different as suggested by previous works (e.g., Sun and Cai, , 2007; Leung and Sun, , 2021).

The paper is structured as follows. Section 2 presents the problem formulation. In Section 3, we introduce $\rho$-values and provide several examples of $\rho$-value based procedures. Sections 4 and 5 present numerical comparisons of the proposed methods and other competing approaches using simulated and real data, respectively. More discussions of the proposed framework are provided in Section 6. Technical proofs are collected in the Appendix.

As introduced by the previous sections, $\rho$-value is an analog of the likelihood ratio, and we now present its rigorous definition below.

It can be seen from the above definition that, the $\rho$-values are broadly defined. In particular, $\rho$-value can be viewed as a generalization of the likelihood ratio. The key difference is that the denominator of $\rho$-value can be any density function and is no longer required to be $f_{1}(X)$.

In addition, $\rho$-value has a close connection with e-value, which is defined below.

Using Markov inequality, it is straightforward to show that the reciprocal of an e-value is a p-value (Wang and Ramdas, , 2022). Note that $\mathbb{E}(1/\rho)=\int\frac{g(x)}{f_{0}(x)}f_{0}(x)dx=\int g(x)dx=1$, which means that $1/\rho$ is an e-value and a $\rho$-value is a special p-value. Therefore, one can directly use $\rho$-values as the inputs for the BH procedure and obtain an e-BH procedure (Wang and Ramdas, , 2022). One attractive feature of the e-BH procedure is that it controls FDR at the target level under arbitrary dependence. However, in practice the e-BH procedure is usually very conservative. To see why, let $c(\cdot)$ be the distribution function of $\rho_{i}\equiv f_{0}(X_{i})/g(X_{i})$ under $H_{0,i}$ and assume that $\rho_{i}$’s are independent. Then if we reject the hypotheses with $\rho_{i}\leq t$, the e-BH procedure uses $mt$ as a conservative estimate of the number of false positives despite the fact that $\rho_{i}$ no longer follows $\mbox{Unif}(0,1)$ under the null. In the following section, we introduce a novel $\rho$-BH procedure which improves e-BH by using a tighter estimate of the number of false positives.

For simplicity, we assume that the $\rho$-values $\rho_{i}\equiv f_{0}(X_{i})/g(X_{i})$ are independent under the null and there are no ties among them. Recall $c(\cdot)$ is the null distribution function of $\rho_{i}$’s. If the null density $f_{0}$ is known, then $c(\cdot)$ can be easily estimated by the empirical distribution of the samples generated from $f_{0}$; the details will be provided in the simulation section. Consider the decision rule that rejects $H_{0,i}$ if and only if $\rho_{i}\leq t$, $i=1,\ldots,m$, then a tight estimate of the number of false positives is $m(1-\pi)c(t)$, where $\pi$ is defined in Model (1). Hence, a reasonable and ideal threshold $t$ should control the estimated false discovery proportion (FDP) at level $\alpha$, i.e., $\frac{m(1-\pi)c(t)}{\#\{i:\rho_{i}\leq t\}\vee 1}\leq\alpha$. To maximize power, we can choose the rejection threshold to be $\rho_{(k)}$, where $k=\text{max}_{1\leq j\leq m}\left\{c(\rho_{(j)})\leq\frac{\alpha j}{m(1-\pi)}\right\}$ with $\rho_{(j)}$ the $j$-th order statistic. We call this procedure the $\rho$-BH procedure and summarize it in Algorithm 1.

It is important to note that, since $c(\cdot)$ is a monotonically increasing function, the rankings produced by $\rho_{i}$ and $c(\rho_{i})$ are identical. Also note that $c(\rho_{i})\sim\mbox{Unif}(0,1)$ under the null. It follows that the $\rho$-BH procedure is in fact equivalent to the original BH procedure with $c(\rho_{i})$ as the p-value inputs and $\alpha/(1-\pi)$ as the target FDR level. Consequently, the $\rho$-BH procedure enjoys the same FDR guarantee as the original BH procedure, as presented in Theorem 1, and has power advantage over e-BH in the meantime.

It is worthwhile to note that, one can use any predetermined $g(\cdot)$ to construct $\rho$-values and Theorem 1 always guarantees FDR control if $X_{i}$’s are independent under the null. Recall the Lfdr statistic is defined as $\mbox{Lfdr}_{i}\equiv\mathbb{P}(\theta_{i}=0|X_{i})=\frac{(1-\pi)f_{0}(X_{i})}{(1-\pi)f_{0}(X_{i})+\pi f_{1}(X_{i})}$, and it is shown in Sun and Cai, (2007); Cai et al., (2019) that a ranking and thresholding rule based on Lfdr is asymptotically optimal among all mFDR control rules. Thus, ideally we would like to adopt $\{f_{0}(X_{i})/f_{1}(X_{i}),i=1,\ldots,m\}$ as the $\rho$-values since the ranking produced by $f_{0}(X_{i})/f_{1}(X_{i})$’s is identical to that produced by Lfdr statistics. In fact, with this choice of $\rho$-values, the $\rho$-BH procedure is asymptotically optimal under some mild conditions as stated by the next theorem.

In practice, $f_{1}$ and $\pi$ are usually unknown and need to be estimated from the data. The problem of estimating non-null proportion has been discussed extensively in the literatures (e.g., Storey, , 2002; Jin and Cai, , 2007; Meinshausen and Rice, , 2006; Chen, , 2019). To ensure valid mFDR control, we require the estimator $\hat{\pi}$ to be conservative consistent, defined as follows.

One possible choice of such $\hat{\pi}$ is the Storey estimator as provided by the following proposition.

The problem of estimating $f_{1}$ is more complicated. If we use the entire sample $\{X_{i}\}_{i=1}^{m}$ to construct $\hat{f}_{1}$ and let $\rho_{i}=f_{0}(X_{i})/\hat{f}_{1}(X_{i})$, then $\rho_{i}$’s are no longer independent even if $X_{i}$’s are. One possible strategy to circumvent this dependence problem is to use sample splitting. More specifically, we can randomly split the data into two disjoint halves and use the first half of the data to estimate the alternative density for the second half, i.e., $\hat{f}_{1}^{\scriptscriptstyle(2)}$ (e.g., we can use the estimator proposed in Fu et al., (2022)), then the $\rho$-values for the second half can be calculated by $f_{0}(X_{i})/\hat{f}_{1}^{\scriptscriptstyle(2)}(X_{i})$. Hence, when testing the second half of the data, $\hat{f}_{1}^{\scriptscriptstyle(2)}$ can be regarded as predetermined and independent of the data being tested. The decisions on the first half of the data can be obtained by switching the roles of the first and the second halves and repeating the above steps. If the FDR is controlled at level $\alpha$ for each half, then the overall mFDR is also controlled at level $\alpha$ asymptotically. We summarize the above discussions in Algorithm 2 and Theorem 3.

We have mentioned in Section 3.1 that a $\rho$-value is the reciprocal of an e-value. Therefore, if $\{\rho_{i}\}_{i=1}^{m}$ are dependent with possibly complicated dependence structures, one can simply employ $\{\rho_{i}\}_{i=1}^{m}$ as the inputs of the BH procedure and the resulting FDR can still be appropriately controlled at the target level.

Alternatively, one can deal with arbitrary dependence by the proposal in Benjamini and Yekutieli, (2001). They have shown that the BH procedure with target level $\alpha$ controls FDR at level $\alpha S(m)$, where $S(m)=\sum_{i=1}^{m}\frac{1}{i}\approx\log(m)$ is known as the Benjamini -Yekutieli (BY) correction. Since the $\rho$-BH procedure is equivalent to the BH procedure with $c(\rho_{i})$ as p-values and target FDR level $\alpha/(1-\pi)$, we also have the following theorem on BY correction.

Note that if one chooses to use BY correction or the e-BH procedure, then the data splitting step in the data-driven procedures is no longer necessary. We also remark that, we can extend our FDR control results under independence in the previous sections to weakly dependent $\rho$-values and same technical tools as in Cai et al., (2022) can be employed to achieve asymptotic error rates control. The current article focuses on the introduction of the new testing framework based on $\rho$-values and skips the discussions of the weakly dependent scenarios.

Similar to incorporating prior information via a p-value weighting scheme (e.g., Benjamini and Hochberg, , 1997; Genovese et al., , 2006; Dobriban et al., , 2015), we can also employ such weighting strategy in the current $\rho$-value framework. Let $\{w_{i}\}_{i=1}^{m}$ be a set of positive weights such that $\sum_{i=1}^{m}w_{i}=m$. The weighted BH procedure proposed in Genovese et al., (2006) uses $p_{i}/w_{i}$’s as the inputs of the original BH procedure. Genovese et al., (2006) proves that, if $p_{i}$’s are independent and $\{w_{i}\}_{i=1}^{m}$ are independent of $\{p_{i}\}_{i=1}^{m}$ conditional on $\{\theta_{i}\}_{i=1}^{m}$, then the weighted BH procedure controls FDR at level less than or equal to $\alpha$.

Following their strategy, we can apply the weighted BH procedure to $c(\rho_{i})$’s and obtain the same FDR control result. However, such procedure might be suboptimal as will be explained in the following section. Alternatively, we derive a weighted $\rho$-BH procedure and the details are presented in Algorithm 3.

Note that $\{\rho_{i}/w_{i}\}$ in Algorithm 3 produces a different ranking than $\{c(\rho_{i})/w_{i}\}$, which may improve the power of the weighted p-value procedure with proper choices of $\rho$-values and weights. On the other hand, the non-linearity of $c(\cdot)$ imposes challenges on the theoretical guarantee for the mFDR control of Algorithm 3 compared to that of Genovese et al., (2006), and we derive the following result based on similar assumptions as in Theorem 3.

It is worthwhile to note that, Genovese et al., (2006) requires $\sum_{i=1}^{m}w_{i}=m$, which makes the weighted p-value procedure conservative. In comparison, Algorithm 3 no longer requires such condition, and it employs a tight estimate of the FDP that leads to a more powerful testing procedure.

When the oracle parameters are unknown, we can similarly construct a data-driven weighted $\rho$-BH procedure with an additional data splitting step as in Algorithm 2. Due to the space limit, the details are presented in Algorithm 6 in Section A of the Appendix.

In this section, we propose a $\rho$-BH procedure that incorporates the auxiliary information while maintaining the proper mFDR control. In many scientific applications, additional covariate informations such as the patterns of the signals and nulls are available. Hence, the problem of multiple testing with side information has received much attention and is becoming an active area of research recently (e.g., Lei and Fithian, , 2018; Li and Barber, , 2019; Ignatiadis and Huber, , 2021; Leung and Sun, , 2021; Cao et al., , 2022; Zhang and Chen, , 2022; Liang et al., , 2023). As studied in the aforementioned works, proper use of such side information can enhance the power and the interpretability of the simultaneous inference methods. Let $X_{i}$ denote the primary statistic and $s_{i}\in\mathbb{R}^{l}$ the side information. Then we model the data generation process as follows.

Upon observing $\{(X_{i},s_{i})\}_{i=1}^{m}$, we would like to construct a decision rule $\boldsymbol{\delta}$ that controls mFDR based on $\rho$-values. As before, we assume the null distributions $f_{0}(\cdot|s_{i})$ are known, and we define the $\rho$-values by $\rho_{i}=\frac{f_{0}(X_{i}|s_{i})}{g(X_{i}|s_{i})}$ for some density function $g(\cdot|s_{i})$. Let $\eta:\mathbb{R}^{l}\rightarrow(0,1)$ be a predetermined function and $c_{i}(\cdot)$ be the null distribution of $\rho_{i}$. Then we incorporate the side information through a $\rho$-value weighting scheme by choosing an appropriate function $\eta(\cdot)$. The details are summarized in Algorithm 4.

In contrast to the $\rho$-BH procedure, Algorithm 4 is no longer equivalent to the BH procedure with $\{c_{i}(\rho_{i})/w_{i}\}$’s as the inputs since the rankings produced by $\{c_{i}(\rho_{i})/w_{i}\}$’s and $\{\rho_{i}/w_{i}\}$’s are different. The ideal choice of $g(\cdot|s_{i})$ is again $f_{1}(\cdot|s_{i})$, while the ideal choice of $\eta(\cdot)$ is $\pi(\cdot)$, the rationale is provided as follows. Define the conditional local false discovery rate (Clfdr, Fu et al., (2022); Cai et al., (2019)) as

Cai et al., (2019) shows that a ranking and thresholding procedure based on Clfdr is asymptotically optimal for FDR control. Note that if we take $g(\cdot|s_{i})$ to be $f_{1}(\cdot|s_{i})$ and $\eta(\cdot)$ to be $\pi(\cdot)$, then the ranking produced by $\rho_{i}/w_{i}$’s is identical to that produced by Clfdr statistics. However, the validity of the data-driven methods proposed in Cai et al., (2019) and Fu et al., (2022) relies on the consistent estimation of Clfdr_i’s. In many real applications, it is extremely difficult to accurately estimate Clfdr even when the dimension of $s_{i}$ is moderate (Cai et al., , 2022). In contrast, the mFDR guarantee of Algorithm 4 does not rely on any of such Clfdr consistency results and our proposal is valid under much weaker conditions as demonstrated by the next theorem.

Note that, the validity of the above theorem allows flexible choices of functions $g(\cdot|s_{i})$ and the weights $w_{i}$. Hence, similarly as the comparison between $\rho$-value and Lfdr, the $\rho$-value framework with side information is again much more flexible than the Clfdr framework that requires the consistent estimation of the Clfdr statistics.

We also remark that, Cai et al., (2022) recommends using $\frac{\pi(s_{i})}{1-\pi(s_{i})}$ to weigh p-values arising from the two sample t-statistics, and the authors only provide a heuristic explanation on the superiority of using $\frac{\pi(s_{i})}{1-\pi(s_{i})}$ over $\frac{1}{1-\pi(s_{i})}$ as weights, while we have the following optimality result which provides a more rigorous justification for $\frac{\pi(s_{i})}{1-\pi(s_{i})}$.

In practice, we need to choose $\eta(\cdot)$ and $g(\cdot|s_{i})$ based on the available data $\{(X_{i},s_{i})\}_{i=1}^{m}$. Again, if the entire sample is used to construct $\eta(\cdot)$ and $g(\cdot|s_{i})$, then the dependence among $w_{i}$’s and $\rho_{i}$’s is complicated. Similar to Algorithm 2, we can use sample splitting to circumvent this problem. The details of the data-driven version of Algorithm 4 is provided in Algorithm 5. To ensure a valid mFDR control, we require a uniformly conservative consistent estimator of $\pi(\cdot)$, whose definition is given below.

The problem of constructing such uniformly conservative consistent estimator $\hat{\pi}(\cdot)$ has been discussed in the literatures; see for example, Cai et al., (2022) and Cai et al., (2019).

The next theorem shows that Algorithm 5 indeed controls mFDR at the target level asymptotically under conditions analogous to those assumed in Theorem 3.