10
$\begingroup$

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?

$\endgroup$

4 Answers 4

6
$\begingroup$

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

$\endgroup$
6
  • $\begingroup$ Thank you. I accept this answer although it is essentially the same as the one by Maciej because of the extra detail. $\endgroup$ Commented Oct 5, 2012 at 12:48
  • $\begingroup$ Thank you - I effectively picked up on the comment of Maciej. I often miss the possibility to accept more than one answer! $\endgroup$
    – dmeu
    Commented Oct 16, 2012 at 13:09
  • $\begingroup$ I doubt that it is correct. In practice very rarely one can get two equal p-values. $\endgroup$
    – Viktor
    Commented Nov 8, 2017 at 15:17
  • $\begingroup$ @dmeu this answer is incomplete and unclear (and incorrect). How did you use multtest to find out these ranks? $\endgroup$ Commented Oct 10, 2022 at 8:44
  • $\begingroup$ @SextusEmpiricus I wish I'd remember, This was ten years ago! If the answer is objectively wrong, could you provide a correction? Do I need to unlock community wiki or similar or can you edit the answer? $\endgroup$
    – dmeu
    Commented Oct 20, 2022 at 8:59
2
$\begingroup$

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

$\endgroup$
1
  • $\begingroup$ I've just exploring my multtest package results and I see that I have a subset of the same FDR corrected values [1] 0.08199 0.08199 0.08199 0.08199 0.08199 0.08199 0.08199 0.08199 0.08199 [10] 0.08199 So maybe it is as I wrote. $\endgroup$ Commented Nov 24, 2011 at 19:36
2
$\begingroup$

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.

$\endgroup$
2
  • 1
    $\begingroup$ Hi, I added an additional part to my answer related to multtest. I hope I managed to understand your question properly. $\endgroup$ Commented Oct 10, 2022 at 7:57
  • 1
    $\begingroup$ I believe your criticism is right. I looked it up as well and it seems indeed that the highest rank counts. I wonder what went different in the other answers. $\endgroup$ Commented Oct 10, 2022 at 8:10
2
$\begingroup$

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.

A visual example of the BH method


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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.