The Benjamini--Hochberg Procedure Can Fail to Control the FDR for Correlated Two-Sided Gaussian Tests

This paper kills a 20-year-old belief at the foundations of multiple testing: the Benjamini–Hochberg (BH) procedure — the default tool for controlling false discoveries in genomics, neuroimaging, A/B testing, everywhere — does not always control the false discovery rate for correlated two-sided Gaussian tests. The counterexample is an explicit three-block factor model, certified by rigorous interval arithmetic to give FDR > 0.0104 at nominal level α = 0.01 as the number of tests grows. The violation is tiny but real. And the meta-story is arguably as notable as the math: the counterexample, the proof strategy, and the numerical certificate were produced by GPT-5.6 Pro in about 90 minutes of reasoning, then checked line-by-line by the author.

The conjecture everyone was leaning on

Quick setup. You test m hypotheses and get p-values P_1,\dots,P_m. BH sorts them, finds the largest r such that P_{(r)} \le \alpha r/m, and rejects the r smallest. The false discovery proportion is \mathrm{FDP} = V/\max(R,1), where R is the number of rejections and V the number of those that are true nulls; the FDR is \mathbb{E}[\mathrm{FDP}].

The classical guarantees form a ladder. Under independence, BH satisfies \mathrm{FDR} \le \pi_0\alpha, where \pi_0 is the fraction of true nulls — so it's actually conservative. The same bound holds under positive regression dependence (PRDS, Benjamini–Yekutieli 2001), which covers one-sided Gaussian tests with nonnegative correlations. Under completely arbitrary dependence, the guarantee degrades badly: you must shrink the level to \alpha/H_m, where H_m = \sum_{r=1}^m 1/r \approx \ln m — a log-factor penalty that makes the BY-corrected procedure painfully conservative, so almost nobody uses it.

The practically crucial middle case: two-sided Gaussian tests with arbitrary correlation. This is what people actually run — z-tests on correlated quantities (neighboring genes, adjacent voxels) where the sign of the effect is unknown. Twenty years of simulations and partial theory (Farcomeni 2006, Reiner-Benaim 2007, Kim & van de Wiel 2008) found no violation. Benjamini himself wrote in 2010 that "convincing simutheoretical evidence indicates that [FDR control] holds for two-sided z-tests with any correlation structure." As recently as 2023, Sarkar called proving this conjecture "an urgent and important undertaking."

It's false.

Why two-sided tests were always the suspicious case

The PRDS machinery for one-sided tests rests on a monotone coupling: if positively correlated nulls drift toward rejection, everything drifts together, and this positive dependence turns out to work in BH's favor. Two-sided tests break the monotonicity in a specific way. The rejection event is

a union of two tails. Conditioning on a correlated coordinate landing in the right tail shifts X_i one way; the left tail shifts it the other way. Folding X into |X| destroys the monotone-regression structure that the one-sided proofs need. Everyone knew the proofs didn't transfer; the question was whether control failed or merely the proof technique did.

Here is the sharper quantitative reason folding is dangerous. If X_i \sim \mathcal{N}(\mu_i, 1) conditional on some latent variable, and the conditional mean is any \mu_i \ne 0, then the two-sided tail \mathbb{P}(|X_i| \ge c) is larger than under \mu_i = 0 — regardless of the sign of \mu_i. So a common factor that jitters null means around zero makes null p-values conditionally super-uniform (enriched near zero) for every nonzero factor value, positive or negative. One-sided tests only get this enrichment on one side. That asymmetry between the one-sided and two-sided worlds is the crack the counterexample drives a wedge into.

The counterexample: a see-saw with one hinge

Take a single latent factor Z \sim \mathcal{N}(0,1) and m = 100N coordinates in three blocks:

The nulls load +0.3 on Z; the two signal blocks load -0.3 and -0.72. Conditional on Z = z, the null block has mean 0.3z (folded tails inflated for any z \ne 0), while the signal means are 2.4 - 0.3z and 4.4 - 0.72z.

Now watch the see-saw. When z is negative, the signals get stronger, producing more rejections — which raises the data-dependent BH threshold t_R = \alpha R/m, letting more nulls through — at the very same moment the nulls' folded tails are inflated by |0.3z|. When z is positive, signals weaken somewhat but nulls inflate again (folding doesn't care about sign). There is no factor value at which the nulls behave better than uniform.

How much inflation is needed? Here's the accounting that makes the whole thing click. With continuous p-values, the FDP at the BH threshold \tau satisfies the exact identity

where \widehat F_0 is the empirical CDF of the null p-values. Under independence, \widehat F_0(\tau) \approx \tau and you get the familiar \pi_0 \alpha = 0.96 \times 0.01 = 0.0096. To break the level 0.01, the nulls need to be super-uniform at the realized threshold by more than a factor 1/\pi_0 \approx 1.042 — on average over Z, at a threshold that itself depends on Z. The loadings and means above are tuned exactly so this self-consistent enrichment holds, clearing the bar with about 4% to spare: certified FDR > 0.0104.

