False Discovery Rate Control For Structured Multiple Testing: Asymmetric Rules And Conformal Q-values

Abstract

The effective utilization of structural information in data while ensuring statistical validity poses a significant challenge in false discovery rate (FDR) analyses. Conformal inference provides rigorous theory for grounding complex machine learning methods without relying on strong assumptions or highly idealized models. However, existing conformal methods have limitations in handling structured multiple testing, as their validity often requires the deployment of symmetric decision rules, which assume the exchangeability of data points and permutation-invariance of fitting algorithms. To overcome these limitations, we introduce the pseudo local index of significance (PLIS) procedure, which is capable of accommodating asymmetric rules and requires only pairwise exchangeability between the null conformity scores. We demonstrate that PLIS offers finite-sample guarantees in FDR control and the ability to assign higher weights to relevant data points. Numerical results confirm the effectiveness and robustness of PLIS and demonstrate improvements in power compared to existing model-free methods in various scenarios. Supplementary materials for this article are available online, including a standardized description of the materials available for reproducing the work.

Keywords:

• Conformal inference
• Generalized e-values
• Local index of significance
• Pairwise exchangeability

1 Introduction

In various applications, data is commonly acquired and organized as ordered sequences or lattices, revealing informative structural patterns that play a crucial role in analysis and interpretation. For example, spatio-temporal data in econometric analyses may display serial or spatial dependence structures, while in genome-wide association studies, single-nucleotide polymorphisms (SNPs) often cluster along biological pathways, indicating functional relationships between genes. To make reliable and meaningful inferences, it is essential to account for the underlying structural patterns in such data. False discovery rate (FDR) methods (Benjamini and Hochberg 1995) are a powerful tool for identifying sparse signals from massive and complex data. In FDR analyses, a significant challenge lies in integrating structural information to enhance statistical power and interpretability, while simultaneously ensuring the validity of FDR control. In this section, we discuss recent advancements in addressing this challenge, as well as the limitations of existing approaches, and present our contributions.

1.1 Model-Based and Model-Free FDR Methods

Structural knowledge, such as the clustering patterns and dependence, can be leveraged to enhance the efficiency of existing FDR methods, as demonstrated by the works of Benjamini and Heller (2007), Sun and Cai (2009), Fan, Han, and Gu (2012), Sun et al. (2015), Rebafka, Roquain, and Villers (2022), and Perrot-Dockès et al. (2023). However, the validity of existing “model-based” FDR methods, which involve constructing new test statistics based on estimated model parameters, typically depend on idealized assumptions, such as correct specification of the data generating models, homogeneity and stationarity of the underlying processes, consistent estimation of unknown parameters, and asymptotic normality of the test statistics. These assumptions may not be fulfilled and may be hard to check in practice. The violation of these assumptions can have serious consequences in FDR analysis, including decreased power and inflated error rate. It remains a significant challenge to develop powerful FDR methods for structured multiple testing that are assumption-lean and provably valid.

The framework of conformal inference (Vovk, Gammerman, and Shafer 2005) provides finite-sample uncertainty guarantees on a flexible class of off-the-shelf machine learning algorithms, only under the assumption of exchangeability. The important connection between machine learning and FDR analysis, as highlighted by Yang et al. (2021), Mary and Roquain (2022), and Bates et al. (2023), indicates that efficient conformity scores, and hence provably valid and powerful conformal p-values, can be constructed based on complex learning algorithms. This insightful perspective leads to “model-free” approaches to FDR control, a research direction with considerable promise. Recent developments along this direction include the BONuS (Yang et al. 2021) and AdaDetect (Marandon et al. 2024) procedures, which convincingly demonstrate that combining conformal p-values with the Benjamini-Hochberg (BH) algorithm enables the use of highly effective machine learning models while ensuring finite-sample FDR control, even in the presence of model misspecification. However, currently available conformal methods are not well-equipped to handle structured multiple testing. The limitation is illustrated next.

1.2 Non-Exchangeable Scores and Asymmetric Decision Rules

Suppose we are interested in testing m null hypotheses $H_{0,i}$, for $1\le i\le m$, based on summary statistics $\mathbf{X}=(X_i:1\le i\le m)$. Denote $\mathbf{\Theta }=({\theta }_i:1\le i\le m)\in {\{0,1\}}^m$ the true states of nature, where ${\theta }_i=0/1$ indicates that $H_{0,i}$ is true/false. Let ${\mathcal{H}}_0=\{i\in [m]:H_{0,i}\text{ is true}\}$. The decisions are represented by a binary vector $\mathit{\delta }(\mathbf{X})=({\delta }_i:1\le i\le m)\in {\{0,1\}}^m$, where ${\delta }_i=1$ indicates that $H_{0,i}$ is rejected and ${\delta }_i=0$ otherwise. We allow ${\delta }_i$ to depend on the entire dataset $\mathbf{X}$. We call $\mathit{\delta }$a symmetric decision rule (Copas 1974) if $\mathit{\delta }\{\Pi (\mathbf{X})\}=\Pi \{\mathit{\delta }(\mathbf{X})\}$ for all permutation operators $\Pi$.

To provide context, consider a hidden Markov model (HMM) where the underlying states $\mathbf{\Theta }$ are unknown, and the observations are conditionally independent given $\mathbf{\Theta }$. In this toy example, we assume that the hidden states form a binary Markov chain and the non-null cases tend to appear in clusters. We further assume that the null and non-null distributions are respectively given by $(X_i|{\theta }_i=0)\sim \mathcal{N}(0,1)$ and $(X_i|{\theta }_i=1)\sim \mathcal{N}(2,1).$ Suppose we have observed $X_j=X_k=2.5$, $j\ne k$, where $X_j$ is surrounded by small observations and $X_k$ is surrounded by large observations. Intuitively, $X_k$ is more likely to be a non-null case compared to $X_j$ due to the clustering effects from the Markovian structure; hence, symmetric rules are inappropriate. Sun and Cai (2009) demonstrated that the optimal FDR procedure in HMMs is an asymmetric rule that relies on thresholding the local index of significance (LIS), whose value depends on the entire sequence, with higher weights given to neighboring locations.

The presence of structured patterns within the data and the use of asymmetric rules pose a significant challenge to the conformal inference framework. To see this, we tentatively consider the following definition, which generalizes the conformal p-value in Bates et al. (2023) in a nuanced (and possibly improper) way: (1) $$\begin{array}{c}\begin{array}{cc}{p}_i\equiv p_i(X_i)\hfill & =\frac{1}{1+\left|{\mathcal{D}}^{\mathit{\text{cal}}}\right|}\left[1+\sum _{j\in {\mathcal{D}}^{\mathit{\text{cal}}}}\mathbb{I}\{s_i(X_i)>s_j(Y_j)\}\right],\hfill \\ \hfill & \hspace{1em}i\in \left[m\right]=\{1,\dots ,m\}.\hfill \end{array}\hfill \end{array}$$(1)

Here $s_i(·)$ is the conformity score function that is allowed to vary across i, ${\mathcal{D}}^{\mathit{\text{cal}}}$ is the index set for calibration data containing null samples $Y_j$, with $|·|$ denoting the cardinality of a set. Bates et al. (2023) showed that if the score functions are identical and permutation-invariant, and the data points $\{X_i,i\in {\mathcal{H}}_0;Y_j,j\in D^{\mathit{\text{cal}}}\}$ are jointly exchangeable, then the conformal p-values in (1) are super-uniform and PRDS.¹ Hence, according to Benjamini and Yekutieli (2001) and Sarkar (2002), applying BH with these conformal p-values is valid for FDR control. However, score functions that leverage structural patterns are typically not permutation-invariant, leading to a violation of the joint exchangeability assumption. Consequently, the conformal p-values defined by (1) may become improper and compromise the PRDS property, making it problematic to apply the BH procedure.

1.3 A Preview of Our Method and Contributions

To develop a provably valid FDR procedure that can effectively leverage structural information and accommodate asymmetric rules, we propose the pseudo local index of significance (PLIS) procedure that consists of three steps:

1. constructing baseline data that preserve useful structural patterns;

2. calculating conformity scores based on user-specified working models;

3. constructing a mirror process that emulates the true process for decision-making.

The proposed algorithm features two critical components: an innovative algorithm for calculating the conformity score through a working model, which captures useful structural patterns and enhances the efficiency of the FDR analysis; a generic framework for constructing a mirror process for decision-making, which bypasses the p-value inference framework, eliminates the need for joint exchangeability assumptions, and establishes finite-sample FDR theory for asymmetric rules.

The proposed research makes several contributions. First, PLIS provides a provably valid model-free approach to multiple testing. Compared to model-based FDR procedures, it achieves comparable power when the underlying model is correctly specified and accurately estimated, while maintaining finite-sample FDR control in scenarios where the model is mis-specified or estimated poorly. Second, PLIS demonstrates superior power compared to existing conformal methods due to its ability to leverage asymmetric rules that exploit structural information. Finally, we have developed a novel theoretical framework that relies only on the pairwise exchangeability between $s_i(X_i)$ and $s_i(Y_i)$ under the null, which is much weaker than the joint exchangeability requirement stipulated in existing theories for conformal inference.

1.4 Organization

