7
$\begingroup$

Is there a possibility to calculate or estimate the overall rejection threshold of the Benjamini–Hochberg procedure (BH)?

For the correction of the FWER using the Bonferroni method, the significance threshold is adjusted to the number of evaluated hypotheses $m$ as follows $\bar{\alpha}= \frac{\alpha}{m}$. But since the BH-procedure produces an individual $q$-value for each independent hypothesis that is compared to an apriori defined FDR, I am not sure how this can be done.

$\endgroup$
6
$\begingroup$

As you sense, there is no fixed p-value cutoff for the Benjamini-Hochberg control of false discovery rate. The cutoff depends on the specific distribution of p-values among the $m$ hypotheses that you are evaluating together. You put them in increasing order and count up in $k$ from the lowest p-value $(k=1)$. You agree to "reject the null hypothesis" for hypotheses up through this value of $k$:

For a given $\alpha$, find the largest $k$ such that $P_{(k)} \leq \frac{k}{m} \alpha.$

If the null hypotheses all hold so there is a uniform distribution of p-values in [0,1], the p-value cutoff will be close to $\alpha$. How much below that you go if some null hypotheses don't hold depends on how non-uniform the distribution of p-values is.

$\endgroup$
6
  • $\begingroup$ Did you mean "the q-value cutoff will be close to α"? $\endgroup$ – Sextus Empiricus Nov 5 '20 at 15:01
  • $\begingroup$ @SextusEmpiricus the terminology gets so confusing on this. I intended that to be the cutoff among the hypothesis p-values to achieve the desired q-value for FDR. The q-value is called $\alpha$ in the quote from Wikipedia on the method. $\endgroup$ – EdM Nov 5 '20 at 15:17
  • $\begingroup$ but this cutoff will be $\frac{k}{m} \alpha$, right? $\endgroup$ – Sextus Empiricus Nov 5 '20 at 15:19
  • $\begingroup$ @SextusEmpiricus yes. So does one call the original hypothesis p-values "q-values" when you put them in order to do the B-H FDR control? Or is the "q-value" the FDR? Or is the "q-value" $\frac{k}{m} \alpha$? I don't know if there is a universally accepted terminology for this, would appreciate a reference if there is. $\endgroup$ – EdM Nov 5 '20 at 15:26
  • $\begingroup$ I am actually not sure what the OP means by overall rejection threshold. I don't believe it applies to this situation because there is no single threshold according to which each hypothesis is treated independently. This confuses me about your statement "the p-value cutoff will be...". $\endgroup$ – Sextus Empiricus Nov 5 '20 at 17:07
4
$\begingroup$

I'm not certain of the formal validity of this approach, but you could calculate the corresponding FWER as given by the Hochberg method.

The Benjamini-Hochberg procedure for controlling the False Discovery Rate is (I'm going to quote Wikipedia)

... we have $H_1 \ldots H_m$ null hypotheses tested and $P_1 \ldots P_m$ their corresponding p-values. We list these p-values in ascending order and denote them by $P_{(1)} \ldots P_{(m)}$. ...

  1. For a given $\alpha$, find the largest $k$ such that $P_{(k)} \leq \frac{k}{m} \alpha.$
  2. Reject the null hypothesis (i.e., declare discoveries) for all $H_{(i)}$ for $i = 1, \ldots, k$.

That method sets the FDR at $\alpha$, i.e., out of the rejected hypotheses, we expect the fraction of Type I Errors to be $\alpha$.

The Family-Wise Error Rate, on the other hand, is the probability of at least one Type I Error in the set of rejected hypotheses. The Hochberg method accomplishes this with a calculation similar to the BH FDR method (again, quoting Wikipedia),

  • Start by ordering the p-values (from lowest to highest) $P_{(1)} \ldots P_{(m)}$ and let the associated hypotheses be $H_{(1)} \ldots H_{(m)}$
  • For a given $\alpha$, let $R$ be the largest $k$ such that $P_{(k)} \leq \frac{\alpha}{m-k+1}$
  • Reject the null hypotheses $H_{(1)} \ldots H_{(R)}$

You could put these together to 1) define an FDR $\alpha$, 2) determine the largest rejected p-value $p^*$ and the number of rejected hypotheses $k$ out of total $m$, 3) calculate the corresponding Hochberg FWER $\tilde{\alpha}$ as $$ \tilde{\alpha} = p^* \times(m-k+1) $$

$\endgroup$
1
  • $\begingroup$ I am not sure about Wikipedia or myself, but when I try to simulate and compute the FDR and FWER of the Benjamini-Hochberg procedure then I get that it is the FWER that is equal to $\alpha$. Do you know anything about that stats.stackexchange.com/questions/495211/… ? $\endgroup$ – Sextus Empiricus Nov 5 '20 at 17:00
0
$\begingroup$

After some thought, I believe that the unadjusted p-value of the last (rank-wise) significant test after the BH-procedure comes closest to a significance threshold.

An example:

Do the BH-procedure:

  1. Some p-values: $0.0001,0.0234,0.3354,0.0021,0.5211,0.9123,0.0008,0.0293,0.0500, 1.000$

  2. Order them: $0.0001, 0.0008, 0.0021, 0.0234, 0.0293, 0.0500, 0.3354, 0.5211, 0.9123, 1.0000$

  3. Compute q-values for all 10 ranks: $q_i = \frac{i}{m}\cdot \alpha$, for $i=1,2,..,m$.

  4. Find the largest ranked p-value that is smaller than its corresponding q-value.

Results: $$\begin{array}{} \textbf{Rank} & \textbf{q-value} & \textbf{p-value} & \textbf{Significance (BH)} \\ \hline 1 & 0.005 & 0.0001 & True \\ \hline 2 & 0.01 & 0.0008 & True \\ \hline 3 & 0.015 & 0.0021 & True \\ \hline 4 & 0.02 & 0.0234 & False \\ \hline 5 & 0.025 & 0.0293 & False \\ \hline 6 & 0.03 & 0.05 & False \\ \hline 7 & 0.035 & 0.3354 & False \\ \hline 8 & 0.04 & 0.5211 & False \\ \hline 9 & 0.045 & 0.9123 & False \\ \hline 10 & 0.05 & 1 & False \\ \hline \end{array}$$

On the table, we can see that all tests above Rank 3 are non-significant, thus we can conclude that 0.0021 acts as our significance threshold. In comparison, the Bonferroni correction has a threshold of $\frac{\alpha}{m}=0.005$.

Here is the R-code I used for this example:

# generate p-values
pValues <- c(0.0001,0.0234,0.3354,0.0021,0.5211,0.9123,0.0008,0.0293,0.0500, 1)

# order the p-values
pValues <- sort(pValues)

# BH-procedure
alpha <- 0.05
m <- length(pValues)
qValues <- c()


for (i in 1:m){
  qV <- (i/m)*alpha
  qValues <- append(qValues, qV)
}

# find the largest p-value that satisfies p_i < q_i  
BH_test <- qValues > pValues

# largest k is 3, thus threshold is 0.0021
threshold <- p[sum(BH_test)];threshold
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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