Benjamini–Hochberg procedure in R

[raphael123]

Tell me if I am wrong but this is the Benjamini–Hochberg procedure is that simple ?

# produce a table of FDR using with pval already sorted in increasing order
fdr = function(pval){

 size=length(pval)

 if(size<2) return(pval)

 FDR=c(pval[size])

 for( i in 1:(size-1)) FDR=c(FDR,min(pval[size-i]*(size)/(size-i),FDR[i]))

 return(rev(FDR))
}

[dpryan]

Essentially, though the exact algorithm isn't the same. This is what p.adjust actually does when computing using the BH method. While this is faster (no for loop in R, cummin is presumably written in C where loops aren't a problem), the output is more or less the same (within rounding error at least):

n <- length(p)
lp <- length(p)
i <- lp:1L
o <- order(p, decreasing = TRUE)
ro <- order(o) #This is the index of the initial p-values, highest to lowest
pmin(1, cummin(n/i * p[o]))[ro]

[raphael123]

Thank you ! We have the same output ! Your code is more elegant but harder to understand
equivalent of :

p.adjust(pval,method="BH")

[adiallo]

Hello,
The method:
n <- length(p)
lp <- length(p)
i <- lp:1L
o <- order(p, decreasing = TRUE)
ro <- order(o) #This is the index of the initial p-values, highest to lowest
pmin(1, cummin(n/i * p[o]))[ro]

is given different results than p.adjust.

I did some tests with (300000 pvalue) and nothing was siginificant with p.adjust but 77 was with the upper method.

Cheers,

[dpryan]

If you don't otherwise specify a method, p.adjust() will return the Holm's method adjusted p-value rather than the BH method that I outlined. At least that's what I recall being the case. You should get the same results if you specify "p.adjust(some_p_values, method='BH')".