A Worked Benjamini-Yekutieli Calculation

Why Work Through a Benjamini–Yekutieli Calculation?

In genomic analysis, a single experiment may test thousands of hypotheses. The Benjamini–Yekutieli (BY) procedure is a step-up method for selecting discoveries while controlling the false discovery rate under arbitrary dependence among tests, provided the null p-values are valid. Its adjustment is more conservative than the Benjamini–Hochberg (BH) adjustment: it divides the BH rank thresholds by a harmonic factor. That factor can noticeably change the number of reported discoveries, so it is important to calculate it and apply the step-up rule carefully.

The earlier material on Benjamini–Yekutieli control under arbitrary dependence establishes the procedure’s control result. Here the emphasis is its numerical use: ordering a p-value list, calculating the thresholds, locating the largest passing rank, and identifying exactly which hypotheses are rejected.

Definition: For \(m\) tested hypotheses, let \(p_{(1)}\le p_{(2)}\le\cdots\le p_{(m)}\) denote the p-values in increasing order, and define the harmonic correction \(c_m=\sum_{j=1}^{m}1/j\). At target level \(\alpha\), the BY threshold for rank \(i\) is \(i\alpha/(m c_m)\). Let \(k\) be the largest rank satisfying \(p_{(k)}\le k\alpha/(m c_m)\). If such a rank exists, reject the hypotheses corresponding to \(p_{(1)},\ldots,p_{(k)}\); if none exists, reject none.

Assumptions and the Step-Up Rule

The procedure is applied to a defined family of \(m\) hypotheses, with one p-value for each test. A null p-value must be valid: under its null hypothesis, \(P(p\le t)\le t\) for every \(t\) between 0 and 1. The p-values may be dependent in any way; independence is not required for the BY guarantee. With these conditions, the procedure controls the false discovery rate at level \(\alpha\):

\operatorname{FDR}=E\left[\frac{V}{\max(R,1)}\right]\le\alpha,

where \(R\) is the number of rejected hypotheses and \(V\) is the number of rejected true null hypotheses. This is a bound on an expected error rate across repeated applications, not a statement that a particular rejection set contains a known fraction of false discoveries.

The calculation has two parts. First, compute the same harmonic correction for every rank. Second, compare each ordered p-value with its rank-specific threshold. The rule is step-up because the largest rank that passes determines the rejection set: if rank \(k\) passes, all ranks below \(k\) are rejected, even if an isolated lower-ranked comparison would not pass its own threshold. The p-values must be sorted before applying the comparisons.

Worked Calculations

Worked Example: Five Tests at a 5% Target Level

Suppose five hypothetical genomic tests have p-values \(0.001, 0.006, 0.012, 0.019,\) and \(0.040\). Assume the family is these five tests, the null p-values are valid, and the desired level is \(\alpha=0.05\). The values are already in increasing order.

1. Calculate the harmonic correction. For \(m=5\),

c_5=1+\frac12+\frac13+\frac14+\frac15 =\frac{137}{60}=2.283333\ldots

Thus \(m c_m=5(2.283333\ldots)=11.416666\ldots\approx11.4167\). Keeping the unrounded harmonic sum through the multiplication avoids a rounding inconsistency.

2. Calculate the thresholds. The threshold at rank \(i\) is \(i(0.05)/(11.416666\ldots)\). The resulting thresholds, rounded to four significant figures, are:

\begin{array}{c|ccccc} i & 1 & 2 & 3 & 4 & 5\\ \hline i\alpha/(m c_m) & 0.004380 & 0.008759 & 0.01314 & 0.01752 & 0.02190 \end{array}

3. Find the largest passing rank. At rank 1, \(0.001\le0.004380\), so it passes. At rank 2, \(0.006\le0.008759\), so it passes. At rank 3, \(0.012\le0.01314\), so it passes. At rank 4, \(0.019>0.01752\), so it fails; rank 5 also fails because \(0.040>0.02190\). The largest passing rank is \(k=3\).

Conclusion. Reject the hypotheses associated with the three smallest p-values, \(0.001, 0.006,\) and \(0.012\). Under the stated validity conditions, the BY procedure controls the FDR at no more than 0.05 across repeated applications; it does not guarantee that exactly 5% of these three rejections are false.

Worked Example: Four Tests and a Different Target Level

Consider the hypothetical ordered p-values \(0.004, 0.020, 0.030,\) and \(0.080\), with \(m=4\) tests and \(\alpha=0.10\). Suppose again that each null p-value is valid.

The harmonic correction is

c_4=1+\frac12+\frac13+\frac14 =\frac{25}{12}=2.083333\ldots, \qquad m c_m=4(2.083333\ldots)=8.333333\ldots\approx8.3333.

The rank thresholds are \(i(0.10)/(25/3)\), or \(0.012, 0.024, 0.036,\) and \(0.048\). Each comparison is: \(0.004\le0.012\), \(0.020\le0.024\), \(0.030\le0.036\), and \(0.080>0.048\). The largest passing rank is \(k=3\), so reject the hypotheses corresponding to the three smallest p-values. Notice that the target level has changed the thresholds, but the harmonic correction still uses the number of tests in the family.

Worked Example: Comparing BY and BH on the Same List

For six hypothetical tests, let the ordered p-values be \(0.001, 0.004, 0.009, 0.015, 0.030,\) and \(0.070\), with \(\alpha=0.05\). Assume valid null p-values. This example compares the two procedures on the same family; it does not change the BY calculation.

