1#  File src/library/stats/R/p.adjust.R
2#  Part of the R package, https://www.R-project.org
3#
4#  Copyright (C) 1995-2020 The R Core Team
5#
6#  This program is free software; you can redistribute it and/or modify
7#  it under the terms of the GNU General Public License as published by
8#  the Free Software Foundation; either version 2 of the License, or
9#  (at your option) any later version.
10#
11#  This program is distributed in the hope that it will be useful,
12#  but WITHOUT ANY WARRANTY; without even the implied warranty of
13#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14#  GNU General Public License for more details.
15#
16#  A copy of the GNU General Public License is available at
17#  https://www.R-project.org/Licenses/
18
19p.adjust.methods <-
20    c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")
21
22p.adjust <- function(p, method = p.adjust.methods, n = length(p))
23{
24    ## Methods 'Hommel', 'BH', 'BY' and speed improvements
25    ## contributed by Gordon Smyth
26    method <- match.arg(method)
27    if(method == "fdr") method <- "BH"	# back compatibility
28    nm <- names(p)
29    p <- as.numeric(p)
30    p0 <- setNames(p, nm)
31    if(all(nna <- !is.na(p))) nna <- TRUE
32    else p <- p[nna]
33    lp <- length(p)
34    stopifnot(n >= lp)
35    if (n <= 1) return(p0)
36    if (n == 2 && method == "hommel") method <- "hochberg"
37
38    p0[nna] <-
39        switch(method,
40               bonferroni = pmin(1, n * p),
41               holm = {
42                   i <- seq_len(lp)
43                   o <- order(p)
44                   ro <- order(o)
45                   pmin(1, cummax( (n+1L - i) * p[o] ))[ro]
46               },
47               hommel = { ## needs n-1 >= 2 in for() below
48                   if(n > lp) p <- c(p, rep.int(1, n-lp))
49                   i <- seq_len(n)
50                   o <- order(p)
51                   p <- p[o]
52                   ro <- order(o)
53                   q <- pa <- rep.int( min(n*p/i), n)
54                   for (j in (n-1L):2L) {
55                       ij <- seq_len(n-j+1L)
56                       i2 <- (n-j+2L):n
57                       q1 <- min(j*p[i2]/(2L:j))
58                       q[ij] <- pmin(j*p[ij], q1)
59                       q[i2] <- q[n-j+1L]
60                       pa <- pmax(pa, q)
61                   }
62                   pmax(pa, p)[if(lp < n) ro[1L:lp] else ro]
63               },
64               hochberg = {
65                   i <- lp:1L
66                   o <- order(p, decreasing = TRUE)
67                   ro <- order(o)
68                   pmin(1, cummin( (n+1L - i) * p[o] ))[ro]
69               },
70               BH = {
71                   i <- lp:1L
72                   o <- order(p, decreasing = TRUE)
73                   ro <- order(o)
74                   pmin(1, cummin( n / i * p[o] ))[ro]
75               },
76               BY = {
77                   i <- lp:1L
78                   o <- order(p, decreasing = TRUE)
79                   ro <- order(o)
80                   q <- sum(1/(1L:n))
81                   pmin(1, cummin(q * n / i * p[o]))[ro]
82               },
83               none = p)
84    p0
85}
86