P-value correction: Difference between revisions

m
mNo edit summary
Line 2,817:
 
=={{header|R}}==
The '''p.adjust''' function is built-in, see [https://stat.ethz.ch/R-manual/R-devel/library/stats/html/p.adjust.html R manual].
<lang R>
#p.adjust's source code is below
#> p.adjust
#function (p, method = p.adjust.methods, n = length(p))
#{
# method <- match.arg(method)
# if (method == "fdr")
# method <- "BH"
# nm <- names(p)
# p <- as.numeric(p)
# p0 <- setNames(p, nm)
# if (all(nna <- !is.na(p)))
# nna <- TRUE
# p <- p[nna]
# lp <- length(p)
# stopifnot(n >= lp)
# if (n <= 1)
# return(p0)
# if (n == 2 && method == "hommel")
# method <- "hochberg"
# p0[nna] <- switch(method, bonferroni = pmin(1, n * p), holm = {
# i <- seq_len(lp)
# o <- order(p)
# ro <- order(o)
# pmin(1, cummax((n - i + 1L) * p[o]))[ro]
# }, hommel = {
# if (n > lp) p <- c(p, rep.int(1, n - lp))
# i <- seq_len(n)
# o <- order(p)
# p <- p[o]
# ro <- order(o)
# q <- pa <- rep.int(min(n * p/i), n)
## for (j in (n - 1):2) {
# ij <- seq_len(n - j + 1)
# i2 <- (n - j + 2):n
# q1 <- min(j * p[i2]/(2:j))
# q[ij] <- pmin(j * p[ij], q1)
# q[i2] <- q[n - j + 1]
# pa <- pmax(pa, q)
# }
# pmax(pa, p)[if (lp < n) ro[1:lp] else ro]
# }, hochberg = {
# i <- lp:1L
# o <- order(p, decreasing = TRUE)
# ro <- order(o)
# pmin(1, cummin((n - i + 1L) * p[o]))[ro]
# }, BH = {
# i <- lp:1L
# o <- order(p, decreasing = TRUE)
# ro <- order(o)
# pmin(1, cummin(n/i * p[o]))[ro]
# }, BY = {
# i <- lp:1L
# o <- order(p, decreasing = TRUE)
# ro <- order(o)
# q <- sum(1L/(1L:n))
# pmin(1, cummin(q * n/i * p[o]))[ro]
# }, none = p)
# p0
#}
#<bytecode: 0x3a61b88>
#<environment: namespace:stats>
 
<lang R>p <- c(4.533744e-01, 7.296024e-01, 9.936026e-02, 9.079658e-02, 1.801962e-01,
8.752257e-01, 2.922222e-01, 9.115421e-01, 4.355806e-01, 5.324867e-01,
4.926798e-01, 5.802978e-01, 3.485442e-01, 7.883130e-01, 2.729308e-01,
Line 2,908 ⟶ 2,847:
 
p.adjust(p, method = 'hommel')
writeLines("Hommel\n")</lang>
</lang>
 
{{out}}
1,336

edits