How to deal with identical p-values with the Benjamini-Hochberg method for correcting for multiple testing

Question

I'm testing for mutations in DNA and I'm using the Benjamini-Hochberg method to modify the threshold I'm testing my p-values against. The method is basically to rank the p-values and compare them to a new threshold defined by

$q(k) = k/n*\alpha$,

where $\alpha$ is your original threshold, commonly 0.05, $k$ is the rank of the p-value you compare with, and $n$ is the total amount of p-values (== total amount of tests).

Due to the nature of my data, many p-values are identical. The following ordered set of p-values could be an example:

$p_1$ = 0.01 $p_2$ = 0.03 $p_3$ = 0.03 $p_4$ = 0.03 $p_5$ = 0.09

Assuming a threshold of $\alpha$ = 0.05, what I have to compare my values to are:

$q_1$ = 0.01 $q_2$ = 0.02 $q_3$ = 0.03 $q_4$ = 0.04 $q_5$ = 0.05

And thus I accept the second test and reject the third and fourth even if they have the same original p-value. I tried to search the literature for a solution to this, but nothing came up. Is there an established method for this? And, if not, what is your preferred ad hoc solution?

Answer 1

Verifying also in the multtest package from Bioconductor, I would suggest to give them the same rank - and very importantly - increment the rank by one for the following p-value(s) rather than using their index+1 in an array! This would have the following result:

considering your examplemulttest's BH would rank $r_1$: 1, $r_2$: 2, $r_3$: 2, $r_4$: 2, $r_5$: 3 rather than of $r_2$: 2, $r_3$: 2, $r_4$: 2, $r_5$: 5

Answer 2

Possible ad hoc solution is to give repeated p-values the same rank.

Answer 3

All answers here are misleading:

considering your p-values are:

p1 = 0.001 p2 = 0.03 p3 = 0.03 p4 = 0.03 p5 = 0.09

Then, BH would rank them: r1 = 1, r2 = 4, r3 = 4, r4 = 4, r5 = 5

Means the ranks for identical p-values will be the max index of them in the sorted list of p-values.

Verified with R using the p.adjust function:

p.adjust(c(0.001, 0.03, 0.03, 0.03, 0.09), method = 'BH')

Yielding the next adjusted p-values:

0.0050 0.0375 0.0375 0.0375 0.0900

Note that for the identical p-values, we got an adjusted p-value of:

0.03 * n / k where n = 5 (as the number of p-values) and k = 4 which is the max index of the identical p-values in the sorted list of p-values...

Accordingly, the adjusted values of alpha are:

alpha * k / n yielding:

q1 = 0.001 q2 = 0.04 q3 = 0.04 q4 = 0.04 q5 = 0.05

Added this part after I was asked:

multtest does the same thing using the same ranks as I mentioned:

using mt.rawp2adjp(c(0.001, 0.03, 0.03, 0.03, 0.09), 'BH') yields the same p-values as the p.adjust function, therefore it ranks identical p-values by the max index of them in the sorted list of p-values. The output of the relevant function from multtest is:

      rawp     BH
[1,] 0.001 0.0050
[2,] 0.030 0.0375
[3,] 0.030 0.0375
[4,] 0.030 0.0375
[5,] 0.090 0.0900

Where the rawp column stands for the original p-values and BH stands for the adjusted p-values. From the adjusted p-values you can calculate the respective alphas as I did above.

Answer 4

These procedures can be confusing. This nice post by spätzle explains the Benjamin Hochberg procedure very well

The Benjamini-Hochberg method is as follows:

1. Order the p-values $p_{(1)},...,p_{(m)}$ and then respectively the hypotheses $H_{0,(1)},...,H_{0,(m)}$
2. Mark as $i_0$ the largest $i$ for which $p_{(i)}\le \frac{i}{m}\alpha$
3. Reject $H_{0,(1)},...,H_{0,(i_0)}$

So if you have several $p_{(i)}$ with the same value, then all that counts is the largest rank. If $p_{(i)}\le \frac{i}{m}\alpha$ for some rank $i$ then also $p_{(j)}\le \frac{j}{m}\alpha$ for some rank $j>i$.

See in the image below for a visual explanation of the method. What counts is the highest ranked p-value that is still below the line. In this case it is the 10-th p-value. All the previous hypothesis tests will be rejected (even if they are above the line). So when you have identical p-values, then what counts is the highest rank.

---

The following ordered set of p-values could be an example:

$p_1$ = 0.01 $p_2$ = 0.03 $p_3$ = 0.03 $p_4$ = 0.03 $p_5$ = 0.09

Assuming a threshold of $\alpha$ = 0.05, what I have to compare my values to are:

$q_1$ = 0.01 $q_2$ = 0.02 $q_3$ = 0.03 $q_4$ = 0.04 $q_5$ = 0.05

And thus I accept the second test and reject the third and fourth even if they have the same original p-value.

You can not have the situation where you accept the second ranked hypothesis while rejecting the third and fourth. You decide on some boundary and all the tests below it are rejected and all the test above it are accepted.

In the BH procedure the second test for which you got $p_2 > q_2$ will be rejected as well.