| p.adjust {stats} | R Documentation |
Adjust P-values for Multiple Comparisons
Description
Given a set of p-values, returns p-values adjusted using one of several methods.
Usage
p.adjust(p, method = p.adjust.methods, n = length(p))
p.adjust.methods
# c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY",
# "fdr", "none")
Arguments
p |
numeric vector of p-values (possibly with
|
method |
correction method, a |
n |
number of comparisons, must be at least |
Details
The adjustment methods include the Bonferroni correction
("bonferroni") in which the p-values are multiplied by the
number of comparisons. Less conservative corrections are also
included by
Holm (1979) ("holm"),
Hochberg (1988) ("hochberg"),
Hommel (1988) ("hommel"),
Benjamini and Hochberg (1995) ("BH" or its alias
"fdr"), and
Benjamini and Yekutieli (2001) ("BY"), respectively.
A pass-through option ("none") is also included.
The set of methods are contained in the p.adjust.methods vector
for the benefit of methods that need to have the method as an option
and pass it on to p.adjust.
The first four methods are designed to give strong control of the family-wise error rate. There seems no reason to use the unmodified Bonferroni correction because it is dominated by Holm's method, which is also valid under arbitrary assumptions.
Hochberg's and Hommel's methods are valid when the hypothesis tests are independent or when they are non-negatively associated (Sarkar 1998;Sarkar and Chang 1997). Hommel's method is more powerful than Hochberg's, but the difference is usually small and the Hochberg p-values are faster to compute.
The "BH" (aka "fdr") and "BY" methods of
Benjamini, Hochberg, and Yekutieli control the false discovery rate,
FDR, the expected proportion of false discoveries amongst the
rejected hypotheses.
The "BY" correction modifies the BH procedure by replacing the
target level q with q / \sum_{i=1}^{m} 1/i, where m is
the number of tests (Theorem 1.3 in the reference), which controls the
FDR under the most general form of dependence structure. This will
be more conservative than "BH", for small p even more than
Bonferroni, see the example.
The FDR as implemented by the "BH" method is a less stringent
condition than the family-wise error rate, so it is typically more
powerful than the others.
Note that you can set n larger than length(p) which
means the unobserved p-values are assumed to be greater than all the
observed p for "bonferroni" and "holm" methods and equal
to 1 for the other methods.
Value
A numeric vector of corrected p-values (of the same length as
p, with names copied from p).
References
Benjamini Y., Hochberg Y. (1995). “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.” Journal of the Royal Statistical Society Series B: Statistical Methodology, 57(1), 289–300. doi:10.1111/j.2517-6161.1995.tb02031.x.
Benjamini Y., Yekutieli D. (2001). “The Control of the False Discovery Rate in Multiple Testing under Dependency.” The Annals of Statistics, 29(4). doi:10.1214/aos/1013699998.
Hochberg Y. (1988). “A Sharper Bonferroni Procedure for Multiple Tests of Significance.” Biometrika, 75(4), 800–802. doi:10.1093/biomet/75.4.800.
Holm S. (1979). “A Simple Sequentially Rejective Multiple Test Procedure.” Scandinavian Journal of Statistics, 6, 65–70. https://www.jstor.org/stable/4615733.
Hommel G. (1988). “A Stagewise Rejective Multiple Test Procedure Based on a Modified Bonferroni Test.” Biometrika, 75(2), 383–386. doi:10.1093/biomet/75.2.383.
Sarkar S. K. (1998).
“Some Probability Inequalities for Ordered \mathrm{MTP}_2 Random Variables: A Proof of the Simes Conjecture.”
The Annals of Statistics, 26(2).
doi:10.1214/aos/1028144846.
Sarkar S. K., Chang C. (1997). “The Simes Method for Multiple Hypothesis Testing with Positively Dependent Test Statistics.” Journal of the American Statistical Association, 92(440), 1601–1608. doi:10.1080/01621459.1997.10473682.
Shaffer J P (1995).
“Multiple Hypothesis Testing.”
Annual Review of Psychology, 46(1), 561–584.
doi:10.1146/annurev.ps.46.020195.003021.
(An excellent review of the area.)
Wright S. P. (1992).
“Adjusted P-Values for Simultaneous Inference.”
Biometrics, 48(4), 1005.
doi:10.2307/2532694.
(Explains the adjusted P-value approach.)
See Also
pairwise.* functions such as pairwise.t.test.
Examples
require(graphics)
set.seed(123)
x <- rnorm(50, mean = c(rep(0, 25), rep(3, 25)))
p <- 2*pnorm(sort(-abs(x)))
round(p, 3)
round(p.adjust(p), 3)
round(p.adjust(p, "BH"), 3)
## or all of them at once (dropping the "fdr" alias):
p.adjust.M <- p.adjust.methods[p.adjust.methods != "fdr"]
p.adj <- sapply(p.adjust.M, function(meth) p.adjust(p, meth))
p.adj.60 <- sapply(p.adjust.M, function(meth) p.adjust(p, meth, n = 60))
stopifnot(identical(p.adj[,"none"], p), p.adj <= p.adj.60)
round(p.adj, 3)
## or a bit nicer:
head(round(100 * p.adj[,c(7,1:6)], 2), n=21) # in [percent]:
## none holm hochberg hommel bonferroni BH BY
## [1,] 0.00 0.00 0.00 0.00 0.00 0.00 0.01 *)
## [2,] 0.00 0.10 0.10 0.10 0.11 0.04 0.19 *)
## [3,] 0.00 0.12 0.12 0.12 0.13 0.04 0.19 *)
## [4,] 0.01 0.46 0.46 0.42 0.49 0.09 0.43
## [5,] 0.01 0.48 0.48 0.45 0.53 0.09 0.43
## .... .......... ............ ...............
## .... .......... ............ ................
## [18,] 0.88 29.06 29.06 27.30 44.03 2.45 11.00
## [19,] 0.94 30.08 30.08 29.14 47.01 2.47 11.13
## [20,] 1.13 35.02 35.02 33.89 56.49 2.82 12.71
## [21,] 2.12 63.45 63.45 57.11 100.00 5.04 22.66
##
## *) The smallest 3 Bonferroni values are smaller than the "BY" ones,
## (John Maindonald, PR#17136)
## number of rejected H0 ("P" < 0.05):
colSums(p.adj < 0.05)
## holm hochberg hommel bonferroni BH BY none
## 11 11 11 11 20 12 22
## visual comparison
matplot(p, p.adj, ylab="p.adjust(p, meth)", type = "l", asp = 1, lty = 1:6,
col = 1:7, main = "P-value adjustments")
legR <- function() {
legend("bottomright", p.adjust.M, col = 1:7, lty = 1:6, bty = "n", inset = 0.05)
rug(p) }
legR()
## zoom in & log scale
lim <- c(7e-7, .20)
matplot(p, p.adj, ylab="p.adjust(p, meth)", type = "l", asp = 1, lty = 1:6, col = 1:7,
main = "P-value adjustments [log-log]", log = "xy", xlim=lim, ylim=lim, las=1)
legR()
## Can work with NAs:
pN <- p; iN <- c(46, 47); pN[iN] <- NA
pN.a <- sapply(p.adjust.M, function(meth) p.adjust(pN, meth))
## The smallest 20 P-values all affected by the NAs :
round((pN.a / p.adj)[1:20, ] , 4)