The article is organized as follows. Section 2 introduces the PLIS procedure for structured multiple testing and establishes its theoretical properties. Section 3 explores extensions of PLIS and its relationship to existing concepts. Section 4 presents simulation results to assess the numerical performance of PLIS and to compare it with existing methods. An illustration of the proposed method is provided in Section 5 through an application. The supplementary material contains proofs, extensions and additional numerical results.

2 Structured Multiple Testing: Conformal Inference with Asymmetric Rules

In this section, we first present the structured probabilistic model (Section 2.1) and the corresponding problem formulation (Section 2.2). Next, we introduce the PLIS procedure (Section 2.3) and establish its theoretical properties (Section 2.4). Concrete examples and guidelines are provided in Section 2.5 to illustrate the PLIS framework. Sections 2.6 and 2.7, respectively discuss the semi-supervised PLIS algorithm and the conformal q-value notion.

2.1 A Class of Structured Probabilistic Models

Consider a multiple testing problem with binary-valued unknown states $\mathbf{\Theta }=({\theta }_i:i\in \mathcal{G})$ that form a graph $\mathcal{G}$. Let $m=\left|\mathcal{G}\right|$ be the number of nodes/hypotheses, where ${\theta }_i=0$ indicates that node i is a null case, and ${\theta }_i=1$ otherwise. Our study focuses on a class of structured probabilistic models (e.g., Goodfellow, Bengio, and Courville 2016), where the correlations between random variables $\mathbf{X}=(X_i:i\in \mathcal{G})$ are captured by the interdependence structure between corresponding latent states $\mathbf{\Theta }=({\theta }_i:i\in \mathcal{G})$. The inference units may be conceptualized as the nodes of the graph, and the interdependence structures between test statistics are encoded as edges that connect the nodes. The observations $(X_i:i\in \mathcal{G})$ are conditionally independent given $\mathbf{\Theta }$, obeying (2) $$\mathbb{P}\left(\mathbf{X},\mathbf{\Theta }\right)=\mathbb{P}\left(\mathbf{\Theta }\right)\prod _{i\in \mathcal{G}}f\left(X_i|{\theta }_i\right),\hspace{1em}{X}_i|{\theta }_i\sim (1-{\theta }_i)f_0(x)+{\theta }_i{f}_{1i}(x),$$(2) where $f_0$ and $f_{1i}$ are the null and non-null densities, respectively. Although the model assumes conditional independence, it still remains highly flexible, as we do not impose any distributional constraints on the unknown states $\mathbf{\Theta }$, and allow $f_{1i}$ to vary across different inference units. Conditional independence is closely related to the widely employed exchangeability assumption in the conformal inference literature (Schervish 2012; Barber et al. 2023). As established by de Finetti’s Theorem (Heath and Sudderth 1976; Diaconis and Freedman 1980; Durrett 2019), if a set of random elements is exchangeable, there must exist a latent $\eta$ such that the random elements are iid conditional on $\eta$. For a more detailed and in-depth exploration of the generalization of model (2) and the exchangeability assumption, please refer to Section 2.6 and Section C of the supplementary material.

We make two assumptions without loss of generality: (a) $f_0(-x)=f_0(x)$; and (b) a larger magnitude of $\left|X_i\right|$ provides stronger evidence against the null hypothesis. In situations where the assumptions do not hold, we may transform $X_i$ into z-values using the formula $z_i={\Phi }^{-1}\{F_0(X_i)\}$, where $F_0$ represents the cumulative distribution function (CDF) of $X_i$ under the null, and $\Phi$ represents the CDF of a standard Gaussian variable.

We discuss several special cases of Model (2). First, we examine the scenario where the latent states ${\theta }_i$ are iid and the non-null densities $f_{1i}$ are identical (and hence denoted as $f_1$). Under this assumption, Model (2) reduces to Efron’s two-group model (Efron et al. 2001): (3) $$X_i\stackrel{\text{iid}}{\sim }(1-\pi )f_0(x)+\pi f_1(x),$$(3) where $\pi =\mathbb{P}({\theta }_i=1)$ is the proportion of non-null cases. While the independent case is not the primary focus of this work, it is worth noting that our proposed methodology remains applicable in this scenario. In Section 2.5, we present a tailored version of our proposal for Model (3), and discuss its connection to the AdaDetect procedure (Marandon et al. 2024).

When the latent states $\mathbf{\Theta }$ form a homogeneous binary Markov chain, a two-state hidden Markov model (HMM, Rabiner 1989) can be recovered from (2). HMMs have been extensively studied in the context of multiple testing (e.g., Sun and Cai 2009; Sesia, Sabatti, and Candès 2018; Perrot-Dockès et al. 2023). By properly accounting for the dependence between test statistics and leveraging the HMM structure, these methods can effectively control the FDR and provide more accurate inference. Section 2.5 presents a customized algorithm for HMMs based on our proposal and illustrates its superiority over existing methods.

The class of models in (2) also includes useful models such as the Ising model (Onsager 1944) and conditional random field model (Lafferty, McCallum, and Pereira 2001), which have been applied in multiple testing for climate change analysis and network analysis (Liang and Nettleton 2010; Shu, Nan, and Koeppe 2015; Liu, Zhang, and Page 2016; Rebafka, Roquain, and Villers 2022). In contrast to HMMs which assume homogeneous transition probabilities and identical emission probabilities, model (2) offers greater flexibility to accommodate deviations from these assumptions. Moreover, in contrast to the stationary hidden states model proposed by Wu (2008), Model (2) allows the underlying process for $\mathbf{\Theta }$ to be nonstationary. As many real-world phenomena exhibit nonstationary behaviors, Model (2) is better equipped to provide a more accurate representation of the underlying dynamics.

2.2 Problem Formulation and Basic Setup

We consider a multiple testing problem that arises from Model (2), which can be complex in nature and difficult to learn accurately from data. The hypotheses of interest are $H_{0,i}:{\theta }_i=0\text{ versus }{H}_{1,i}:{\theta }_i=1,$ for $i\in \mathcal{G}$. A key feature of the models in (2) is the exchangeability of null observations $(X_i:i\in \mathcal{G},{\theta }_i=0)$, which remains true regardless of the model complexity. This important characteristic provides a strong foundation for conformal inference and has been widely assumed. If we make the additional assumption that our working model has learned a score function that is permutation-invariant, then we can properly define the conformal p-values in (1) for FDR analysis.

As previously indicated, we want to employ a score tailored specifically to index i. The score does not require permutation invariance regarding i. This attribute provides a substantial degree of flexibility, enabling the construction of powerful scores that can effectively leverage useful structural patterns. However, the violation of permutation invariance poses a significant challenge to the conformal inference framework, as the p-value is no longer super-uniform under the null. To address this issue, Section 2.3 develops a new framework that eliminates the need of defining p-values as required by existing methods.

In conformal inference, the calibration dataset, denoted as $\mathbf{Y}$, plays a key role in assessing the relative significance of hypotheses. Various approaches can be employed to obtain the calibration data $\mathbf{Y}$. In the classical setting of multiple testing, $\mathbf{Y}$ can be sampled from the known null distribution $F_0$. Conversely, in certain machine learning applications, such as the outlier detection problem, $\mathbf{Y}$ may be sampled from a set of labeled null observations. For the present, we assume that the null distribution $F_0$ is known. In Section 2.6, we turn to scenarios where $F_0$ is unknown, yet labeled null samples are available.

Let $s_i^X=s_i(X_i)$ denote the test score of unit i. The superscript X is used to distinguish $s_i^X$ from the calibration score, denoted as $s_i^Y$, which are computed based on the calibration data. Let ${\mathbf{S}}_X=\{s_i^X:i\in \mathcal{G}\}$ and ${\mathbf{S}}_Y=\{s_i^Y:i\in \mathcal{G}\}$. In constructing the scores, it is crucial to guarantee the pairwise exchangeability between $s_i^X$ and $s_i^Y$ for $i\in {\mathcal{H}}_0$. This important notion is rigorously defined in Section 2.4.

We denote the decision for test unit i as ${\delta }_i\in \{0,1\}$, where ${\delta }_i=1$ signifies the rejection of $H_{0,i}$ and ${\delta }_i=0$ otherwise. Denote $\mathit{\delta }=({\delta }_i:1\le i\le m)$. The false discovery proportion (FDP) and true discovery proportion (TDP) are respectively defined as (4) $$\text{FDP}(\mathit{\delta })=\frac{\sum _{j\in \mathcal{G}}(1-{\theta }_j){\delta }_j}{(\sum _{j\in \mathcal{G}}{\delta }_j)\vee 1},\hspace{1em}\text{TDP}(\mathit{\delta })=\frac{\sum _{j\in \mathcal{G}}{\theta }_j{\delta }_j}{(\sum _{j\in \mathcal{G}}{\theta }_j)\vee 1}.$$(4)

The (frequentist) FDR and average power (AP) can be defined accordingly: $\text{FDR}=\mathbb{E}[\text{FDP}(\mathit{\delta })]\text{ and AP}=\mathbb{E}[\text{TDP}(\mathit{\delta })],$ where the expectation operator is taken over the observed and calibration data while holding the hidden states $\Theta$ fixed. The objective is to develop a decision rule that effectively controls the FDR while maximizing the AP.

2.3 The PLIS Procedure for FDR Control

