FDR Control Under Dependence

Why Dependence Matters for False-Discovery Control

A genome-wide association analysis may test many variants for association with a phenotype. The resulting p-values are not necessarily independent: variants in linkage disequilibrium can carry overlapping information, and tests may also share participants or correlated phenotypes. Dependence does not make every multiple-testing procedure invalid, but it matters because a procedure’s error guarantee depends on how the test statistics relate to one another.

The Benjamini–Hochberg procedure controls the false discovery rate under independence and under some forms of positive dependence. When that structure cannot be justified, the Benjamini–Yekutieli procedure provides a guarantee that permits arbitrary dependence. It does so by making the Benjamini–Hochberg thresholds more stringent.

Definition: Suppose a family contains $m$ hypothesis tests. Let $R$ be the number of rejected hypotheses and $V$ the number of rejected hypotheses that are true nulls. The false discovery proportion is $V/\max(R,1)$, and the false discovery rate (FDR) is its expected value: $$ \operatorname{FDR}=E\left[\frac{V}{\max(R,1)}\right]. $$

FDR is an expectation over repeated applications of the testing procedure. It is not the probability that any particular rejected association is false, nor is it the probability that the null hypothesis is true. FDR control at level $\alpha$ means that the expected fraction of false discoveries among the rejections is at most $\alpha$, under the procedure’s assumptions. By contrast, family-wise error concerns the probability of making at least one false rejection.

The Benjamini–Yekutieli Step-Up Procedure

Write the p-values in increasing order as $p_{(1)}\le p_{(2)}\le\cdots\le p_{(m)}$. Let $\alpha$ be the chosen FDR level, and define the $m$th harmonic number $c_m=\sum_{j=1}^{m}1/j$. The Benjamini–Yekutieli procedure scales the Benjamini–Hochberg thresholds by $c_m$.

Definition: For $k=1,\ldots,m$, the Benjamini–Yekutieli step-up threshold is $$ t_k=\frac{k\alpha}{m c_m},\qquad c_m=\sum_{j=1}^{m}\frac{1}{j}. $$ Find the largest $k$ for which $p_{(k)}\le t_k$. If such a $k$ exists, reject the hypotheses corresponding to $p_{(1)},\ldots,p_{(k)}$. If no p-value meets its threshold, reject none.

The formal guarantee assumes that each true-null p-value is valid: under its null hypothesis, it is super-uniform, meaning $P(p_i\le u)\le u$ for every $0\le u\le1$. The tests can otherwise have arbitrary dependence. The family of $m$ tests and the target level $\alpha$ must be specified for the analysis. In a genetic association setting, dependence among variants does not itself prevent use of the procedure, provided the individual null p-values are valid.

Worked Example: Applying the Procedure to Correlated Variant Tests

Suppose a hypothetical GWAS tests $m=6$ variants. The variants may be correlated, but assume each true-null p-value is valid. The ordered p-values are $0.0004, 0.0015, 0.0040, 0.0090, 0.0300,$ and $0.2000$. Apply the Benjamini–Yekutieli procedure at $\alpha=0.05$.

First calculate the harmonic factor: $c_6=1+1/2+1/3+1/4+1/5+1/6=2.450$. Thus the threshold for rank $k$ is $k(0.05)/(6\times2.450)$. The common multiplier is $0.05/14.7=0.003401$, so the rank-specific thresholds are:

(t_1,\ldots,t_6)=(0.003401,\ 0.006803,\ 0.01020,\ 0.01361,\ 0.01701,\ 0.02041).

Compare each ordered p-value with its threshold. The first three comparisons pass: $0.0004\le0.003401$, $0.0015\le0.006803$, and $0.0040\le0.01020$. The fourth does too: $0.0090\le0.01361$. The fifth and sixth fail: $0.0300>0.01701$ and $0.2000>0.02041$. The largest passing rank is therefore $k=4$, so the procedure rejects the four hypotheses with the smallest p-values. Under the assumptions above, the procedure controls the FDR at no more than 0.05 despite dependence among the tests.

Why the Guarantee Holds

The harmonic factor is not just a numerical precaution: it allows an upper bound on the expected false-discovery proportion without assuming independence. A short derivation shows how the adjustment enters. Let $I_0$ be the set of true-null tests, with size $m_0$, and let $R$ be the number of rejections. If the procedure rejects $R=r>0$ hypotheses, every rejected p-value is at most $t_r=\alpha r/(m c_m)$. The false discoveries can therefore be counted through the events that each true-null p-value falls below the threshold determined by $R$.

Fix one true-null test $i$. Define $A_r=\{R=r,\ p_i\le t_r\}$, for $r=1,\ldots,m$. These events are disjoint. For any $r$, if $A_s$ occurs for some $s\le r$, then $p_i\le t_s\le t_r$. Consequently, the cumulative probability $Q_r=\sum_{s=1}^{r}P(A_s)$ satisfies $Q_r\le P(p_i\le t_r)\le t_r$, where the last inequality is the valid-p-value condition.

The expected contribution of this true null to the false discovery proportion is bounded by a weighted sum. Summation by parts gives the following expression and bound:

\begin{aligned} E\left[\frac{\mathbf{1}\{i\text{ is rejected}\}}{\max(R,1)}\right] &=\sum_{r=1}^{m}\frac{P(A_r)}{r}\\ &=\frac{Q_m}{m}+\sum_{r=1}^{m-1}Q_r\left(\frac{1}{r}-\frac{1}{r+1}\right)\\ &\le \frac{t_m}{m}+\sum_{r=1}^{m-1}\frac{t_r}{r(r+1)}\\ &=\frac{\alpha}{m c_m}\left(1+\sum_{r=1}^{m-1}\frac{1}{r+1}\right) =\frac{\alpha}{m}. \end{aligned}

The equality in the last line uses $t_r=\alpha r/(m c_m)$ and $c_m=1+1/2+\cdots+1/m$. Summing the bound over the $m_0$ true-null tests yields $\operatorname{FDR}\le m_0\alpha/m\le\alpha$. This argument uses only the marginal validity of each true-null p-value; it does not require the p-values to be independent.

Worked Comparisons and Adjusted P-Values

Worked Example: Comparing Benjamini–Hochberg and Benjamini–Yekutieli

Consider a hypothetical family of $m=5$ valid association p-values, ordered as $0.001, 0.004, 0.008, 0.020,$ and $0.040$, with $\alpha=0.05$. The Benjamini–Hochberg thresholds are $k\alpha/m$, namely $0.010, 0.020, 0.030, 0.040,$ and $0.050$. Each p-value is at or below its corresponding threshold, so this procedure rejects all five.

For Benjamini–Yekutieli, $c_5=1+1/2+1/3+1/4+1/5=2.2833$, giving the common multiplier $0.05/(5\times2.2833)=0.004380$. Its thresholds are therefore approximately $0.004380, 0.008759, 0.01314, 0.01752,$ and $0.02190$. The first three ordered p-values pass, while $0.020>0.01752$ and $0.040>0.02190$. Benjamini–Yekutieli rejects three tests. The results differ because the more permissive independence-based guarantee is not the arbitrary-dependence guarantee.

The step-up decision can also be reported through Benjamini–Yekutieli adjusted p-values. For the p-value of rank $i$, define the adjusted value as the smallest adjusted threshold level that would allow that rank to be rejected, enforcing monotonicity across ranks:

p^{\mathrm{BY}}_{(i)} =\min\left\{1,\ \min_{j\ge i}\frac{m c_m p_{(j)}}{j}\right\}.

At target level $\alpha$, reject the hypotheses whose adjusted p-values are at most $\alpha$. The adjusted values are not posterior probabilities that individual discoveries are false; they are thresholds for decisions under this multiple-testing procedure.

Worked Example: Calculating Benjamini–Yekutieli Adjusted P-Values

Suppose a hypothetical analysis tests $m=4$ variants and obtains ordered, valid p-values $0.001, 0.006, 0.012,$ and $0.040$. Here $c_4=2.0833$, so $m c_m=8.333$. The values $m c_m p_{(j)}/j$, calculated rank by rank, are approximately $0.008333, 0.02500, 0.03333,$ and $0.08333$.

For rank 1, the minimum over ranks $j\ge1$ is $0.008333$. For rank 2, the minimum over ranks $j\ge2$ is $0.02500$; for rank 3 it is $0.03333$, and rank 4 has adjusted value $0.08333$. At $\alpha=0.05$, the first three adjusted p-values are at most 0.05, but the fourth is not. Thus the procedure rejects the first three hypotheses. This agrees with the step-up calculation: the rank-3 threshold is $3(0.05)/(4\times2.0833)=0.01800$, which exceeds $0.012$, while the rank-4 threshold is $0.02400$, below $0.040$.

Conservatism and Practical Interpretation

The adjustment becomes more pronounced as the number of tests grows because $c_m$ increases with $m$. For large $m$, it grows approximately like $\log(m)+0.5772$, where $\log$ is the natural logarithm. Thus, at the same $\alpha$, the Benjamini–Yekutieli thresholds are the Benjamini–Hochberg thresholds divided by $c_m$. This can reduce the number of discoveries, especially in a large family, and may reduce power to detect genuine associations.

That conservatism is the cost of a guarantee that does not require a particular dependence structure. It does not mean that every correlated dataset needs this correction: if a suitable dependence condition is justified, another procedure may be appropriate. But a researcher should not treat correlation alone as proof that a given condition holds. Conversely, the arbitrary-dependence guarantee does not repair invalid input p-values. Poor calibration of single-test p-values can undermine the guarantee before the multiple-testing correction is applied.

Common pitfall: “The variants are correlated, so Benjamini–Hochberg must fail; Benjamini–Yekutieli then guarantees that at most 5% of the rejected variants are false.” The first claim is too broad: the Benjamini–Hochberg guarantee can hold under some positive dependence structures. The second confuses an expected proportion with a fixed bound in one dataset. Benjamini–Yekutieli controls the expected false discovery proportion at no more than the chosen level when true-null p-values are valid; it does not ensure that the realized fraction in every analysis is below that level.

Check Your Understanding

Use the definitions and procedure above to answer the following questions.

In a family of $m=3$ tests, calculate $c_3$ and the Benjamini–Yekutieli threshold at rank 2 for $\alpha=0.05$.
What condition must each true-null p-value satisfy for the arbitrary-dependence guarantee to apply?
Why can linkage disequilibrium make association-test p-values dependent, and does the Benjamini–Yekutieli procedure require independence?
Explain why an FDR level of 0.05 is not a guarantee that exactly 5% or fewer of the rejected associations in one dataset are false.
How does the harmonic factor affect the thresholds compared with Benjamini–Hochberg, and what practical consequence can that have?

What You'll Learn