For BY, \(c_6=1+1/2+1/3+1/4+1/5+1/6=2.45\), so \(m c_m=6(2.45)=14.7\). The thresholds at ranks 1 through 6 are approximately \(0.003401, 0.006803, 0.01020, 0.01361, 0.01701,\) and \(0.02041\). The first three p-values pass their respective thresholds: \(0.001\le0.003401\), \(0.004\le0.006803\), and \(0.009\le0.01020\). Rank 4 fails because \(0.015>0.01361\), and ranks 5 and 6 fail as well. BY therefore rejects the first three hypotheses.

For comparison, the BH thresholds are \(i(0.05)/6\): approximately \(0.008333, 0.01667, 0.02500, 0.03333, 0.04167,\) and \(0.05000\). The largest passing rank is 5, since \(0.030\le0.04167\), while rank 6 fails because \(0.070>0.05000\). BH rejects the first five hypotheses in this example. This difference illustrates the effect of BY’s harmonic correction; it is not evidence that the two methods have different target levels.

Why the Harmonic Correction Appears

The BY thresholds are the BH thresholds divided by \(c_m\). Since \(c_m\ge1\), each BY threshold is no larger than the corresponding BH threshold. The factor is designed to preserve FDR control even when test statistics, and hence p-values, have arbitrary dependence. As the number of tests increases, \(c_m\) increases, though slowly; for large \(m\), it is approximately \(\log(m)+\gamma\), where \(\gamma\) is Euler’s constant.

The rank adjustment itself follows directly from the step-up definition. For a fixed rank \(i\), the BH threshold is \(i\alpha/m\). Replacing the target level by \(\alpha/c_m\) gives

\frac{i(\alpha/c_m)}{m}=\frac{i\alpha}{m c_m}.

Therefore the BY computation can be carried out using either the explicit denominator \(m c_m\) or the adjusted level \(\alpha/c_m\), provided the same choice is used consistently at every rank. The largest passing rank is then selected, and all hypotheses at or below that rank in the ordered list are rejected.

Checks, Alternatives, and Limitations

Before reporting a rejection set, verify that the p-values are sorted, \(m\) is the number of tests in the intended family, and the harmonic sum includes every integer from 1 through \(m\). Keep adequate precision in intermediate calculations; rounding the correction early can change a comparison when a p-value lies near its threshold. A useful final check is to confirm that the selected rank passes and that every higher rank fails.

BY is a conservative alternative to BH when arbitrary dependence is a concern. BH uses thresholds \(i\alpha/m\) and can be more powerful under independence or suitable positive dependence conditions. BY reduces each threshold by \(c_m\) to obtain its arbitrary-dependence guarantee. This protection can substantially reduce discoveries, especially for a large family. Choosing a procedure therefore requires attention to the dependence assumptions and the scientific aim, not just which method produces more rejections.

The guarantee does not repair invalid p-values: if null p-values are not superuniform, BY’s stated control need not hold. Nor does the procedure establish that rejected hypotheses are biologically important, identify causal effects, or provide a posterior probability that any one null is true. The family definition also matters: adding or omitting tests changes \(m\), the correction, and potentially the rejection set.

Common Mistakes

Using the BH threshold by mistake. Comparing \(p_{(i)}\) with \(i\alpha/m\) omits the BY correction. Use \(i\alpha/(m c_m)\) for BY.
Choosing the first passing rank instead of the largest. The step-up rule selects the largest rank that passes, then rejects all hypotheses ranked up to it.
Changing the family size after sorting. The value of \(m\) is the number of hypotheses in the specified family, not the number of small p-values or the number that pass.
Mixing rounded and unrounded arithmetic. For example, \(5(2.2833)=11.4165\), whereas \(5(2.283333\ldots)=11.416666\ldots\approx11.4167\). State which precision is being used.
Calling the target level a realized false-discovery proportion. The guarantee concerns the expected false discovery proportion across repetitions, not the known fraction of false rejections in one result.

Key takeaway: Compute \(c_m=\sum_{j=1}^{m}1/j\), compare each ordered p-value with \(i\alpha/(m c_m)\), and reject through the largest passing rank. BY’s harmonic correction permits FDR control under arbitrary dependence when the null p-values are valid.

Check Your Understanding

For each problem, assume the p-values are valid under their null hypotheses. Solutions are included after the questions.

For \(m=3\), calculate \(c_m\) and \(m c_m\).
At \(m=3\) and \(\alpha=0.05\), the ordered p-values are \(0.002, 0.010,\) and \(0.020\). Find the BY rejection set.
For \(m=4\), \(c_4=25/12\). At \(\alpha=0.05\), calculate the rank-2 threshold, retaining enough precision to make a comparison.
Explain why a p-value list cannot be analyzed with BY by simply comparing every value with \(\alpha/m\).

Solutions. (1) \(c_3=1+1/2+1/3=11/6=1.833333\ldots\), so \(m c_m=3(11/6)=11/2=5.5\). (2) The thresholds are \(0.05/5.5=0.009091\), \(0.10/5.5=0.01818\), and \(0.15/5.5=0.02727\). The first two p-values pass and the third also passes, since \(0.020\le0.02727\); the largest passing rank is 3, so reject all three. (3) \(m c_m=4(25/12)=25/3\), and the rank-2 threshold is \(2(0.05)/(25/3)=0.012\). (4) The BY rule uses rank-specific thresholds \(i\alpha/(m c_m)\), not a single threshold \(\alpha/m\); it requires sorting and identifying the largest passing rank.