Tutorials › High-Dimensional Genomic Data Analysis › Benjamini-Yekutieli Control Under Arbitrary Dependence

False Discovery Rate · Tutorial 160 of 795

Benjamini-Yekutieli Control Under Arbitrary Dependence

Learn why dividing the BH thresholds by a harmonic number protects false discovery rate under arbitrary dependence, and what that guarantee does—and does not—require.

Intermediate 11 min read

What You'll Learn

  • State the marginal validity condition required of every true-null p-value.
  • Derive the general arbitrary-dependence bound for a step-up procedure.
  • Obtain the Benjamini–Yekutieli harmonic-number correction and its FDR guarantee.
  • Calculate BY thresholds and interpret the guarantee in a genomic testing family.
  • Recognize why dependence correction cannot repair invalid p-values.
  • Compare BY with BH under independence or PRDS and with Bonferroni familywise control.

Why a Dependence-Robust FDR Procedure Is Useful

Genomic analyses often test thousands of hypotheses whose test statistics are related. Genes may be correlated because of shared pathways, nearby variants may be associated through linkage disequilibrium, and tests may share samples or model components. The Benjamini–Hochberg (BH) procedure controls the false discovery rate under independence and under positive regression dependence on a subset (PRDS). When neither condition is defensible, a procedure with a guarantee under arbitrary dependence is useful.

The Benjamini–Yekutieli (BY) procedure keeps the BH step-up structure but makes its critical values smaller by a harmonic-number factor. This correction can be substantial when the number of tests is large. Its benefit is a guarantee that does not require independence or a particular dependence pattern; its cost is usually fewer rejections.

Setup and Assumptions

Suppose \(m\) hypotheses are tested, with p-values \(p_1,\ldots,p_m\). Let \(m_0\) be the number of true null hypotheses. Order the p-values as \(p_{(1)}\le\cdots\le p_{(m)}\). For a nondecreasing sequence of critical values \(a_1,\ldots,a_m\), a step-up procedure finds \(R=\max\{k:p_{(k)}\le a_k\}\), taking \(R=0\) if there is no such \(k\), and rejects the hypotheses corresponding to \(p_{(1)},\ldots,p_{(R)}\). The false discovery proportion is \(V/\max(R,1)\), where \(V\) is the number of rejected true nulls; the false discovery rate is its expectation, \(E[V/\max(R,1)]\).

Definition: A p-value for a true null is marginally valid (or super-uniform) if \(P(p_i\le t)\le t\) for every \(t\in[0,1]\). The BY guarantee requires this condition for each true-null p-value. It places no independence, positive-dependence, or other joint-distribution condition on the p-values.

Marginal validity is a statement about each null p-value separately; it does not say that p-values are independent. The test models and sampling procedures must support valid p-values for their null hypotheses. Dependence adjustment at the multiple-testing stage cannot compensate for p-values that are systematically too small under true nulls.

The General Step-Up Bound

The key result is a bound for any nondecreasing step-up critical values. It explains why the relevant quantity is the increment in each threshold, rather than the threshold itself.

Theorem: If \(0\le a_1\le\cdots\le a_m\le1\), and every true-null p-value is marginally valid, then under arbitrary dependence, $$ \operatorname{FDR}\le m_0\sum_{k=1}^{m}\frac{a_k-a_{k-1}}{k}, \qquad a_0=0. $$

Here is a justification. For a particular true null \(i\), let \(B_k\) be the event that hypothesis \(i\) is rejected and the total number of rejections \(R\) is at most \(k\). On \(B_k\), the step-up rule gives \(p_i\le a_R\le a_k\). Thus marginal validity implies \(P(B_k)\le P(p_i\le a_k)\le a_k\). Write \(x_k=P(i\text{ is rejected and }R=k)\), so \(P(B_k)=\sum_{j=1}^k x_j\). Summation by parts gives

$$ \begin{aligned} E\left[\frac{\mathbf{1}\{i\text{ is rejected}\}}{R}\right] &=\sum_{k=1}^{m}\frac{x_k}{k}\\ &=\frac{P(B_m)}{m}+ \sum_{k=1}^{m-1}P(B_k)\left(\frac{1}{k}-\frac{1}{k+1}\right)\\ &\le \frac{a_m}{m}+ \sum_{k=1}^{m-1}a_k\left(\frac{1}{k}-\frac{1}{k+1}\right)\\ &=\sum_{k=1}^{m}\frac{a_k-a_{k-1}}{k}. \end{aligned} $$

