[R] Benjamini-Hochberg (BH) Q value (false discovery rate)

[David L Carlson]

R is open source so you can look at the code. 
In this case just typing the function name gets it:

> p.adjust
function (p, method = p.adjust.methods, n = length(p)) 
{
    method <- match.arg(method)
    if (method == "fdr") 
        method <- "BH"
    . . . . . Stuff deleted . . . .  
# The code for Benjamini and Hochberg:
    }, BH = {
        i <- lp:1L
        o <- order(p, decreasing = TRUE)
        ro <- order(o)
        pmin(1, cummin(n/i * p[o]))[ro]

# A quick example
p <- c(0.7723, 0.1992, 0.1821, 0.6947, 0.0932, 0.1484, 0.0006, 0.0004)
round(p.adjust(p, method="BH"), 4)
#  [1] 0.7723  0.2656  0.2656  0.7723  0.2485  0.2656  0.0024  0.0024

n <- length(p)
i <- n:1L
o <- order(p, decreasing=TRUE)    # Order from largest to smallest
ro <- order(o)     # Reverse the ordering to return the adjusted values
pmin(1, cummin(n/i * p[o]))[ro]   # Compute the adjusted value

The adjustment orders the p-values from largest to smallest and multiplies them by n/i, in this case

> n/i
[1] 1.000000 1.142857 1.333333 1.600000 2.000000 2.666667 4.000000 8.000000

So the largest is multiplied by 1 and the smallest by 8.

Then cummin() takes the cumulative minimum so that the multiplication does not change the decreasing order of the values. In this example the last value is the same as the second to last instead of 8 times the original value (.0032). The pmin() function ensures that the adjusted value never exceeds 1.

----------------------------------------
David L Carlson
Department of Anthropology
Texas A&M University
College Station, TX 77843-4352


-----Original Message-----
From: R-help <r-help-bounces using r-project.org> On Behalf Of David Bars
Sent: Wednesday, March 6, 2019 3:16 AM
To: r-help using r-project.org
Subject: [R] Benjamini-Hochberg (BH) Q value (false discovery rate)

Dear everybody,

I'm using stats package (version 3.5.2) and in detail their p.adjust()
function in order to control the false discovery rate in multiple
comparisons. In particular I used the Benjamini-Hochberg method.

Nevertheless,  the help of the stats package indicates that the adjusting
method is based on:
Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery
rate: a practical and powerful approach to multiple testing. Journal of the
Royal Statistical Society Series B, 57, 289–300.
http://www.jstor.org/stable/2346101.

But, in the function mentioned above (p.adjust()) there are no possibiity
to define the Q value (the false discovery rate). I understand that the Q
value is calculated automatically but through which method???

For example, in another package (sgof v.2.3 also deposited on CRAN)
indicates that the the false discovery rate is estimated by the method
developed by:

Dalmasso C, Broet P and Moreau T (2005). A simple procedure for estimating
the false discovery rate. *Bioinformatics* 21:660--668.


Many thanks for your help,

David Bars.

	[[alternative HTML version deleted]]

______________________________________________
R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.