A robust and powerful method for assessing replicability of high dimensional data

1. Problem Statement

The paper addresses the challenge of identifying signals that replicate across multiple independent studies in high-dimensional settings (e.g., Genome-Wide Association Studies with thousands to millions of hypotheses).

• The Core Challenge: The "replicability null hypothesis" is composite. A feature is considered non-replicable if it is a null effect in at least one of the $n$ studies.
• Limitations of Existing Methods:
  • Ad hoc approaches (e.g., applying Benjamini-Hochberg separately to each study and taking the intersection) fail to control the False Discovery Rate (FDR) under the composite null.
  • Conservative methods (e.g., MaxP) guarantee FDR control but suffer from low statistical power.
  • Parametric/Bayesian methods often rely on restrictive modeling assumptions or require tuning parameters.
  • Scalability: Directly modeling $n$ studies requires estimating $2^n$ latent configurations, which is computationally infeasible for large $n$.
  • Heterogeneity: Many methods assume a common signal configuration or ignore between-study heterogeneity in non-null distributions.

2. Methodology

The authors propose a general Empirical Bayes framework that jointly models summary-level $p$-values while explicitly accounting for between-study heterogeneity.

A. Two-Study Framework ($n=2$)

1. Model Setup:
   • Let $\theta_{ij} \in \{0, 1\}$ be the latent state (null/non-null) for hypothesis $i$ in study $j$.
   • The joint latent configuration $\tau_i = (\theta_{i1}, \theta_{i2})$ follows a multinomial distribution with probabilities $\xi_{uv}$.
   • The Replicability Null ($H_0$) is defined as $\theta_{i1}\theta_{i2} = 0$ (i.e., configurations $(0,0), (0,1), (1,0)$).
   • The $p$-values follow a two-group mixture model where the null density $f_0$ is uniform, and non-null densities $f_1, f_2$ are study-specific.
2. Nonparametric Density Estimation:
   • Instead of assuming parametric forms for $f_1$ and $f_2$, the authors estimate them nonparametrically under a monotonicity constraint (non-increasing density ratios $f_k/f_0$). This ensures robustness against model misspecification and eliminates the need for tuning parameters.
   • Estimation is performed using the EM algorithm combined with the Pool-Adjacent-Violators Algorithm (PAVA).
3. Test Statistic:
   • The Local False Discovery Rate (Lfdr) is used as the test statistic, defined as the posterior probability that a feature is non-replicable:
     $$\text{Lfdr}(p_1, p_2) = P(\theta_{i1}\theta_{i2} = 0 \mid p_1, p_2)$$
4. Thresholding:
   • A data-driven threshold $\hat{\lambda}$ is selected to control the FDR at a nominal level $\alpha$ using a step-up procedure based on the estimated Lfdr.

B. Multi-Study Extension ($n > 2$)

To overcome the computational infeasibility of estimating $2^n$ parameters:

1. Pairwise Rejection Strategy:
   • The method decomposes the exponentially large composite null space into disjoint components based on the number of null studies.
   • It utilizes pairwise Lfdr statistics ($Lfdr_{jj'}$) for all study pairs $(j, j')$.
   • A feature is declared replicable only if all pairwise Lfdrs fall below their respective thresholds.
2. Scalability:
   • This approach reduces the computational complexity from exponential ($O(2^n)$) to quadratic ($O(n^2)$) with respect to the number of studies.
3. Parameter Estimation:
   • The mixture proportions (prior probabilities of latent configurations) are estimated using a method of moments approach based on the probability of observing small $p$-values simultaneously across studies, leveraging the fact that true signals tend to have small $p$-values in all studies.

3. Key Theoretical Contributions

The paper provides rigorous theoretical guarantees for the proposed method:

• Identifiability: Proves that the model parameters and nonparametric density functions are uniquely identifiable under mild conditions (specifically, that non-null densities vanish at $p=1$).
• Convergence Rates:
  • Establishes cubic-rate convergence ($O(m^{-1/3})$) for the nonparametric maximum likelihood estimator (NPMLE) of the joint $p$-value density under the Hellinger distance.
  • Proves minimax optimality, showing that the proposed estimator achieves the optimal worst-case risk rate of $O(m^{-2/3})$ for squared Hellinger loss.
• FDR Control: Proves that the proposed procedure asymptotically controls the FDR at the prespecified level $\alpha$ under mild regularity conditions.
• Optimal Power: Demonstrates that the Lfdr-based thresholding rule is asymptotically power-optimal among all procedures controlling the marginal FDR.

4. Results

Simulation Studies

Extensive simulations were conducted comparing the proposed method against state-of-the-art alternatives (Ad hoc BH, MaxP, IDR, MaRR, radjust, AdaFilter, JUMP).

• FDR Control: The proposed method maintained valid FDR control across independent and dependent (correlated) data scenarios, whereas methods like Ad hoc BH and IDR failed to control FDR under dependence.
• Statistical Power: The proposed method consistently achieved the highest statistical power across various scenarios (sparse vs. dense signals, varying effect sizes, and different correlation structures). It significantly outperformed conservative methods like MaxP and radjust.
• Robustness: The method remained robust under local dependence structures (block correlations), a common feature in genomic data.

Real Data Application

The method was applied to Type 2 Diabetes (T2D) GWAS data from European and East Asian cohorts.

• Findings: The proposed method identified 622 replicable SNPs, including 280 SNPs that were not detected by any competing method.
• Biological Validation: Many of the uniquely identified SNPs mapped to genes with established links to T2D (e.g., GCKR, BCL11A, GRB14), validating the biological relevance of the new discoveries.
• Comparison: Competing methods were either too conservative (missing true signals) or failed to detect the specific cross-ancestry replicability patterns captured by the proposed framework.

5. Significance and Impact

• Methodological Advancement: This work bridges the gap between robust nonparametric inference and high-dimensional replicability analysis. By avoiding parametric assumptions and tuning parameters, it offers a "plug-and-play" solution for diverse data types.
• Scalability: The pairwise decomposition strategy makes rigorous replicability analysis feasible for large-scale meta-analyses involving many studies (e.g., global consortia), a task previously computationally prohibitive.
• Scientific Utility: The application to T2D GWAS demonstrates the method's ability to uncover biologically meaningful signals that other methods miss, directly addressing the "replicability crisis" by providing a statistically rigorous framework for cross-study validation.
• Theoretical Foundation: The establishment of minimax optimality and cubic-rate convergence for nonparametric MLEs in this specific context sets a new theoretical benchmark for high-dimensional inference.

In summary, the paper presents a robust, powerful, and scalable Empirical Bayes framework that outperforms existing methods in both theoretical properties and practical application, offering a new standard for assessing replicability in high-dimensional biomedical research.