Rejection threshold of the Benjamini-Hochberg procedure

Question

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.

Answer 1

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.

Answer 2

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) $$

Answer 3

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