Skip to content

Commit a581735

Browse files
author
maechler
committed
say more about "BY" -- fixing PR#17136
git-svn-id: https://svn.r-project.org/R/trunk@89312 00db46b3-68df-0310-9c12-caf00c1e9a41
1 parent 5481640 commit a581735

File tree

2 files changed

+116
-76
lines changed

2 files changed

+116
-76
lines changed

‎src/library/stats/man/p.adjust.Rd‎

Lines changed: 48 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
% File src/library/stats/man/p.adjust.Rd
22
% Part of the R package, https://www.R-project.org
3-
% Copyright 1995-2025 R Core Team
3+
% Copyright 1995-2026 R Core Team
44
% Distributed under GPL 2 or later
55

66
\name{p.adjust}
@@ -55,10 +55,18 @@ p.adjust.methods
5555
5656
The \code{"BH"} (aka \code{"fdr"}) and \code{"BY"} methods of
5757
\I{Benjamini}, \I{Hochberg}, and \I{Yekutieli} control the false discovery rate,
58-
the expected proportion of false discoveries amongst the rejected
59-
hypotheses. The false discovery rate is a less stringent condition
60-
than the family-wise error rate, so these methods are more powerful
61-
than the others.
58+
\I{FDR}, the expected proportion of false discoveries amongst the
59+
rejected hypotheses.
60+
61+
The \code{"BY"} correction modifies the \I{BH} procedure by replacing the
62+
target level \eqn{q} with \eqn{q / \sum_{i=1}^{m} 1/i}, where \eqn{m} is
63+
the number of tests (Theorem 1.3 in the reference), which controls the
64+
\I{FDR} under the most general form of dependence structure. This will
65+
be more conservative than \code{"BH"}, for small \code{p} even more than
66+
\I{Bonferroni}, see the example.
67+
The \I{FDR} as implemented by the \code{"BH"} method is a less stringent
68+
condition than the family-wise error rate, so it is typically more
69+
powerful than the others.
6270
6371
Note that you can set \code{n} larger than \code{length(p)} which
6472
means the unobserved p-values are assumed to be greater than all the
@@ -99,15 +107,44 @@ p.adjust.M <- p.adjust.methods[p.adjust.methods != "fdr"]
99107
p.adj <- sapply(p.adjust.M, function(meth) p.adjust(p, meth))
100108
p.adj.60 <- sapply(p.adjust.M, function(meth) p.adjust(p, meth, n = 60))
101109
stopifnot(identical(p.adj[,"none"], p), p.adj <= p.adj.60)
110+
102111
round(p.adj, 3)
103112
## or a bit nicer:
104-
noquote(apply(p.adj, 2, format.pval, digits = 3))
105-
106-
107-
## and a graphic:
113+
head(round(100 * p.adj[,c(7,1:6)], 2), n=21) # in [percent]:
114+
## none holm hochberg hommel bonferroni BH BY
115+
## [1,] 0.00 0.00 0.00 0.00 0.00 0.00 0.01 *)
116+
## [2,] 0.00 0.10 0.10 0.10 0.11 0.04 0.19 *)
117+
## [3,] 0.00 0.12 0.12 0.12 0.13 0.04 0.19 *)
118+
## [4,] 0.01 0.46 0.46 0.42 0.49 0.09 0.43
119+
## [5,] 0.01 0.48 0.48 0.45 0.53 0.09 0.43
120+
## .... .......... ............ ...............
121+
## .... .......... ............ ................
122+
## [18,] 0.88 29.06 29.06 27.30 44.03 2.45 11.00
123+
## [19,] 0.94 30.08 30.08 29.14 47.01 2.47 11.13
124+
## [20,] 1.13 35.02 35.02 33.89 56.49 2.82 12.71
125+
## [21,] 2.12 63.45 63.45 57.11 100.00 5.04 22.66
126+
##
127+
## *) The smallest 3 Bonferroni values are smaller than the "BY" ones,
128+
## (John Maindonald, PR#17136)
129+
130+
## number of rejected H0 ("P" < 0.05):
131+
colSums(p.adj < 0.05)
132+
## holm hochberg hommel bonferroni BH BY none
133+
## 11 11 11 11 20 12 22
134+
135+
## visual comparison
108136
matplot(p, p.adj, ylab="p.adjust(p, meth)", type = "l", asp = 1, lty = 1:6,
109-
main = "P-value adjustments")
110-
legend(0.7, 0.6, p.adjust.M, col = 1:6, lty = 1:6)
137+
col = 1:7, main = "P-value adjustments")
138+
legR <- function() {
139+
legend("bottomright", p.adjust.M, col = 1:7, lty = 1:6, bty = "n", inset = 0.05)
140+
rug(p) }
141+
legR()
142+
143+
## zoom in & log scale
144+
lim <- c(7e-7, .20)
145+
matplot(p, p.adj, ylab="p.adjust(p, meth)", type = "l", asp = 1, lty = 1:6, col = 1:7,
146+
main = "P-value adjustments [log-log]", log = "xy", xlim=lim, ylim=lim, las=1)
147+
legR()
111148
112149
## Can work with NAs:
113150
pN <- p; iN <- c(46, 47); pN[iN] <- NA