Bracket, don't solve

The proof architecture is worth understanding because it's a template for computer-certified asymptotic statistics.

Conditional on Z = z, the three blocks are i.i.d. within themselves, so by Glivenko–Cantelli the empirical p-value CDF converges uniformly to an explicit deterministic mixture G_z. BH has a classical geometric reading: its threshold \tau is the rightmost point where the empirical CDF sits above the line t/\alpha. The usual asymptotic approach solves for the limiting crossing point — but that requires uniqueness, transversality, all fragile to certify numerically.

The paper's Lemma 2 sidesteps this entirely: you only need two strict sign conditions. Find one point v where G_z(v) > v/\alpha (a feasible threshold exists at least this far out) and a band [w, \alpha] where G_z(t) < t/\alpha strictly (no threshold survives beyond w). Then v \le \liminf \tau_N \le \limsup \tau_N \le w, and Lemma 3 converts the bracket into a conditional FDP lower bound \alpha \pi_0 F_{0,z}(v)/w.

That reduction is what makes machine certification possible. Both conditional tail functions are monotone — decreasing in the cutoff c, increasing in the absolute conditional mean — so checking the sign conditions over a continuum of (z, c) values reduces to finitely many axis-aligned rectangles: 1,000 bins for z \in [-5,5] crossed with a rational grid of ~7,400 cutoff values. Every input is an exact rational; every Gaussian tail is evaluated with Arb's outward-rounded ball arithmetic, and a strict comparison is accepted only when the entire ball lies on the claimed side of zero. Per-bin FDP lower bounds d_k are then integrated against the Gaussian measure of Z, and Fatou's lemma turns pointwise liminf bounds into the FDR bound. The certified total:

Notice what this chart rules out: the failure is not driven by rare tail values of the factor. The bulk of the certified excess comes from perfectly ordinary |Z| \le 2. The construction is a knife-edge in effect size (0.0104 vs 0.01), but not a rare-event pathology.

Does it show up at finite m?

The theorem is asymptotic ("for all sufficiently large N"), so the paper runs a stratified Monte Carlo — stratifying on Z into 1,000 equiprobable slices, which keeps the estimand exactly the FDR while slashing variance — with 100,000 full datasets per dimension:

At m = 20{,}000 the excess is about 3.6% of the nominal level and clearly resolved. At m = 5{,}000 it isn't — either the simulation needs more replications or the finite-sample FDR genuinely doesn't exceed \alpha there yet. Whether BH provably controls FDR for small m under arbitrary two-sided Gaussian dependence remains open.

The GPT-5.6 Pro angle

The model was given only the mathematical definition of the BH procedure and asked to prove or disprove the conjecture. After ~90 minutes it produced the counterexample, the bracketing proof strategy, and the code for the numerical certificate; Dobriban then verified everything and edited the draft (the ChatGPT conversation is shared publicly, and the certificate code is in the paper's Appendix B and on GitHub).

Two things are striking. First, the model went against twenty years of expert consensus and "simutheoretical evidence" — it didn't pattern-match the literature's expected answer. Second, the proof mode is telling: a disproof-by-construction with a machine-checkable interval-arithmetic certificate plays directly to what these models are good at — searching parameter space, writing verification code — while the human's job compresses to auditing a finite, fully explicit artifact. Expect more open problems whose resolution has this shape (adversarial construction + rigorous numerics) to fall soon; they may be the most AI-tractable class in mathematics.

What changes, and what to be skeptical of

What changes theoretically: the two-sided Gaussian case now provably sits strictly between the PRDS-safe regime and the arbitrary-dependence regime — BH can exceed its level, so the conjecture-driven hope of a clean universal proof is dead, and the dependence-adjusted procedures people have been developing (Fithian–Lei 2022, Sarkar's shifted BH variants) have concrete rather than merely hypothetical motivation. The live question shifts to: how big can the inflation be? The certified excess here is ~4% relative; the paper reports that additional searches over related models found similarly small violations. Meanwhile the only universal upper bound remains the harmonic factor H_m \approx \ln m — roughly 10× at m = 20{,}000. The gap between "4% inflation exists" and "10× inflation is possible" is enormous, and a universal inflation bound would be a far more consequential result than this one. The paper explicitly flags it as open.

What to be skeptical of, practically: if you run BH at \alpha = 0.05 on real data, a few percent relative inflation is dwarfed by every other violated assumption in your pipeline. The construction is adversarial — signals anti-correlated with nulls through a common factor, magnitudes tuned to a knife edge — though note that shared latent factors with opposite-sign loadings (batch effects, market factors) are not exotic in real data, and this construction shows the certified excess accrues at typical factor values, not in the tails. Also, the rigorous result is asymptotic; the finite-m evidence is Monte Carlo, not certified.

Where to spend your time: Sections 4–5 are the heart — the conditional mixture CDF G_z, the see-saw tuning, and the rectangle reduction (Figure 3 in the paper) — and they're short. If you're mainly here for the AI story, the linked ChatGPT conversation is a primary source of a kind we'll be seeing much more of.