The conformity score that best captures structural patterns from probabilistic models (2) is $\mathbb{P}({\theta }_i=0|\mathbf{X})$, assuming perfect knowledge of the model structure. However, accurately estimating this score is challenging. To address this, we propose to compute a score through a user-specified working model that efficiently extracts structural information and predicts unknown states. Specifically, we choose a working model $\mathcal{M}$ and corresponding algorithm $\mathcal{A}$ to compute pseudo scores ${\mathbf{S}}_X=\{s_i^X:i\in \mathcal{G}\}$, bearing in mind that the working model may deviate from the complex data-generating model.

To effectively leverage the structural information, the size of the calibration set needs to match that of the test set, ensuring that each test unit is paired with a corresponding calibration data point. Suppose the test and calibration data points have been paired and denoted as $\{(X_i,Y_i),i\in \mathcal{G}\}$. In our methodology, the initial step involves constructing the baseline dataset, denoted as $\mathbf{W}=(W_i\equiv h(X_i,Y_i):i\in \mathcal{G})$, where (5) $$\begin{array}{c}h(x,y)=\{\begin{array}{cc}x& \text{if }\left|x\right|\ge \left|y\right|\\ y& \text{otherwise}\end{array}.\hfill \end{array}$$(5)

The construction of $\mathbf{W}$ ensures that $W_i$ depends equally on $X_i$ and $Y_i$, while satisfying the flipping-coin property, that is, $\mathbb{P}(W_i=X_i)=\mathbb{P}(W_i=Y_i)=1/2$ if $i\in {\mathcal{H}}_0$.

Remark 1

Any symmetric function satisfying $h(x,y)=h(y,x)$ can be employed to generate the baseline data $\mathbf{W}$. However, our choice to use (5) enables us to retain data with large effect sizes, thereby preserving the structural patterns among the non-null cases. This is particularly advantageous for datasets that display local dependency structures. Section E.7 of the supplement shows that (5) leads to substantial power gain compared to alternative choices such as $h^{\prime }(x,y)=x+y$.

We proceed to generate two parallel datasets $\{{\tilde{\mathbf{X}}}_i:i\in \mathcal{G}\}$ and $\{{\tilde{\mathbf{Y}}}_i:i\in \mathcal{G}\}$ by respectively substituting $X_i$ and $Y_i$ for $W_i$ in $\mathbf{W}$. For example, if $\mathbf{X}$ is an ordered sequence $(X_1,\dots ,X_m)$, then ${\tilde{\mathbf{X}}}_i=(W_1,\dots ,W_{i-1},X_i,W_{i+1},\dots ,W_m)$. Next, we calculate the pseudo score $s_i^X\equiv {\mathbb{P}}_{\mathcal{M},\mathcal{A}}({\theta }_i=0|{\tilde{\mathbf{X}}}_i)$, referred to as the pseudo local index of significance (PLIS). The corresponding calibration score is given by $s_i^Y\equiv {\mathbb{P}}_{\mathcal{M},\mathcal{A}}({\theta }_i=0|{\tilde{\mathbf{Y}}}_i)$, applying the same working model $\mathcal{M}$ and computational algorithm $\mathcal{A}$ to the calibrated data ${\tilde{\mathbf{Y}}}_i$. Figure 1 visually illustrates the process of data generation and score computation, which guarantees the pairwise exchangeability between $s_i^X$ and $s_i^Y$, a key property that will be rigorously defined and established in Section 2.4.

Figure 1: A graphical illustration of constructing the baseline data and calculating the scores. To create ${\tilde{\mathbf{X}}}_i$ and ${\tilde{\mathbf{Y}}}_i$, each $W_i$ in $\mathbf{W}$ is replaced by $X_i$ and $Y_i$ in turn, for $i\in \mathcal{G}$, with the values in remaining nodes unchanged. The same function $s_k(·)$ is used for computing $s_k^X$ and $s_k^Y$, ensuring that the two scores are pairwise exchangeable if $k\in {\mathcal{H}}_0$.

Suppose that a smaller value of $s_i^X$ indicates stronger evidence against the null hypothesis. Intuitively, the rejection of $H_{0,i}$ requires two pieces of evidence: first, $s_i^X$ must be sufficiently small, and second, $s_i^X$ should be smaller than its null counterpart $s_i^Y$. Therefore, we only select small $s_i^X$ from the candidate rejection set ${\mathcal{G}}^r=\{i\in \mathcal{G}:s_i^X<s_i^Y\}$. Consider a class of thresholding rules $\mathit{\delta }=({\delta }_i:i\in \mathcal{G})$, where ${\delta }_i\equiv 0$ for $i\notin {\mathcal{G}}^r$, and ${\delta }_i=\mathbb{I}(s_i^X\le t)$ for $i\in {\mathcal{G}}^r$. Correspondingly, we introduce the calibration set ${\mathcal{G}}^c=\{i\in \mathcal{G}:s_i^Y<s_i^X\}$, which mirrors ${\mathcal{G}}^r$. This mirror set of null scores $\{s_i^Y:i\in {\mathcal{G}}^c\}$ allows us to characterize the joint behavior of their counterparts $\{s_i^X:i\in {\mathcal{G}}^r\}$ under the null. It follows that the true FDP process and its mirror process, denoted $Q(t)$, are respectively given by (6) $$\begin{array}{cc}\text{FDP}(t)\hfill & =\frac{\sum _{j\in {\mathcal{G}}^r\cap {\mathcal{H}}_0}\mathbb{I}\left\{s_j^X\le t\right\}}{\left(\sum _{j\in {\mathcal{G}}^r}\mathbb{I}\left\{s_j^X\le t\right\}\right)\vee 1},\hfill \\ Q(t)\hfill & =\frac{1+\sum _{j\in {\mathcal{G}}^c}\mathbb{I}\left\{s_j^Y\le t\right\}}{\left(\sum _{j\in {\mathcal{G}}^r}\mathbb{I}\left\{s_j^X\le t\right\}\right)\vee 1},\hspace{1em}t>0.\hfill \end{array}$$(6)

To gain a better understanding of our methodology, it is helpful to note that $Q(t)$ provides a conservative estimate of $\text{FDP}(t)$. This conservativeness arises from two factors: (a) $\sum _{j\in {\mathcal{G}}^c}\mathbb{I}\{s_j^Y\le t\}$ is always greater than or equal to $\sum _{j\in {\mathcal{G}}^c\cap {\mathcal{H}}_0}\mathbb{I}\{s_j^Y\le t\}$; and (b) $\sum _{j\in {\mathcal{G}}^c\cap {\mathcal{H}}_0}\mathbb{I}\{s_j^Y\le t\}$ can be seen as a mirror process that resembles the unobserved process $\sum _{j\in {\mathcal{G}}^r\cap {\mathcal{H}}_0}\mathbb{I}\{s_j^X\le t\}$ since, as will be established in the next section (Property 3), $s_j^X$ and $s_j^Y$ are pairwise exchangeable for $j\in {\mathcal{H}}_0$. We choose the largest threshold, denoted as $\tau$, such that the conservative estimate of the FDP falls below $\alpha$: $$\tau =\mathrm{\text{sup}}\left\{t\in {\mathbf{S}}_X\cup {\mathbf{S}}_Y:Q(t)\le \alpha \right\}.$$

The proposed PLIS procedure, summarized in Algorithm 1, rejects $H_{0,i}$ if ${\delta }_i\equiv \mathbb{I}(i\in {\mathcal{G}}^r,s_i^X\le \tau )=1$. By mathematical conventions, the supremum of an empty set is defined as $-\infty$. Therefore, if the set $\{t\in {\mathbf{S}}_X\cup {\mathbf{S}}_Y:Q(t)\le \alpha \}$ is empty, then $\tau$ is set to $-\infty$, and Algorithm 1 does not reject any hypotheses.

Algorithm 1

The PLIS procedure

Input: The observations $\mathbf{X}=(X_i:i\in \mathcal{G})$, the null distribution $F_0$, a pre-specified FDR level $\alpha$, a working model $\mathcal{M}$ and a computational algorithm $\mathcal{A}$.

Output: A decision rule $\mathit{\delta }=({\delta }_i:i\in \mathcal{G})$.

1: Generate $\mathbf{Y}=(Y_i:i\in \mathcal{G})$, where $Y_i$ are iid observations from $F_0$.

2: Create $\mathbf{W}=(W_i:i\in \mathcal{G})$, where $W_i=h(X_i,Y_i)$ and $h(x,y)$ is given by (5). Construct ${\tilde{\mathbf{X}}}_i$ and ${\tilde{\mathbf{Y}}}_i$ by respectively substituting $X_i$ and $Y_i$ in place of $W_i$ in $\mathbf{W}$ for $i\in \mathcal{G}$.

3: Calculate $s_i^X\equiv {\mathbb{P}}_{\mathcal{M},\mathcal{A}}({\theta }_i=0|{\tilde{\mathbf{X}}}_i)$ and $s_i^Y\equiv {\mathbb{P}}_{\mathcal{M},\mathcal{A}}({\theta }_i=0|{\tilde{\mathbf{Y}}}_i)$ for $i\in \mathcal{G}$ with the user-specified working model $\mathcal{M}$ and algorithm $\mathcal{A}$. Let ${\mathcal{G}}^r=\{i:s_i^X<s_i^Y\}$.

4: Let $\tau =\mathrm{\text{sup}}\{t\in {\mathbf{S}}_X\cup {\mathbf{S}}_Y:Q(t)\le \alpha \}$, where $Q(t)$ is defined in (2.3).

5: Let ${\delta }_i=\mathbb{I}\{s_i^X\le \tau \}$ for $i\in {\mathcal{G}}^r$, otherwise ${\delta }_i=0$.

6: Return $\mathit{\delta }=({\delta }_i:i\in \mathcal{G})$.

We conclude this section with several remarks. First, in order to ensure pairwise exchangeability between $s_i^X$ and $s_i^Y$ under the null, it is crucial to estimate the working model $\mathcal{M}$ from $\mathbf{W}$ before calculating the conformity scores. Detailed explanations and practical guidelines for the deployment of PLIS under two widely-used working models can be found in Section 2.5. Second, we can define an alternative estimate of $\text{FDP}(t)$, denoted as $Q^{\prime }(t)$, by substituting $\mathcal{G}$ in place of both ${\mathcal{G}}^r$ and ${\mathcal{G}}^c$ in (2.3). However, employing $Q^{\prime }(t)$ tends to be overly conservative; further elaborations on this are provided in Sections 3.3. Third, Algorithm 1 closely resembles the operation of the conformal BH, albeit with several modifications. Finally, as revealed by an insightful reviewer, $Q(t)$ in (2.3) is closely related to the Selective SeqStep + algorithm (Barber and Candès 2015). Further discussions on these connections can be found in Section D of the supplement.

2.4 Theory

We start by establishing several exchangeability properties, first for the data points and subsequently for the conformity scores. The random elements $\{Z_i:i\in [K]\}$ are said to be (jointly) exchangeable if the distribution of $(Z_1,\dots ,Z_K)$ is the same as that of $(Z_{{\Pi }_1},\dots ,Z_{{\Pi }_K})$ for any permutation $({\Pi }_1,\dots ,{\Pi }_K)$ of the indices $\{1,\dots ,K\}$. We write $(Z_1,\dots ,Z_K)\stackrel{d}{=}(Z_{{\Pi }_1},\dots ,Z_{{\Pi }_K}),$ where $\stackrel{d}{=}$ denotes equality in distribution. We will soon introduce a weaker notion, the pairwise exchangeability, which is more relevant to our theory. The following exchangeability properties are proven in Section A of the supplement.

The first property can be easily deduced from the definition of Model (2) and the construction of $\mathbf{W}$ outlined in Algorithm 1.

Property 1 (Exchangeability and conditional independence between data points)

Suppose $\mathbf{X}=(X_i:i\in \mathcal{G})$ are observations from Model (2), and $(Y_i:i\in \mathcal{G})$ are randomly drawn from the null distribution $F_0$. Then we have: (a) The random variables $\left\{Y_i,i\in \mathcal{G};X_j,j\in {\mathcal{H}}_0\right\}$ are jointly exchangeable conditional on $\Theta$; (b) $\{Z_j=\{X_j,Y_j\}:j\in \mathcal{G}\}$ are conditionally independent given $\mathbf{W}$ and $\mathbf{\Theta }$; (c) For $i\in {\mathcal{H}}_0$, $X_i$ and $Y_i$ are pairwise exchangeable conditional on $\mathbf{W}$ and $\mathbf{\Theta }$, that is, $(X_i,Y_i|\mathbf{W},\mathbf{\Theta })\stackrel{d}{=}(Y_i,X_i|\mathbf{W},\mathbf{\Theta }).$

The next property characterizes the pairwise exchangeability between ${\tilde{\mathbf{X}}}_j$ and ${\tilde{\mathbf{Y}}}_j$.

Property 2 (Pairwise exchangeability and conditional independence between datasets)

Consider ${\tilde{\mathbf{X}}}_j$ and ${\tilde{\mathbf{Y}}}_j$ as generated according to Algorithm 1. Then we have: (a) $\{{\mathcal{D}}_j=\{{\tilde{\mathbf{X}}}_j,{\tilde{\mathbf{Y}}}_j\}:j\in \mathcal{G}\}$ are conditionally independent given $\mathbf{W}$ and $\mathbf{\Theta }$; (b) For $i\in {\mathcal{H}}_0$, ${\tilde{\mathbf{X}}}_i$ and ${\tilde{\mathbf{Y}}}_i$ are pairwise exchangeable, that is $({\tilde{\mathbf{X}}}_i,{\tilde{\mathbf{Y}}}_i|\mathbf{W},\mathbf{\Theta })\stackrel{d}{=}({\tilde{\mathbf{Y}}}_i,{\tilde{\mathbf{X}}}_i|\mathbf{W},\mathbf{\Theta })$.

The next property is concerned with the conditional independence and pairwise exchangeability between null scores.

Property 3 (Pairwise exchangeability and conditional independence between scores)

Consider scores $s_j^X$ and $s_j^Y$ computed according to Algorithm 1, denote ${\mathbf{s}}_j=\{s_j^X,s_j^Y\}$ for $j\in \mathcal{G}$. Then we have: (a) $\{{\mathbf{s}}_j:j\in \mathcal{G}\}$ are conditionally independent given $\mathbf{W}$ and $\mathbf{\Theta }$; (b) Denote ${\mathbf{S}}_{-i}={\cup }_{j\in \mathcal{G},j\ne i}{\mathbf{s}}_j$. For $i\in {\mathcal{H}}_0$, $s_i^X$ and $s_i^Y$ are pairwise exchangeable given ${\mathbf{S}}_{-i}$ and $\mathbf{\Theta }$: (7) $$\begin{array}{c}(s_i^X,s_i^Y|{\mathbf{S}}_{-i},\mathbf{\Theta })\stackrel{d}{=}(s_i^Y,s_i^X|{\mathbf{S}}_{-i},\mathbf{\Theta }),\text{ or equivalently }\hfill \\ (s_i^X,s_i^Y,{\mathbf{S}}_{-i}|\mathbf{\Theta })\stackrel{d}{=}(s_i^Y,s_i^X,{\mathbf{S}}_{-i}|\mathbf{\Theta }).\hfill \end{array}$$(7)

The pairwise exchangeability (7), which is crucial in proving the FDR theory for PLIS, is a weaker assumption than the joint exchangeability among all null scores $\{s_i^X:i\in {\mathcal{H}}_0\}\cup \{s_i^Y:i\in \mathcal{G}\}$ required in most existing conformal methods. Now we state our main theorem that establishes the validity of PLIS for controlling the FDR.

Theorem 1

. Consider Model (2) and the PLIS procedure outlined in Algorithm 1. Then Properties 1–3 hold and PLIS controls the FDR at level $\alpha$.

Our approach sets itself apart from current theories in conformal inference that rely on joint exchangeability (Yang et al. 2021; Mary and Roquain 2022; Bates et al. 2023; Marandon et al. 2024). In contrast, we introduce a framework that uses only pairwise exchangeability and establishes our theory grounded in the Selective SeqStep + algorithm (Barber and Candès 2015).

2.5 Practical Guidelines and Examples

PLIS uses resampling techniques to generate a mirror process that quantifies the FDP in multiple testing. At its core, PLIS can be situated within the conformal inference framework, which facilitates the development of valid FDR rules that exhibit enhanced robustness to model misspecification. This section presents examples that contextualize PLIS as a valuable tool for “conformalizing” some well-known FDR procedures, producing inferences that possess similar operational characteristics as state-of-the-art conformal methods (Yang et al. 2021; Bates et al. 2023; Marandon et al. 2024; Liang, Sesia, and Sun 2024).

Example 1

(Conformalizing the LIS procedure under HMMs)

A variety of real-world applications involve data with HMM-type structures. Several works (Sun and Cai 2009; Perrot-Dockès et al. 2023) have demonstrated that exploiting the HMM structure can greatly enhance the power of FDR analysis. However, their theoretical analysis necessitates certain conditions to be satisfied, such as homogeneous transition probabilities and identical non-null distributions, which may not be strictly met in practical scenarios. The HMMs are characterized by a set of parameters, denoted as $\vartheta =(\mathbf{A},f_0,f_1)$, where $\mathbf{A}$ is the transition probability matrix, and $f_0$ and $f_1$ are the emission densities for the null and non-null observations, respectively.

In our setup, we use ${\mathcal{M}}_{\text{HM}}$ as the HMM working model, where $\vartheta$ is estimated using the EM algorithm on $\mathbf{W}$. To implement the PLIS procedure, we use the forward-backward algorithm, denoted by ${\mathcal{A}}_{\text{FB}}$, to compute the scores $s_i^X={\mathbb{P}}_{{\mathcal{M}}_{\text{HM}},{\mathcal{A}}_{\text{FB}}}({\theta }_i=0|{\tilde{\mathbf{X}}}_i)$ and $s_i^Y={\mathbb{P}}_{{\mathcal{M}}_{\text{HM}},{\mathcal{A}}_{\text{FB}}}({\theta }_i=0|{\tilde{\mathbf{Y}}}_i)$ for $i\in \mathcal{G}$. The decision-making process is based on Algorithm 1. The PLIS procedure can be regarded as a conformalized adaptation of the LIS procedure proposed by Sun and Cai (2009), computed using ${\mathcal{M}}_{\text{HM}}$ and algorithm ${\mathcal{A}}_{\text{FB}}$. PLIS differs from LIS in that it can handle a variety of data structures that deviate from HMMs. This ability ensures valid FDR control in finite-samples, even when the working model is misspecified or the model parameters are estimated poorly.

Example 2

(Conformalizing the adaptive z-value procedure under the two-group model (3) with independent observations). In the context where observations are assumed to be independent and Efron’s two-group model (3), denoted ${\mathcal{M}}_{\text{TG}}$, is adopted, Sun and Cai (2007) showed that the density ratio (DR) $r(x)=\frac{f_0(x)}{f(x)}$ or the local false discovery rate $\text{Lfdr}(x)=(1-\pi )r(x)$ serves as the optimal building block for FDR analysis. The PLIS framework treats DR and Lfdr equivalently since $(1-\pi )$ is a constant across all study units and only the relative ranking contributes to the operation of PLIS. Using the working model ${\mathcal{M}}_{\text{TG}}$, the conformity score $s_i(·)$ can be taken as the DR function $\widehat{r}(·)=f_0(·)/\widehat{f}(·)$, where $f_0$ is known, and $\widehat{f}(·)$ is the standard kernel density estimator constructed based on $\mathbf{W}$. The PLIS procedure operates by employing $s_i^X\equiv \widehat{r}(X_i)$ and $s_i^Y\equiv \widehat{r}(Y_i)$ in Algorithm 1. As carefully explained in Section D.7 of the supplement, although DR is permutation-invariant with respect to $\mathbf{W}$, the null scores are not jointly exchangeable. Nevertheless, PLIS still works since $(s_i^X,s_i^Y)$ are pairwise exchangeable given $\mathbf{W}$ and $\mathbf{\Theta }$.

Sun and Cai (2007) proposed a class of adaptive z-value (AZ) procedures, and showed that AZ outperforms p-value based methods by adapting to the shape of the alternative distribution. PLIS may be viewed as a conformalized adaptation of AZ, which effectively incorporates the structural information into inference while controlling the FDR in finite samples. This represents a notable advantage over the AZ procedure, which only ensures asymptotic FDR control. We note that the AdaDetect method (Marandon et al. 2024), which also uses the two-group model to estimate a permutation-invariant score function, exhibits similar advantages to PLIS. However, when encountering score functions that lack permutation-invariance, particularly when using the HMM as the working model, PLIS significantly deviates from AdaDetect.

Our customized PLIS procedures remain valid as long as the underlying model belongs to the broad class (2), regardless of the working models or algorithms used. Developing conformity scores using sensible models and computational algorithms is essential for improving the practicality and effectiveness of PLIS. Finally, instances of numerical issues such as significant reductions in power and convergence failures frequently serve as warnings for practitioners, highlighting the unsuitability of the current model for the specific application and emphasizing the need for a new working model.

2.6 Semi-Supervised PLIS

This section considers the scenario in which the null distribution is unknown, and users only have access to n labeled null data points $\mathbf{U}=\{U_1,\dots ,U_n\}$ drawn from that distribution. The objective is to assess whether new data points $\{X_i:i\in \mathcal{G}\}$ follow the same unknown distribution as the reference dataset $\mathbf{U}$. This scenario is known as the semi-supervised multiple testing problem (Mary and Roquain 2022), or the outlier detection problem in the conformal inference literature (Bates et al. 2023; Marandon et al. 2024). We introduce a semi-supervised PLIS (Algorithm 2) to solve this problem.

The algorithm first partitions the labeled null data $\mathbf{U}$ into two subsets: the calibration set $\mathbf{Y}=(Y_i:i\in \mathcal{G})$ and the training set ${\mathbf{U}}^{\mathit{\text{tr}}}=\mathbf{U}\setminus \mathbf{Y}$. The test data $\mathbf{X}$ and calibration data $\mathbf{Y}$ are then used to construct the baseline data $\mathbf{W}$ by following the same steps in Algorithm 1. We move on to the estimation issue. In the case of an HMM, an estimator for the emission distribution $f_0$ can be learned from the labeled null data ${\mathbf{U}}^{\mathit{\text{tr}}}$. Subsequently, we use the baseline data $\mathbf{W}$ to estimate $f_1$ and transition probabilities via the EM algorithm, and then calculate the conformity scores via the forward-backward procedure. If the working model is a two-group model, the conformity scores correspond to the density ratios, which can be estimated using both ${\mathbf{U}}^{\mathit{\text{tr}}}$ and $\mathbf{W}$ via the positive unlabeled (PU) learning algorithms. This approach has also been suggested by Marandon et al. (2024). The above steps are summarized in Algorithm 2.

Algorithm 2

The semi-supervised PLIS procedure

Input: The test data $\mathbf{X}=(X_i:i\in \mathcal{G})$, the null samples $\mathbf{U}=\{U_1,\dots ,U_n\}$, a pre-specified FDR level $\alpha$, a working model $\mathcal{M}$ and a computational algorithm $\mathcal{A}$.

Output: A decision rule $\mathit{\delta }=({\delta }_i:i\in \mathcal{G})$.

1: Split $\mathbf{U}$ into the calibration data $\mathbf{Y}=(Y_i:i\in \mathcal{G})$ and the training data ${\mathbf{U}}^{\mathit{\text{tr}}}=\mathbf{U}\setminus \mathbf{Y}$.

2: Create $\mathbf{W}=(W_i:i\in \mathcal{G})$, where $W_i=h(X_i,Y_i)$ and $h(x,y)$ is given by (5). Construct ${\tilde{\mathbf{X}}}_i$ and ${\tilde{\mathbf{Y}}}_i$ by respectively substituting $X_i$ and $Y_i$ in place of $W_i$ in $\mathbf{W}$ for $i\in \mathcal{G}$.

3: Calculate $s_i^X\equiv {\mathbb{P}}_{\mathcal{M},\mathcal{A}}({\theta }_i=0|{\tilde{\mathbf{X}}}_i)$ and $s_i^Y\equiv {\mathbb{P}}_{\mathcal{M},\mathcal{A}}({\theta }_i=0|{\tilde{\mathbf{Y}}}_i)$ for $i\in \mathcal{G}$ with the user-specified working model $\mathcal{M}$ and algorithm $\mathcal{A}$ based on $\mathbf{W}$ and ${\mathbf{U}}^{\mathit{\text{tr}}}$.

4: Let $\tau =\mathrm{\text{sup}}\{t\in {\mathbf{S}}_X\cup {\mathbf{S}}_Y:Q(t)\le \alpha \}$, where $Q(t)$ is defined in (2.3).

5: Let ${\delta }_i=\mathbb{I}\{s_i^X\le \tau \}$ for $i\in {\mathcal{G}}^r=\{j\in \mathcal{G}:s_j^X<s_j^Y\}$, otherwise ${\delta }_i=0$.

6: Return $\mathit{\delta }=({\delta }_i:i\in \mathcal{G})$.

Remark 2

Here we have assigned exactly one null data point to each test unit, but in situations where there is an abundance of null samples, it may be appropriate to employ the de-randomization idea in Section 3.2 to improve stability. In situations where the size of the null dataset is smaller than m, one may conduct an independent screening procedure to narrow down the focus before using PLIS on a smaller subset of hypotheses. These represent interesting issues for future investigation.

The next theorem, established using de Finetti’s theorem (see Section C in the supplement) and techniques employed in proving Properties 1–3, generalizes Theorem 1 in two ways. First, it considers the semi-supervised setup, where we have access to labeled null data instead of explicit knowledge of the null distribution. Second, it relaxes the assumption of conditional independence given $\Theta$ in Model (2) by allowing for (exchangeable) correlated noise in the data generation process.

Theorem 2

Consider Algorithm 2, if the null data $\{U_1,\dots ,U_n,X_i,i\in {\mathcal{H}}_0\}$ are exchangeable conditional on the non-null data $\{X_i:i\notin {\mathcal{H}}_0\}$ and the true states, then (a) for $i\in {\mathcal{H}}_0$, $s_i^X$ and $s_i^Y$ are pairwise exchangeable, satisfying (7); (b) Algorithm 2 controls FDR at $\alpha$.

In Section E.4 of the supplement, we provide numerical results to corroborate Theorem 2 by considering the semi-supervised setup in which correlated noises are introduced to the observations generated from Model (2).

2.7 Conformal q-values

When the joint exchangeability condition fails to hold, conformal p-values can no longer be properly defined. To address this issue, we introduce the concept of conformal q-value as a significance index to measure the risk associated with individual decisions.

We begin by examining $Q(t)$ in (2.3), which offers a conservative estimate of the FDP. Consider the minimum FDR level $\alpha$ at which $H_{0,i}$ can be “just” rejected: (8) $$q_i\equiv \underset{t\in {\mathbf{S}}_X\cup {\mathbf{S}}_Y,t\ge s_i^X}{\mathrm{\text{min}}}Q(t),\text{ for }i\in {\mathcal{G}}^r;\hspace{1em}{q}_i\equiv 1,\text{ for }i\in \mathcal{G}\setminus {\mathcal{G}}^r.$$(8)

The adjusted $q_i$ is referred to as the conformal q-value corresponding to $H_{0,i}$, owing to its resemblance to the q-value idea introduced by Storey (2003). While Storey’s q-value is constructed based on the empirical distribution of p-values, our conformal q-value is derived from a resampling method and a carefully designed mirror process.

The conformal q-value is a valid and user-friendly significance index that provides clear interpretability for individual decisions, and practitioners can directly use it for decision-making by comparing them with a pre-specified $\alpha$. Theorem 1 and Proposition 1 below together establish the validity of using the conformal q-value (8) in FDR analysis.

Proposition 1

Consider the PLIS procedure as outlined in Algorithm 1 and the conformal q-value defined by (8). Then the decision ${\delta }_i$ presented in Algorithm 1 is equivalent to rejecting $H_{0,i}$ if $q_i\le \alpha$.

3 Connections to Existing Works and Extensions

In this section, we first establish the connection between PLIS and the e-BH method (Section 3.1), then introduce several ensuing extensions: derandomized PLIS (Section 3.2), as well as ${\text{PLIS}}_{\text{cbh}}$ and ${\text{PLIS}}_{\text{sym}}$ (Sections 3.3). Finally, we discuss the distinctions of PLIS from related works (Section 3.4).

3.1 Connection to the e-BH Procedure

In hypothesis testing, an e-value (Vovk and Wang 2021) is defined as the observed value of a nonnegative random variable E that satisfies the condition $\mathbb{E}[E]\le 1$ under the null hypothesis. E-values can be constructed using betting scores (Shafer 2021), likelihood ratios, and stopped super-martingales (Grünwald, de Heide, and Koolen 2020). This section demonstrates how the PLIS framework can be used to construct robust and powerful e-values.

Let $(e_i:1\le i\le m)$ be the e-values for testing $(H_{0,i}:1\le i\le m)$. Wang and Ramdas (2022) proposed the e-BH procedure, which involves first ordering the e-values as $e_{(1)}\ge e_{(2)}\ge \cdots \ge e_{(m)}$, and then choosing a cutoff along the ranking using the following step-wise algorithm. Let $\widehat{k}=\mathrm{\text{max}}\left\{i:(i/m)e_{(i)}\ge (1/\alpha )\right\}$, then reject hypotheses in the set ${\mathcal{R}}_{\mathit{\text{ebh}}}=\{1\le j\le m:e_j\ge e_{(\widehat{k})}\}$. We call $\{e_j:1\le j\le m\}$a set of generalized e-values if (9) $$\mathbb{E}\left(\sum _{j\in {\mathcal{H}}_0}{e}_j\right)\le m.$$(9)

Condition (9) is strictly weaker than the condition that $\mathbb{E}[e_j]\le 1$ for all $j\in {\mathcal{H}}_0$. Wang and Ramdas (2022) proved that if $\{e_j:1\le j\le m\}$ are a set of generalized e-values, then the e-BH procedure controls FDR at $\alpha$ under arbitrary dependence. We can construct a set of generalized e-values based on the PLIS framework. Specifically, define (10) $$e_j=\frac{m{\delta }_j}{1+\sum _{i\in {\mathcal{G}}^c}\mathbb{I}\{s_i^Y\le \tau \}},\hspace{1em}\text{for }j\in \mathcal{G},$$(10) where ${\delta }_j$, ${\mathcal{G}}^c$ and $\tau$ are defined by Algorithm 1. The following theorem can be proved using similar techniques for proving Theorem 1.

Theorem 3

The variables $e_j$ defined in (10) constitute a set of generalized e-values.

The following proposition demonstrates that the PLIS procedure is equivalent to implementing the e-BH procedure with generalized e-values as defined by (10). The connection to e-BH offers valuable insights into our theory, providing a new perspective on why PLIS performs robustly in the face of general dependence and model misspecification.

Proposition 2

If we implement the e-BH procedure with e-values in (10), then the rejection set ${\mathcal{R}}_{\mathit{\text{ebh}}}=\mathcal{R}$, where $\mathcal{R}=\{i\in \mathcal{G}:{\delta }_i=1\}$ is the set of rejected hypotheses by Algorithm 1.

3.2 Derandomized PLIS

While the PLIS framework employs resampling methods, which can introduce additional uncertainties and may be viewed as undesirable by practitioners, the fact that the average of e-values is still an e-value (Vovk and Wang 2021) offers motivation to consider a principled derandomization method. This involves aggregating multiple random samples to reduce variability. Similar ideas have been used in contemporaneous works, including Ren and Barber (2023) and Bashari et al. (2023).

The derandomized PLIS procedure involves running the PLIS procedure repeatedly for N times and averaging the results for decision-making. Specifically, for the kth replication, $k\in [N]$, we generate $\{Y_j^{(k)}:j\in \mathcal{G}\}\stackrel{\text{iid}}{\sim }{f}_0$, and construct ${\tilde{\mathbf{X}}}_i^{(k)}$ and ${\tilde{\mathbf{Y}}}_i^{(k)}$ by respectively substituting $X_i^{(k)}$ and $Y_i^{(k)}$ in place of $W_i^{(k)}$ in ${\mathbf{W}}^{(k)}$, where ${\mathbf{W}}^{(k)}$ is created by combining $\mathbf{X}$ and ${\mathbf{Y}}^{(k)}=(Y_i^{(k)}:i\in \mathcal{G})$ via (5). We then compute the scores $\{s_i^{X(k)},s_i^{Y(k)}:i\in \mathcal{G}\}$ and threshold ${\tau }^{(k)}$ based on Algorithm 1. The corresponding e-values $\{e_i^{(k)}:i\in \mathcal{G}\}$ for replication k can be calculated using (10). This process is repeated for $k=1,\dots ,N$. Finally, we average the individual e-values $e_i^{(k)}$ to compute the summary e-values ${\overline{e}}_i=\frac{1}{N}{\sum }_{k=1}^N{e}_i^{(k)}$ for all $i\in \mathcal{G}$, and apply the e-BH procedure with $\{{\overline{e}}_i:i\in \mathcal{G}\}$. The proposed derandomized PLIS procedure is summarized in Algorithm 3. The theoretical FDR control of derandomized PLIS follows from the validity of the e-BH procedure and Theorem 3.

Remark 3

Algorithm 3 has allowed different replications to have varied ${({\alpha }_k)}_{k=1}^N$. As noted by Ren and Barber (2023), derandomization could lead to higher power if the e-values are carefully aggregated through powerful weighting schemes. We have investigated this issue with some preliminary results provided in Section E.6 in the supplement. The development of optimal weights is an interesting direction for future research.

Algorithm 3

The Derandomized PLIS procedure

Input: Observations $\mathbf{X}=(X_i:i\in \mathcal{G})$, the null distribution $F_0$, working model $\mathcal{M}$, algorithm $\mathcal{A}$, number of replications N, weights ${({\alpha }_k)}_{k=1}^N$, nominal FDR level $\alpha$.

Output: Decisions $\mathit{\delta }=({\delta }_i:i\in \mathcal{G})$.

1: for all $k=1,2,\dots ,N$ do

2: Generate ${\left\{Y_j^{(k)}\right\}}_{j=1}^m\stackrel{\text{iid}}{\sim }{F}_0$. Denote ${\mathbf{Y}}^{(k)}=(Y_j^{(k)}:j\in \mathcal{G})$.

3: Create ${\mathbf{W}}^{(k)}=(W_i^{(k)}:i\in \mathcal{G})$ by combining $\mathbf{X}$ and ${\mathbf{Y}}^{(k)}$ via (5). Construct ${\tilde{\mathbf{X}}}_i^{(k)}$ and ${\tilde{\mathbf{Y}}}_i^{(k)}$ by respectively substituting $X_i$ and $Y_i^{(k)}$ in place of $W_i^{(k)}$ in ${\mathbf{W}}^{(k)}$ for $i\in \mathcal{G}$.

4: Calculate $s_i^{X(k)}\equiv {\mathbb{P}}_{\mathcal{M},\mathcal{A}}({\theta }_i=0|{\tilde{\mathbf{X}}}_i^{(k)})$ and $s_i^{Y(k)}\equiv {\mathbb{P}}_{\mathcal{M},\mathcal{A}}({\theta }_i=0|{\tilde{\mathbf{Y}}}_i^{(k)})$ for $i\in \mathcal{G}$ with user-specified working model $\mathcal{M}$ and algorithm $\mathcal{A}$. Let ${\mathcal{G}}^{r(k)}=\{i\in \mathcal{G}:s_i^{X(k)}<s_i^{Y(k)}\}$, ${\mathcal{G}}^{c(k)}=\{i\in \mathcal{G}:s_i^{Y(k)}<s_i^{X(k)}\}$, and ${\mathbf{S}}_X^{(k)}={\left\{s_i^{X(k)}\right\}}_{i\in \mathcal{G}}$, ${\mathbf{S}}_Y^{(k)}={\left\{s_i^{Y(k)}\right\}}_{i\in \mathcal{G}}$.

5: Let $Q^{(k)}(t)=\frac{1+\sum _{i\in {\mathcal{G}}^{c(k)}}\mathbb{I}\{s_i^{Y(k)}\le t\}}{1\vee (\sum _{i\in {\mathcal{G}}^{r(k)}}\mathbb{I}\{s_i^{X(k)}\le t\})}$ and ${\tau }^{(k)}=\mathrm{\text{sup}}\{t\in {\mathbf{S}}_X^{(k)}\cup {\mathbf{S}}_Y^{(k)}:Q^{(k)}(t)\le {\alpha }_k\}$.

6: Calculate $e_i^{(k)}=\frac{m\mathbb{I}\{s_i^{X(k)}\le {\tau }^{(k)}\}}{1+\sum _{j\in {\mathcal{G}}^{c(k)}}\mathbb{I}\{s_j^{Y(k)}\le {\tau }^{(k)}\}}$ for $i\in {\mathcal{G}}^{r(k)}$, otherwise $e_i^{(k)}=0$.

7: end for

8: Let ${\overline{e}}_i=\frac{1}{N}{\sum }_{k=1}^N{e}_i^{(k)}$ for $i\in \mathcal{G}$. Denote the ordered statistics by ${\overline{e}}_{(1)}\ge {\overline{e}}_{(2)}\ge \cdots \ge {\overline{e}}_{(m)}$. Let $\widehat{k}=\mathrm{\text{max}}\{i:(i{\overline{e}}_{(i)}/m)\ge (1/\alpha )\}$.

9: Let ${\delta }_i=\mathbb{I}\{{\overline{e}}_i\ge {\overline{e}}_{(\widehat{k})}\}$ for $i\in \mathcal{G}$.

Return $\mathit{\delta }=({\delta }_i:i\in \mathcal{G})$

3.3 Comparison with Conformal and Knockoff Methods

PLIS is developed based on the principles of conformal inference, which involve the utilization of machine learning algorithms to calculate conformity scores and the use of calibration data to quantify uncertainties in decision-making. Inspired by an insightful comment from a reviewer, Section D.4 of the supplement shows that the mirror process (2.3), originally adapted from the conformal BH algorithm (see Section D.3), falls into the class of Selective SeqStep + algorithms (Barber and Candès 2015). Our innovative construction of baseline data $\mathbf{W}$, development of pairwise exchangeable scores, and utilization of the new mirror process offer several enhancements to existing conformal methods, enabling the extension of the knockoff framework beyond regression settings (Barber and Candès 2015; Candès et al. 2018; Ren and Candès 2023).

This section explores two variations of PLIS, emphasizing its connection to existing works and providing insights into its advantages.

The first variation replaces both ${\mathcal{G}}^r$ and ${\mathcal{G}}^c$ with $\mathcal{G}$ in Algorithm 1. This variation, denoted as ${\text{PLIS}}_{\text{cbh}}$, essentially involves applying the conformal BH algorithm (Mary and Roquain 2022; Bates et al. 2023) to the conformity scores $\{(s_i^X,s_i^Y):i\in \mathcal{G}\}$ generated by Algorithm 1. However, implementing the conformal BH algorithm with pairwise exchangeable scores may result in overly conservative FDR levels. Additionally, using the adaptive BH algorithm (Bates et al. 2023; Marandon et al. 2024) to overcome this conservativeness is inappropriate, as conformal p-values computed from pairwise exchangeable scores do not satisfy the super-uniformity and PRDS conditions.

The second variation involves applying symmetrization techniques, which have shown success in regression contexts (Barber and Candès 2015; Candès et al. 2018; Ren and Candès 2023; Du et al. 2023). Its key challenge lies in identifying efffective anti-symmetric functions that fulfill $T_j(a,b)=-T_j(b,a)$ to transform pairwise exchangeable scores into symmetrized test statistics. In Section D.4 of the supplement, we demonstrate that the mirror process (2.3) corresponds to using the Selective SeqStep + algorithm proposed by Barber and Candès (2015) with a carefully designed anti-symmetric function. To underscore the significance of the symmetrization process, we consider a popular choice of anti-symmetric function $T_j\equiv s_j^Y-s_j^X$, as suggested in Candès et al. (2018) and their later works. This approach, labeled as ${\text{PLIS}}_{\text{sym}}$, rejects $H_{0,j}$ if $T_j=s_j^Y-s_j^X\ge {\tau }_{\text{sym}}$, where ${\tau }_{\text{sym}}$ is determined by the following FDP process: ${\tau }_{\text{sym}}=\mathrm{\text{inf}}\left\{t\in \left\{\right|T_j|:j\in \mathcal{G}\}:\frac{1+\sum _{j\in \mathcal{G}}\mathbb{I}\{T_j\le -t\}}{\sum _{j\in \mathcal{G}}\mathbb{I}\{T_j\ge t\}}\le \alpha \right\}.$ However, using $T_j=s_j^Y-s_j^X$ may lead to information loss because the contrast has a tendency to cancel out the inherent clustering structure among non-null effects.

Section E.8 presents numerical results that illustrate the advantages of PLIS over its two variations ${\text{PLIS}}_{\text{cbh}}$ and ${\text{PLIS}}_{\text{sym}}$.

3.4 Comparison with Covariate-Adaptive Methods

Structured multiple testing is a subject of extensive research that has attracted significant interest. Our PLIS procedure distinguishes itself from existing covariate-adaptive methods in several aspects. First, in the works of Lei and Fithian (2018), Ignatiadis and Huber (2021) and Ren and Candès (2023), the structural information is represented by a covariate sequence and is incorporated through techniques such as weighted p-values or covariate-adaptive thresholds. By contrast, our approach leverages structural knowledge about the underlying data generation process and integrates it into inference via a user-specified working model. Second, the theory in existing methods is developed based on group-wise cross-weighting (Ignatiadis and Huber 2021) or covariate modulation through data masking (Lei and Fithian 2018; Ren and Candès 2023). In contrast, PLIS operates within the theoretical framework of conformal inference, harnessing the flexibility of machine learning algorithms to learn powerful conformity scores, with associated decision risk being assessed by employing carefully calibrated null scores. Thus, PLIS not only facilitates the effective incorporation of structural information but also guarantees the validity of inference, even in scenarios where the machine learning models are mis-specified. Finally, PLIS offers an intuitive, fast, and user-friendly computational algorithm that efficiently computes scores from the baseline data. This differs from the more computationally intensive algorithms proposed in Lei and Fithian (2018) and Ren and Candès (2023), which involve intricate techniques such as adaptive unmasking and boundary updating.

4 Simulation

This section presents simulation results that compare PLIS, which is implemented as described in the examples of Section 2.5, with several competing methods. For all of the numerical experiments, we set $m=2000$ and $\alpha =0.05$. The observations are generated following this model: $X_i|{\theta }_i\stackrel{\text{ind}.}{\sim }(1-{\theta }_i)\mathcal{N}(0,1)+{\theta }_i\mathcal{N}(\mu ,1)$. The latent states $\Theta$ obey various dependence structures in different experiments. The following methods are considered in our analysis: (a) ${\text{PLIS}}_{\text{HM}},$ ² implemented via Algorithm 1, using the working model ${\mathcal{M}}_{\text{HM}}$ as described in Section 2.5; (b) BH (Benjamini and Hochberg 1995); (c) AdaDetect (Marandon et al. 2024), where the density ratio is estimated by $\widehat{r}(·)=f_0(·)/\widehat{f}(·)$, with $f_0$ being known and $\widehat{f}(·)$ being the standard kernel estimator based on $\{\mathbf{X},\mathbf{Y}\}$. (d) ${\text{PLIS}}_{\text{TG}}$, implemented via Algorithm 1 using the working model ${\mathcal{M}}_{\text{TG}}$ as described in Section 2.5. (e) AdaPT (Lei and Fithian 2018), implemented using the R package adaptMT. Throughout all the simulations conducted, the FDR and Average Power (AP) levels are computed by averaging the results obtained from 200 replications.

4.1 Comparisons with Symmetric Rules

We first consider HMMs in which $\mathbf{\Theta }={({\theta }_i)}_{i=1}^m$ is a binary Markov chain with transition matrix $\mathbf{A}={(a_{\mathit{\text{ij}}})}_{i,j=0,1}={(\mathbb{P}({\theta }_{t+1}=j|{\theta }_t=i))}_{i,j=0,1}$. Here, we fix $a_{00}=0.95$, and the initial state of the latent chain is set to be ${\theta }_1=0$.

We apply BH, AdaDetect, AdaPT, ${\text{PLIS}}_{\text{HM}}$, and ${\text{PLIS}}_{\text{TG}}$ to the simulated data, and summarize the results in Figure 2. Our results demonstrate that all of the methods under consideration effectively control the FDR at the nominal level. Both the ${\text{PLIS}}_{\text{TG}}$ and AdaDetect exhibit similar performance and outperform BH in terms of power. Notably, ${\text{PLIS}}_{\text{HM}}$ exhibits the highest power among all five methods in most scenarios. This can be attributed to the use of asymmetric rules in ${\text{PLIS}}_{\text{HM}}$. AdaPT outperforms/underperforms BH when $a_{11}$ becomes large/small. Additional simulation results for HMMs are presented in Section E.2 of the supplement.

Figure 2: The comparison of FDR and AP for HMMs with $\mu =2.6$ and varying $a_{11}$.

4.2 Comparisons beyond HMMs

This section examines the case where the true transition matrices change over time, which is commonly referred to as a heterogeneous HMM. Such models are particularly useful since the structure of real-world data is often nonstationary. Specifically, we define ${\mathbf{A}}^{(k)}={(a_{\mathit{\text{ij}}}^{(k)})}_{i,j=0,1}$ as the transition probability matrix for the kth transition, where $a_{\mathit{\text{ij}}}^{(k)}=\mathbb{P}({\theta }_k=j|{\theta }_{k-1}=i)$, and $a_{00}^{(k)}=0.95$ is fixed for all k. The conditional distributions remain the same as in the preceding simulations. We consider a scenario in which $a_{11}^{(k)}=0.9e^{-k/1000}$ decreases with k exponentially. This decreasing trend indicates that the stochastic system will gradually stabilize over time, with outliers becoming increasingly rare within the clusters. Since the data in this scenario possess a complex structure, accurately modeling them can be challenging. However, these data do exhibit patterns resembling HMMs, with moderate deviations from homogeneity conditions. Hence, we adopt a homogeneous HMM as a working model to capture the underlying structural patterns in the data and subsequently use the PLIS framework for inference. To assess the performance of PLIS and to compare it with other methods, specifically BH, AdaDetect, and AdaPT, we vary the parameter $\mu$ from 2.2 to 3.2 and apply these methods to simulated data. The results are summarized in Figure 3.

Figure 3: FDR and AP comparison for the heterogeneous HMM with exponentially vanishing $a_{11}$.

Our results demonstrate that all methods are effective in controlling the FDR, with PLIS exhibiting higher average power compared to other methods. This highlights the effectiveness of PLIS in leveraging structural information, even when wrong models and algorithms are employed. Additionally, our analysis of more general data structures, such as two-layer dynamic models and structured models generated from a renewal process, demonstrates that PLIS remains more powerful than the other methods in most cases (see Figures 8–12 in Section E of the supplement). These findings underscore the versatility and robustness of PLIS in a variety of settings, suggesting its considerable potential for analyzing complex data structures in practice.

4.3 Comparisons Under a Covariate-Adaptive Model

Let $\{S_i:i\in \mathcal{G}\}$ denote the covariates that encode side information. Consider a covariate-adaptive model that can be described using a hierarchical approach: ${\theta }_i|(S_i=s)\stackrel{\text{ind}.}{\sim }\text{Bernoulli}({\pi }_s),X_i|(S_i=s,{\theta }_i=0)\stackrel{\text{iid}}{\sim }{F}_0(x),X_i|(S_i=s,{\theta }_i=1)\stackrel{\text{ind}.}{\sim }{F}_{1s}(x),$ where $\mathbb{P}({\theta }_i=1|S_i=s)={\pi }_s$ represents the local sparsity level. An equivalent model, motivated from an empirical Bayes perspective, is given by $X_i|(S_i=s)\stackrel{\text{ind}.}{\sim }{F}_s(x)=(1-{\pi }_s)F_0(x)+{\pi }_s{F}_{1s}(x).$

Let $S_i=i$ for $i\in [3000]$ and fix $F_0=\mathcal{N}(0,1)$. We consider the following scenarios: (i) $F_{1s}=\mathcal{N}(\mu +0.2\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}(0.6s),1)$; ${\pi }_s=0.4(1+\hspace{0.17em}\mathrm{\text{sin}}\hspace{0.17em}(0.2s))$ for $s\in [201,500]\cup [801,1100]\cup [1501,1800]\cup [2101,2400]$, and ${\pi }_s=0.02$ otherwise; (ii) $F_{1s}=\mathcal{N}(2.8,1)$; ${\pi }_s=2\pi$ for $s\in [201,350]\cup [1501,1650]$, ${\pi }_s=\pi$ for $s\in [801,1000]\cup [2101,2300]$, and ${\pi }_s=0.02$ otherwise. Hence, the external covariate $S_i$ captures the local dependence structure, with ${\pi }_s$ and $F_{1s}$ exhibiting similar values when $S_i$ values are comparable.