‎tests/Examples/stats-Ex.Rout.save‎

Lines changed: 68 additions & 65 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11

2-
R Under development (unstable) (2025-12-26 r89234) -- "Unsuffered Consequences"
3-
Copyright (C) 2025 The R Foundation for Statistical Computing
2+
R Under development (unstable) (2026-01-21 r89311) -- "Unsuffered Consequences"
3+
Copyright (C) 2026 The R Foundation for Statistical Computing
44
Platform: x86_64-pc-linux-gnu
55

66
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -3353,7 +3353,7 @@ Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm) :
33533353
[1] 2 2 3 6 6 6 5 5 3 3 3
33543354
>
33553355
> ## 'ties' + NAs -- notably NAs for tied x[], situation as PR#17604
3356-
>
3356+
>
33573357
> x <- c(2:3, 3:5, 5:7)
33583358
> y <- c(1,NA,2:4,NA,1,0)
33593359
> ## allapprox() [defined above] for all variants :
@@ -8370,10 +8370,10 @@ Residual Deviance: 0.01267 AIC: 27.03
83708370
> ## tree from the cluster centers.
83718371
> hc <- hclust(dist(USArrests)^2, "cen")
83728372
> memb <- cutree(hc, k = 10)
8373-
> cent <- NULL
8374-
> for(k in 1:10){
8375-
+ cent <- rbind(cent, colMeans(USArrests[memb == k, , drop = FALSE]))
8376-
+ }
8373+
> cent <- matrix(numeric(), 10, 4)
8374+
> for(k in 1:10)
8375+
+ cent[k,] <- colMeans(USArrests[memb == k, , drop = FALSE])
8376+
>
83778377
> hc1 <- hclust(dist(cent)^2, method = "cen", members = table(memb))
83788378
> opar <- par(mfrow = c(1, 2))
83798379
> plot(hc, labels = FALSE, hang = -1, main = "Original Tree")
@@ -11875,6 +11875,7 @@ $objective
1187511875
> p.adj <- sapply(p.adjust.M, function(meth) p.adjust(p, meth))
1187611876
> p.adj.60 <- sapply(p.adjust.M, function(meth) p.adjust(p, meth, n = 60))
1187711877
> stopifnot(identical(p.adj[,"none"], p), p.adj <= p.adj.60)
11878+
>
1187811879
> round(p.adj, 3)
1187911880
holm hochberg hommel bonferroni BH BY none
1188011881
[1,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000
@@ -11928,64 +11929,66 @@ $objective
1192811929
[49,] 1.000 0.944 0.944 1.000 0.930 1.000 0.912
1192911930
[50,] 1.000 0.944 0.944 1.000 0.944 1.000 0.944
1193011931
> ## or a bit nicer:
11931-
> noquote(apply(p.adj, 2, format.pval, digits = 3))
11932-
holm hochberg hommel bonferroni BH BY none
11933-
[1,] 1.18e-05 1.18e-05 1.18e-05 1.18e-05 1.18e-05 5.3e-05 2.35e-07
11934-
[2,] 0.00103 0.00103 0.00101 0.00105 0.000429 0.00193 2.10e-05
11935-
[3,] 0.00124 0.00124 0.00124 0.00129 0.000429 0.00193 2.58e-05
11936-
[4,] 0.00461 0.00461 0.00422 0.00491 0.000947 0.00426 9.81e-05
11937-
[5,] 0.00484 0.00484 0.00453 0.00526 0.000947 0.00426 0.000105
11938-
[6,] 0.00559 0.00559 0.00521 0.00621 0.000947 0.00426 0.000124
11939-
[7,] 0.00583 0.00583 0.00557 0.00663 0.000947 0.00426 0.000133
11940-
[8,] 0.00674 0.00674 0.00659 0.00784 0.000980 0.00441 0.000157
11941-
[9,] 0.00947 0.00947 0.00924 0.01127 0.001253 0.00564 0.000225
11942-
[10,] 0.01556 0.01556 0.01518 0.01898 0.001898 0.00854 0.000380
11943-
[11,] 0.02446 0.02446 0.02446 0.03057 0.002780 0.01251 0.000611
11944-
[12,] 0.06294 0.06294 0.05810 0.08070 0.006725 0.03026 0.001614
11945-
[13,] 0.12549 0.12549 0.10898 0.16512 0.012637 0.05686 0.003302
11946-
[14,] 0.13092 0.13092 0.11677 0.17692 0.012637 0.05686 0.003538
11947-
[15,] 0.18853 0.18853 0.16758 0.26185 0.017457 0.07854 0.005237
11948-
[16,] 0.23912 0.23912 0.21179 0.34160 0.020762 0.09341 0.006832
11949-
[17,] 0.24001 0.24001 0.21884 0.35296 0.020762 0.09341 0.007059
11950-
[18,] 0.29057 0.29057 0.27296 0.44026 0.024459 0.11004 0.008805
11951-
[19,] 0.30083 0.30083 0.29143 0.47005 0.024740 0.11131 0.009401
11952-
[20,] 0.35024 0.35024 0.33894 0.56490 0.028245 0.12708 0.011298
11953-
[21,] 0.63451 0.63451 0.57105 1.00000 0.050358 0.22657 0.021150
11954-
[22,] 1.00000 0.94379 0.94379 1.00000 0.111880 0.50337 0.049227
11955-
[23,] 1.00000 0.94379 0.94379 1.00000 0.130463 0.58698 0.060533
11956-
[24,] 1.00000 0.94379 0.94379 1.00000 0.130463 0.58698 0.062622
11957-
[25,] 1.00000 0.94379 0.94379 1.00000 0.147903 0.66545 0.073952
11958-
[26,] 1.00000 0.94379 0.94379 1.00000 0.159252 0.71651 0.082811
11959-
[27,] 1.00000 0.94379 0.94379 1.00000 0.159877 0.71932 0.086333
11960-
[28,] 1.00000 0.94379 0.94379 1.00000 0.212617 0.95661 0.119065
11961-
[29,] 1.00000 0.94379 0.94379 1.00000 0.325999 1.00000 0.189080
11962-
[30,] 1.00000 0.94379 0.94379 1.00000 0.343082 1.00000 0.205849
11963-
[31,] 1.00000 0.94379 0.94379 1.00000 0.356325 1.00000 0.220921
11964-
[32,] 1.00000 0.94379 0.94379 1.00000 0.446250 1.00000 0.285600
11965-
[33,] 1.00000 0.94379 0.94379 1.00000 0.461954 1.00000 0.304889
11966-
[34,] 1.00000 0.94379 0.94379 1.00000 0.683577 1.00000 0.466068
11967-
[35,] 1.00000 0.94379 0.94379 1.00000 0.683577 1.00000 0.483081
11968-
[36,] 1.00000 0.94379 0.94379 1.00000 0.683577 1.00000 0.492175
11969-
[37,] 1.00000 0.94379 0.94379 1.00000 0.718845 1.00000 0.531945
11970-
[38,] 1.00000 0.94379 0.94379 1.00000 0.741435 1.00000 0.575155
11971-
[39,] 1.00000 0.94379 0.94379 1.00000 0.741435 1.00000 0.578319
11972-
[40,] 1.00000 0.94379 0.94379 1.00000 0.762606 1.00000 0.618589
11973-
[41,] 1.00000 0.94379 0.94379 1.00000 0.762606 1.00000 0.636362
11974-
[42,] 1.00000 0.94379 0.94379 1.00000 0.762606 1.00000 0.644859
11975-
[43,] 1.00000 0.94379 0.94379 1.00000 0.762606 1.00000 0.655841
11976-
[44,] 1.00000 0.94379 0.94379 1.00000 0.782487 1.00000 0.688588
11977-
[45,] 1.00000 0.94379 0.94379 1.00000 0.798874 1.00000 0.718986
11978-
[46,] 1.00000 0.94379 0.94379 1.00000 0.880265 1.00000 0.817954
11979-
[47,] 1.00000 0.94379 0.94379 1.00000 0.880265 1.00000 0.827449
11980-
[48,] 1.00000 0.94379 0.94379 1.00000 0.930478 1.00000 0.897130
11981-
[49,] 1.00000 0.94379 0.94379 1.00000 0.930478 1.00000 0.911868
11982-
[50,] 1.00000 0.94379 0.94379 1.00000 0.943789 1.00000 0.943789
11983-
>
11984-
>
11985-
> ## and a graphic:
11932+
> head(round(100 * p.adj[,c(7,1:6)], 2), n=21) # in [percent]:
11933+
none holm hochberg hommel bonferroni BH BY
11934+
[1,] 0.00 0.00 0.00 0.00 0.00 0.00 0.01
11935+
[2,] 0.00 0.10 0.10 0.10 0.11 0.04 0.19
11936+
[3,] 0.00 0.12 0.12 0.12 0.13 0.04 0.19
11937+
[4,] 0.01 0.46 0.46 0.42 0.49 0.09 0.43
11938+
[5,] 0.01 0.48 0.48 0.45 0.53 0.09 0.43
11939+
[6,] 0.01 0.56 0.56 0.52 0.62 0.09 0.43
11940+
[7,] 0.01 0.58 0.58 0.56 0.66 0.09 0.43
11941+
[8,] 0.02 0.67 0.67 0.66 0.78 0.10 0.44
11942+
[9,] 0.02 0.95 0.95 0.92 1.13 0.13 0.56
11943+
[10,] 0.04 1.56 1.56 1.52 1.90 0.19 0.85
11944+
[11,] 0.06 2.45 2.45 2.45 3.06 0.28 1.25
11945+
[12,] 0.16 6.29 6.29 5.81 8.07 0.67 3.03
11946+
[13,] 0.33 12.55 12.55 10.90 16.51 1.26 5.69
11947+
[14,] 0.35 13.09 13.09 11.68 17.69 1.26 5.69
11948+
[15,] 0.52 18.85 18.85 16.76 26.18 1.75 7.85
11949+
[16,] 0.68 23.91 23.91 21.18 34.16 2.08 9.34
11950+
[17,] 0.71 24.00 24.00 21.88 35.30 2.08 9.34
11951+
[18,] 0.88 29.06 29.06 27.30 44.03 2.45 11.00
11952+
[19,] 0.94 30.08 30.08 29.14 47.01 2.47 11.13
11953+
[20,] 1.13 35.02 35.02 33.89 56.49 2.82 12.71
11954+
[21,] 2.12 63.45 63.45 57.11 100.00 5.04 22.66
11955+
> ## none holm hochberg hommel bonferroni BH BY
11956+
> ## [1,] 0.00 0.00 0.00 0.00 0.00 0.00 0.01 *)
11957+
> ## [2,] 0.00 0.10 0.10 0.10 0.11 0.04 0.19 *)
11958+
> ## [3,] 0.00 0.12 0.12 0.12 0.13 0.04 0.19 *)
11959+
> ## [4,] 0.01 0.46 0.46 0.42 0.49 0.09 0.43
11960+
> ## [5,] 0.01 0.48 0.48 0.45 0.53 0.09 0.43
11961+
> ## .... .......... ............ ...............
11962+
> ## .... .......... ............ ................
11963+
> ## [18,] 0.88 29.06 29.06 27.30 44.03 2.45 11.00
11964+
> ## [19,] 0.94 30.08 30.08 29.14 47.01 2.47 11.13
11965+
> ## [20,] 1.13 35.02 35.02 33.89 56.49 2.82 12.71
11966+
> ## [21,] 2.12 63.45 63.45 57.11 100.00 5.04 22.66
11967+
> ##
11968+
> ## *) The smallest 3 Bonferroni values are smaller than the "BY" ones,
11969+
> ## (John Maindonald, PR#17136)
11970+
>
11971+
> ## number of rejected H0 ("P" < 0.05):
11972+
> colSums(p.adj < 0.05)
11973+
holm hochberg hommel bonferroni BH BY none
11974+
11 11 11 11 20 12 22
11975+
> ## holm hochberg hommel bonferroni BH BY none
11976+
> ## 11 11 11 11 20 12 22
11977+
>
11978+
> ## visual comparison
1198611979
> matplot(p, p.adj, ylab="p.adjust(p, meth)", type = "l", asp = 1, lty = 1:6,
11987-
+ main = "P-value adjustments")
11988-
> legend(0.7, 0.6, p.adjust.M, col = 1:6, lty = 1:6)
11980+
+ col = 1:7, main = "P-value adjustments")
11981+
> legR <- function() {
11982+
+ legend("bottomright", p.adjust.M, col = 1:7, lty = 1:6, bty = "n", inset = 0.05)
11983+
+ rug(p) }
11984+
> legR()
11985+
>
11986+
> ## zoom in & log scale
11987+
> lim <- c(7e-7, .20)
11988+
> matplot(p, p.adj, ylab="p.adjust(p, meth)", type = "l", asp = 1, lty = 1:6, col = 1:7,
11989+
+ main = "P-value adjustments [log-log]", log = "xy", xlim=lim, ylim=lim, las=1)
11990+
> legR()
11991+
Warning in rug(p) : some values will be clipped
1198911992
>
1199011993
> ## Can work with NAs:
1199111994
> pN <- p; iN <- c(46, 47); pN[iN] <- NA
@@ -19682,7 +19685,7 @@ Number of Fisher Scoring iterations: 6
1968219685
> cleanEx()
1968319686
> options(digits = 7L)
1968419687
> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
19685-
Time elapsed: 5.797 0.16 5.961 0 0
19688+
Time elapsed: 3.598 0.325 4.098 0 0
1968619689
> grDevices::dev.off()
1968719690
null device
1968819691
1

0 commit comments

Comments
 (0)