The term on the left is the expected contribution of null \(i\) to the false discovery proportion; it is defined as zero when there are no rejections. Summing this bound over the \(m_0\) true nulls proves the theorem, because \(V/\max(R,1)\) is the sum of those individual contributions. The argument uses only each null p-value’s marginal validity, which is why arbitrary dependence is permitted.

From the General Bound to BY Thresholds

For a target FDR level \(q\in(0,1)\), define the \(m\)th harmonic number \(c_m=\sum_{j=1}^{m}1/j\). BY uses the critical values \(a_k=kq/(m c_m)\), for \(k=1,\ldots,m\). These are nondecreasing and lie in \([0,1]\), since \(a_m=q/c_m\le1\) for \(q\le1\). Each increment is \(a_k-a_{k-1}=q/(m c_m)\). Substituting into the theorem gives

$$ \operatorname{FDR} \le m_0\sum_{k=1}^{m}\frac{q/(m c_m)}{k} =\frac{m_0}{m}q \le q. $$

The factor \(c_m\) is the dependence correction: compared with BH, every critical value is divided by the same harmonic number. The bound is at most \(q\) even if the p-values have arbitrary dependence, provided the null p-values are valid. It can be smaller than \(q\) when \(m_0<m\), but the procedure does not need to know \(m_0\).

Worked Example: How the Harmonic Factor Changes Thresholds

Suppose a hypothetical genomic analysis tests \(m=5\) hypotheses at target level \(q=0.05\). Calculate the harmonic factor and the BY critical values, then compare them with the BH values.

First, \(c_5=1+1/2+1/3+1/4+1/5=2.283\) to four significant figures. The BY thresholds are \(a_k=k(0.05)/(5(2.283))\). The denominator is approximately \(11.42\), so the thresholds for ranks 1 through 5 are approximately \(0.004380,\ 0.008759,\ 0.01314,\ 0.01752,\ 0.02190\).

For BH, the corresponding thresholds are \(k(0.05)/5\), or \(0.01000,\ 0.02000,\ 0.03000,\ 0.04000,\ 0.05000\). Each BY threshold is the BH threshold divided by \(2.283\). If exactly \(m_0=3\) nulls are true, the general bound gives \(\operatorname{FDR}\le(3/5)(0.05)=0.030\); without knowing \(m_0\), the procedure still guarantees \(\operatorname{FDR}\le0.05\). The thresholds are decision cutoffs for ordered p-values, not confidence limits or probabilities that individual hypotheses are true.

Worked Example: Checking the Increment Bound

Consider \(m=3\) tests with step-up thresholds \(a_1=a_2=a_3=0.30\). One null p-value \(U\) is Uniform\((0,1)\), and the two other p-values are fixed at 1. These are hypothetical values chosen to illustrate the bound; the null p-value is valid, and no independence issue arises because the other p-values are constant.

If \(U\le0.30\), the ordered p-values are \(U,1,1\). Only rank 1 meets its threshold, so the null is rejected and the false discovery proportion is 1. If \(U>0.30\), there are no rejections and the false discovery proportion is 0. Therefore \(\operatorname{FDR}=P(U\le0.30)=0.30\).

The increment bound gives \(m_0[(0.30-0)/1+(0.30-0.30)/2+(0.30-0.30)/3]=0.30\), as required. This example also shows why a bound formed by summing \(a_k/[k(m-k+1)]\) is not valid in general: here it would give \(0.30(1/3+1/4+1/3)=0.275\), below the actual FDR. The correct general step-up bound uses \((a_k-a_{k-1})/k\).

Worked Example: A Small Ordered List Under BY

Suppose a hypothetical analysis tests \(m=4\) hypotheses at \(q=0.10\), with ordered p-values \(0.008,\ 0.020,\ 0.041,\ 0.200\). Assume each true-null p-value is valid, and allow arbitrary dependence. The number of true nulls is not needed to apply the procedure.

Here \(c_4=1+1/2+1/3+1/4=2.083\), so the BY cutoffs are \(k(0.10)/(4(2.083))\): approximately \(0.01200,\ 0.02400,\ 0.03600,\ 0.04800\). The ordered p-values meet the cutoffs at ranks 1 and 2; rank 3 fails because \(0.041>0.03600\), and rank 4 fails because \(0.200>0.04800\). Thus \(R=2\), and the first two hypotheses are rejected.

