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