Why Compare Procedures Under Block Dependence?
Genomic tests are often correlated because nearby variants, genes, or other features share biological or technical structure. A useful way to study this dependence is to arrange tests into blocks: tests within a block are correlated, while different blocks may be independent. The question here is practical and specific: when we apply the Benjamini–Hochberg (BH) and Benjamini–Yekutieli (BY) procedures to such simulated data, how do their false discovery rates and discovery counts compare?
A simulation can show how procedures behave under a stated model. It cannot, by itself, prove a general error-control guarantee. Keep that distinction in view: BH has an established FDR control result under independence and under positive regression dependence on a subset (PRDS), while BY provides control under arbitrary dependence when the individual null p-values are valid. Block correlation gives a concrete setting in which to compare their behavior, not a replacement for those conditions.
BH and BY Thresholds
Write the ordered p-values as \(p_{(1)}\le \cdots \le p_{(m)}\). At significance level \(\alpha\), BH finds the largest rank \(k\) for which \(p_{(k)}\le \alpha k/m\), then rejects the hypotheses with p-values no larger than \(p_{(k)}\). If no rank meets the condition, it makes no rejections.
BY uses the same step-up rule after reducing the rank-specific critical values by the harmonic factor \(c_m=\sum_{j=1}^m 1/j\). Thus its critical value at rank \(i\) is \(\alpha i/(m c_m)\). For \(m>1\), \(c_m>1\), so the BY critical values are strictly smaller than BH’s. For \(m=1\), \(c_1=1\), and the two procedures have identical critical values. Smaller critical values can lead BY to make fewer discoveries, although the actual counts depend on the p-values.
Worked Example: Comparing BH and BY on One Ranked List
Suppose six tests have ordered p-values \(0.002,\ 0.009,\ 0.018,\ 0.040,\ 0.200,\ 0.700\), and set \(\alpha=0.05\). The values are hypothetical. Here \(c_6=1+1/2+1/3+1/4+1/5+1/6=2.45\).
The BH critical values at ranks 1 through 6 are \(0.00833,\ 0.01667,\ 0.02500,\ 0.03333,\ 0.04167,\ 0.05000\), since the rank-\(i\) value is \(0.05i/6\). The first three comparisons pass: \(0.002\le0.00833\), \(0.009\le0.01667\), and \(0.018\le0.02500\). Rank 4 fails because \(0.040>0.03333\), and ranks 5 and 6 fail as well. The largest passing rank is 3, so BH rejects the first three hypotheses.
For BY, divide each BH critical value by \(2.45\). At ranks 1 through 6 the values are approximately \(0.00340,\ 0.00680,\ 0.01020,\ 0.01361,\ 0.01701,\ 0.02041\). Rank 1 passes, since \(0.002\le0.00340\). Rank 2 fails, and the p-values at ranks 3 through 6 also exceed their BY critical values. BY therefore rejects only the first hypothesis. This example shows the consequence of the more stringent thresholds for this list; it does not establish what either procedure will do on every data set.
A Block Simulation With Valid Null P-Values
A simulation comparison is meaningful only if its data-generating assumptions are explicit. One illustrative design has \(m=1000\) tests in 50 blocks of 20. In every block, 16 hypotheses are null and 4 are non-null, giving 800 nulls and 200 non-nulls overall. The same placement of null and non-null tests is used in each repetition.
For test \(j\) in block \(b\), generate a test statistic
where all \(U_b\) and \(\varepsilon_{bj}\) are mutually independent standard normal random variables. Set \(\delta_{bj}=0\) for null tests and \(\delta_{bj}=2.8\) for non-null tests. Simulate three settings, \(\rho=0,\ 0.3,\ 0.7\). Within a block, distinct statistics have correlation \(\rho\); statistics in different blocks are independent. The variance of each statistic is \(\rho+(1-\rho)=1\). Under a null, \(Z_{bj}\sim N(0,1)\), so the one-sided p-value \(p_{bj}=1-\Phi(Z_{bj})\) is uniform on \([0,1]\). Under a non-null, the statistic has a positive mean shift. This establishes valid marginal null p-values in the design; it does not make tests within a block independent.
For each simulated data set, apply BH and BY at \(\alpha=0.05\). The BH result can be justified under independence and, for appropriate one-sided Gaussian testing with nonnegative dependence, PRDS. BY is included as the procedure with arbitrary-dependence FDR control for valid marginal null p-values. The simulation itself is not the proof of either result.
Reading Simulation Results
The following hypothetical output illustrates how results might be summarized over \(B=5000\) repetitions. The figures are illustrative rather than a claim about an observed genomic study. Empirical FDR is averaged across repetitions, and mean discoveries are averaged across the same repetitions.
| Within-block correlation \(\rho\) | Procedure | Empirical FDR | Mean discoveries |
|---|---|---|---|
| 0 | BH | 0.049 | 96.2 |
| 0 | BY | 0.008 | 66.4 |
| 0.3 | BH | 0.047 | 93.8 |
| 0.3 | BY | 0.007 | 63.7 |
| 0.7 | BH | 0.043 | 88.9 |
| 0.7 | BY | 0.005 | 58.2 |
In this illustrative run, BH has more discoveries and an empirical FDR near, but below, the target \(0.05\) in all three settings. BY has fewer discoveries and a lower empirical FDR. One possible reason for the gap is that the harmonic factor for \(m=1000\) is about \(c_{1000}=7.485\), substantially reducing BY’s critical values. The differences in the table are outcomes of the specified simulation, not universal numerical effects of block correlation.
Notice also what the mean discovery count does not say: it is not the number of discoveries in every repetition, nor does it indicate how many are true. A fuller report can include mean true discoveries \(E(R-V)\), the frequency of zero rejections, and the distribution of FDP across repetitions. The empirical FDR is an average; it can conceal variability between simulated data sets.
Worked Example: Calculating Empirical FDR From Repetitions
Suppose five hypothetical repetitions of one procedure produce these pairs \((V,R)\): \((2,20),\ (0,15),\ (1,10),\ (0,0),\ (3,30)\). The FDP values are \(2/20=0.10,\ 0/15=0,\ 1/10=0.10,\ 0/\max(0,1)=0,\) and \(3/30=0.10\).
The empirical FDR is their average:
The discovery counts sum to \(20+15+10+0+30=75\), so the mean discovery count is \(75/5=15.0\). The total number of true discoveries is \(18+15+9+0+27=69\), giving a mean of \(69/5=13.8\). These three summaries answer different questions: average false discovery proportion, average total discoveries, and average true discoveries.
Checks, Failure Modes, and Alternatives
Before interpreting a simulation, check the design and its implementation. Under the null, the generated p-values should be marginally uniform; their histogram or empirical quantiles can reveal a coding error or an incorrectly calibrated test. Check that simulated within-block correlations are near the intended \(\rho\), between-block correlations are near zero, and that the intended number of null and non-null tests is used in every repetition. Confirm that each procedure uses the ordered p-values and the same \(\alpha\), and that \(V\) is counted using the known simulated truth.
Monte Carlo variation matters. A different random run can produce different empirical FDR estimates and discovery counts. Increasing \(B\) generally stabilizes the averages, but does not repair a misspecified dependence model or invalid null p-values. A useful report states the number of repetitions and the data-generating assumptions; when estimates are close to a decision boundary, report Monte Carlo uncertainty rather than treating rounded values as exact.
The main choice in this comparison is between BH and BY. BH is less conservative and often yields more discoveries when its dependence conditions apply. BY has smaller critical values for \(m>1\), often yielding fewer discoveries, but its general dependence guarantee requires valid marginal null p-values. If \(m=1\), the harmonic factor is exactly one, so neither threshold is more stringent. Neither procedure repairs p-values that are invalid under the null, and a simulation under one block structure cannot represent every genomic dependence pattern.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- For \(m=1\), what is \(c_m\), and how do the BH and BY critical values compare?
- In a simulation with four repetitions, the FDP values are \(0.2,\ 0,\ 0.1,\ 0\). What is the empirical FDR?
- Why does a valid uniform marginal p-value for each null test not imply independence among tests in a block?
- A simulation reports mean discoveries but does not report \(V\). Can it estimate empirical FDR from that summary alone? Explain.
- Why should a low empirical FDR in one block simulation not be presented as proof that BH controls FDR under every dependence pattern?
Solutions. (1) \(c_1=1\), so the BH and BY critical values are identical. (2) The empirical FDR is \((0.2+0+0.1+0)/4=0.075\). (3) Marginal validity describes each null p-value’s distribution individually; shared block random effects can make the values jointly dependent. (4) No. FDP requires the number of false rejections \(V\), as well as \(R\), in each repetition; mean \(R\) alone does not supply it. (5) The result is conditional on the simulated model and finite number of repetitions. It cannot establish behavior under dependence structures that were not simulated or replace the assumptions of a theoretical guarantee.