David Bars
2019-Mar-06 09:16 UTC
[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]]
David L Carlson
2019-Mar-06 15:36 UTC
[R] Benjamini-Hochberg (BH) Q value (false discovery rate)
R is open source so you can look at the code. In this case just typing the function name gets it:> p.adjustfunction (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 at r-project.org> On Behalf Of David Bars Sent: Wednesday, March 6, 2019 3:16 AM To: r-help at 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 at 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.