That conclusion describes this particular list. The FDR guarantee is about the expected false discovery proportion over repeated applications under the assumed data-generating conditions, not a claim that exactly 10% of these two rejections are false. If, for example, \(m_0=2\), the theorem bounds the procedure’s FDR by \((2/4)(0.10)=0.050\); this does not identify which rejected hypotheses are true.

Choosing BY, Diagnosing Problems, and Understanding Limits

ProcedureDependence condition for its standard FDR guaranteePractical implication
BHIndependence or PRDS dependence, with valid null p-valuesUsually less conservative than BY when its dependence condition is justified
BYArbitrary dependence, with valid null p-valuesUses BH thresholds divided by \(c_m\); robust guarantee can cost power
BonferroniValid individual null p-values; no dependence condition for familywise controlTargets familywise error rather than FDR and can be substantially more conservative for many discoveries

Under independence or PRDS, the earlier material on BH control explains why the unmodified BH thresholds can control FDR. BY is a conservative alternative when that condition is not justified and an arbitrary-dependence guarantee is desired. Its penalty grows slowly: \(c_m\) is approximately \(\log(m)+\gamma\) for large \(m\), where \(\gamma\) is Euler’s constant. Thus a large testing family can have BY cutoffs much smaller than BH cutoffs. Bonferroni answers a different error-control question: it controls the probability of one or more false rejections, rather than the expected false discovery proportion.

Useful checks include confirming that the stated family of tests is the family being controlled, that every true-null p-value is calibrated under its null, and that the sorted values and rank-specific cutoffs have been compared consistently. Strong dependence is not itself evidence that BY assumptions fail; arbitrary dependence is allowed. In contrast, model misspecification, unaccounted selection, or invalid null p-values can break marginal validity. A harmonic correction does not repair those problems. Nor does a BY guarantee establish that the chosen tests, biological interpretations, or follow-up claims are scientifically sound.

Common Mistakes

  • Using the BH thresholds while claiming arbitrary-dependence control. The standard BH guarantee needs independence or an appropriate positive-dependence condition. For the general BY guarantee, divide each BH threshold by \(c_m\).
  • Using a false general step-up bound. Under arbitrary dependence, the valid bound stated here involves threshold increments, \(\sum_k(a_k-a_{k-1})/k\), not an unverified rank-weighted sum of the thresholds.
  • Calling \(q\) the probability that a particular rejection is false. The guarantee concerns the expected false discovery proportion across repetitions. It does not give a posterior probability for any one finding.
  • Assuming BY makes invalid p-values safe. The proof requires marginal validity for every true-null p-value; multiplying thresholds by a harmonic factor does not ensure that condition.
  • Expecting BY to produce the same number of discoveries as BH. Since \(c_m\ge1\), BY thresholds are no larger than BH thresholds and can be much smaller in large families.
Key takeaway: For valid true-null p-values and any dependence structure, BY controls FDR at no more than \(q\) by applying the step-up thresholds \(kq/(m c_m)\), where \(c_m=\sum_{j=1}^{m}1/j\). The guarantee is robust to dependence, not to invalid p-values.

Check Your Understanding

Try each question, then check the solutions below.

  1. What assumption about a true-null p-value is needed for the BY guarantee, and what dependence assumption is needed?
  2. For \(m=3\) and \(q=0.06\), calculate \(c_3\) and the BY critical value at rank 2.
  3. For step-up thresholds \(a_1=0.02,\ a_2=0.05,\ a_3=0.09\), calculate the increment sum \(\sum_{k=1}^{3}(a_k-a_{k-1})/k\).
  4. In words, what does an FDR bound of 0.05 mean, and what does it not mean?

Solutions. (1) Each true-null p-value must satisfy \(P(p_i\le t)\le t\) for every \(t\in[0,1]\); no dependence assumption is required. (2) \(c_3=1+1/2+1/3=1.833\), so \(a_2=2(0.06)/(3(1.833))\approx0.02182\). (3) With \(a_0=0\), the sum is \(0.02/1+0.03/2+0.04/3=0.04833\). (4) Across repeated applications under the assumptions, the expected proportion of rejections that are false is at most 0.05. It does not mean that each rejection has a 5% chance of being false or that exactly 5% of a particular discovery list is false.

Clinical Biostats