We apply PLIS (employing the HMM as its working model), LIS (Sun and Cai 2009; using the HMM as its underlying true model), BH and AdaPT (Lei and Fithian 2018; using the R package adaptMT with ${(S_i)}_{i=1}^m$ as its covariate sequence) at FDR level $\alpha =0.05$. The results are presented in Figure 4. The following patterns can be observed. First, AdaPT demonstrates increased efficiency as $S_i$ becomes more informative and may outperform both BH and PLIS in certain settings. Second, LIS, which mistakenly assumes that the data are generated from an HMM, suffers from severe FDR inflation. In contrast, PLIS, despite using the same EM algorithm as LIS, remains effective in FDR control even when the true model significantly diverges from the working model.

Figure 4: FDR and AP comparison under covariate-adaptive models.

5 Application

The identification of eye states has wide-ranging applications across various fields. By monitoring changes in eye states, such as eyelid closures and eye movements, systems can alert drivers when they exhibit signs of drowsiness, thereby preventing potential accidents. In mental health diagnosis, analyzing eye states can offer valuable insights into conditions such as schizophrenia or autism spectrum disorder. Eye state information can also effectively improve facial recognition performance and human–computer interactions.

This section illustrates semi-supervised PLIS for outlier detection using the Electroencephalogram (EEG) Eye State dataset, which is accessible at the UCI Machine Learning Repository (https://archive.ics.uci.edu). By conceptualizing the eye state sequence in a visual task as a Markovian process, the HMM provides a useful tool for modeling eye movement patterns, primarily due to its ability in capturing temporal dependencies within nearby latent states. In recent works such as Chuk et al. (2020) and Hsiao et al. (2022), it is shown that HMMs can effectively represent the dynamic characteristics of eye movement behaviors during a range of cognitive activities, offering valuable insights into the interplay between eye movement patterns and individual differences in cognitive styles or abilities. Furthermore, these studies indicate that the HMM is an appropriate model for describing the local dependence structure within eye states.

During the data preparation process, the eye states of one individual were continuously recorded in chronological order, with detections and manual labeling denoted as 0/1 to signify the eye-open/eye-closed states. The analysis uses data collected from an experiment focusing on measuring the Frontal left Hemisphere (F7). The primary goal is to pinpoint any eye-closed states (outliers) within a sequence of 1500 consecutive records (test samples). After collecting the test dataset, the remaining samples of this measurement with label “0” serve as the null samples. Concretely, we set aside the test samples and proceeded to generate the calibration samples (consisting of 1500 data points) and training samples by splitting the remaining data points with a label of “0”. Subsequently, we applied semi-supervised/conformal multiple testing methods, employing the null samples in both training and calibrating stages. The training data is also used to evaluate the null distribution and calculate p-values when implementing BH and AdaPT.

For the anomaly detection task, we implement various methods at the nominal FDR level of 0.05. The methods being considered include BH, AdaPT, AdaDetect, as well as the proposed semi-supervised PLIS (Algorithm 2) with two distinct working models: the HMM and Efron’s two-group model, respectively denoted as ${\text{PLIS}}_{\text{HM}}$ and ${\text{PLIS}}_{\text{TG}}$. The counts of discovered outliers and corresponding FDPs, computed based on the true labels of test samples and presented in parentheses, are provided for each method below: ${\text{PLIS}}_{\text{HM}}$: 69 (0.014), AdaPT: 38 (0), AdaDetect/${\text{PLIS}}_{\text{TG}}$: 37 (0.054), and BH: 12 (0.083). Notably, the detected outliers are identical for AdaDetect and ${\text{PLIS}}_{\text{TG}}$, corroborating the concluding remark in Example 2 of Section 2.5 regarding the similarity of AdaDetect and ${\text{PLIS}}_{\text{TG}}$ when employing Efron’s two-group model as the working model.

To further examine the comparison results, we visualize the outliers identified by ${\text{PLIS}}_{\text{HM}}$, AdaPT, and AdaDetect in Figure 5. The following patterns can be observed. First, AdaPT shows notable efficacy in detecting large outlier clusters compared to AdaDetect, which may be attributed to its successful integration of structural patterns, as encoded in the external covariate “Time”, in the decision-making process. Second, AdaPT exhibits reduced sensitivity in cases involving smaller clusters, leading to the omission of the final segment where structural information is relatively weak. Third, AdaDetect seems to demonstrate superior effectiveness in identifying isolated outliers. However, it turns out that the discovery at the single location of 711 is indeed a false positive, with the true label being “0”. This aligns with our intuition: if a person closes the eyes, the closed state will persist for a certain period of time. Conversely, by leveraging the structural information, ${\text{PLIS}}_{\text{HM}}$ and AdaPT successfully filter out this seemingly plausible signal. Finally, ${\text{PLIS}}_{\text{HM}}$ outperforms both AdaPT and AdaDetect by detecting the highest number of outliers while keeping the FDP well below the nominal level. This is made possible through the utilization of the HMM working model, which effectively exploits the informative structures within the data.

Figure 5: The discoveries of ${\text{PLIS}}_{\text{HM}}$ (red, top), AdaPT (green, middle) and AdaDetect (blue, bottom) in detecting eye-closed states in 1500 consecutive records according to F7 measurements.

Supplementary Materials

The supplement contains the proofs of theorems and propositions, some extensions of the proposed PLIS procedure, a comparative review of PLIS and related methods, and additional numerical results.

Acknowledgments

The authors thank the associate editor and two referees for their valuable comments, which have greatly enhanced the article. The comment from an insightful referee on the Selective SeqStep + algorithm and anti-symmetric functions has contributed significantly to enhance our understanding of related works.

Disclosure Statement

No potential conflict of interest was reported by the author(s).

Notes

1 $(p_1,\dots ,p_m)$ have the positive regression dependency on each subset (PRDS) property on ${\mathcal{H}}_0$ if $\mathbb{P}\left\{(p_j,j\in [m])\in D|p_i=u\right\}$ is nondecreasing in u for any $i\in {\mathcal{H}}_0$ and any nondecreasing set $D\subset {\mathbb{R}}^m$.

2 In subsequent simulations, ${\text{PLIS}}_{\text{HM}}$ is referred to as PLIS if ${\text{PLIS}}_{\text{TG}}$ is not involved.