c.mco.utf8

Multiple Comparisons

Simulate 50 independent p-values from the null

Some p-values are less than the significance threshold of 0.05. Remember we expect 5% do be below 0.05 by random chance. But you have at most a 5% a priori chance to get a p-value below the Bonferroni threshold.

set.seed(1)
n <- 50
p <- runif(n)
p

##  [1] 0.26550866 0.37212390 0.57285336 0.90820779 0.20168193 0.89838968 0.94467527 0.66079779 0.62911404
## [10] 0.06178627 0.20597457 0.17655675 0.68702285 0.38410372 0.76984142 0.49769924 0.71761851 0.99190609
## [19] 0.38003518 0.77744522 0.93470523 0.21214252 0.65167377 0.12555510 0.26722067 0.38611409 0.01339033
## [28] 0.38238796 0.86969085 0.34034900 0.48208012 0.59956583 0.49354131 0.18621760 0.82737332 0.66846674
## [37] 0.79423986 0.10794363 0.72371095 0.41127443 0.82094629 0.64706019 0.78293276 0.55303631 0.52971958
## [46] 0.78935623 0.02333120 0.47723007 0.73231374 0.69273156

sort(p)

##  [1] 0.01339033 0.02333120 0.06178627 0.10794363 0.12555510 0.17655675 0.18621760 0.20168193 0.20597457
## [10] 0.21214252 0.26550866 0.26722067 0.34034900 0.37212390 0.38003518 0.38238796 0.38410372 0.38611409
## [19] 0.41127443 0.47723007 0.48208012 0.49354131 0.49769924 0.52971958 0.55303631 0.57285336 0.59956583
## [28] 0.62911404 0.64706019 0.65167377 0.66079779 0.66846674 0.68702285 0.69273156 0.71761851 0.72371095
## [37] 0.73231374 0.76984142 0.77744522 0.78293276 0.78935623 0.79423986 0.82094629 0.82737332 0.86969085
## [46] 0.89838968 0.90820779 0.93470523 0.94467527 0.99190609

p<0.05

##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [18] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [35] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE

table(p<0.05)

## 
## FALSE  TRUE 
##    48     2

which(p<0.05)

## [1] 27 47

p[which(p<0.05)]

## [1] 0.01339033 0.02333120

# Bonferroni
0.05/n

## [1] 0.001

which(p<0.05/n)

## integer(0)

Using 47 random and 3 very small p-values

set.seed(2)
psmall <- runif(3,0,1e-3)
psmall

## [1] 0.0001848823 0.0007023740 0.0005733263

p <- c(psmall,runif(47))
p

##  [1] 0.0001848823 0.0007023740 0.0005733263 0.1680519204 0.9438393388 0.9434749587 0.1291589767 0.8334488156
##  [9] 0.4680185155 0.5499837417 0.5526740670 0.2388947594 0.7605133131 0.1808201007 0.4052821812 0.8535484530
## [17] 0.9763984894 0.2258254611 0.4448092291 0.0749794247 0.6618987585 0.3875495426 0.8368891769 0.1505014407
## [25] 0.3472722489 0.4887732314 0.1492468629 0.3570625901 0.9626440483 0.1323720033 0.0104145254 0.1646422420
## [33] 0.8101921445 0.8688610368 0.5142817630 0.6271962866 0.8444290031 0.2848705743 0.6672256477 0.1504697520
## [41] 0.9817278623 0.2970107351 0.1150840754 0.1632008718 0.9440418081 0.7948638236 0.9746878976 0.3490884106
## [49] 0.5019698809 0.8103972627

table(p<0.05)

## 
## FALSE  TRUE 
##    46     4

which(p<0.05)

## [1]  1  2  3 31

p[which(p<0.05)]

## [1] 0.0001848823 0.0007023740 0.0005733263 0.0104145254

# Bonferroni
which(p<0.05/n)

## [1] 1 2 3

p.b <- p.adjust(p,method="bonferroni")
p.b

##  [1] 0.009244113 0.035118702 0.028666317 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000
##  [9] 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000
## [17] 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000
## [25] 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 0.520726270 1.000000000
## [33] 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000
## [41] 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000
## [49] 1.000000000 1.000000000

which(p.b<0.05)

## [1] 1 2 3

Using 40 random and 10 moderately small p-values

set.seed(3)
psmall <- runif(10,0,0.01)
psmall

##  [1] 0.001680415 0.008075164 0.003849424 0.003277343 0.006021007 0.006043941 0.001246334 0.002946009
##  [9] 0.005776099 0.006309793

p <- c(psmall,runif(40))
p

##  [1] 0.001680415 0.008075164 0.003849424 0.003277343 0.006021007 0.006043941 0.001246334 0.002946009
##  [9] 0.005776099 0.006309793 0.512015898 0.505023914 0.534035353 0.557249436 0.867919488 0.829708693
## [17] 0.111449153 0.703688359 0.897488264 0.279732554 0.228201881 0.015329893 0.128981559 0.093381929
## [25] 0.236885007 0.791147409 0.599731566 0.910147711 0.560424554 0.755704769 0.379171891 0.373280977
## [33] 0.170290643 0.453307324 0.258413960 0.336265952 0.889583035 0.201946299 0.579186043 0.207632030
## [41] 0.281468792 0.786281204 0.173019349 0.570747518 0.419282963 0.267622165 0.047809442 0.103493054
## [49] 0.314031456 0.800641062

table(p<0.05)

## 
## FALSE  TRUE 
##    38    12

which(p<0.05)

##  [1]  1  2  3  4  5  6  7  8  9 10 22 47

p[which(p<0.05)]

##  [1] 0.001680415 0.008075164 0.003849424 0.003277343 0.006021007 0.006043941 0.001246334 0.002946009
##  [9] 0.005776099 0.006309793 0.015329893 0.047809442

# Bonferroni
which(p<0.05/n)

## integer(0)

# Benjamini-Hochberg
psort <- sort(p)
psort

##  [1] 0.001246334 0.001680415 0.002946009 0.003277343 0.003849424 0.005776099 0.006021007 0.006043941
##  [9] 0.006309793 0.008075164 0.015329893 0.047809442 0.093381929 0.103493054 0.111449153 0.128981559
## [17] 0.170290643 0.173019349 0.201946299 0.207632030 0.228201881 0.236885007 0.258413960 0.267622165
## [25] 0.279732554 0.281468792 0.314031456 0.336265952 0.373280977 0.379171891 0.419282963 0.453307324
## [33] 0.505023914 0.512015898 0.534035353 0.557249436 0.560424554 0.570747518 0.579186043 0.599731566
## [41] 0.703688359 0.755704769 0.786281204 0.791147409 0.800641062 0.829708693 0.867919488 0.889583035
## [49] 0.897488264 0.910147711

(1:n)/n*0.05

##  [1] 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.010 0.011 0.012 0.013 0.014 0.015 0.016 0.017
## [18] 0.018 0.019 0.020 0.021 0.022 0.023 0.024 0.025 0.026 0.027 0.028 0.029 0.030 0.031 0.032 0.033 0.034
## [35] 0.035 0.036 0.037 0.038 0.039 0.040 0.041 0.042 0.043 0.044 0.045 0.046 0.047 0.048 0.049 0.050

table(psort<(1:n)/n*0.05)

## 
## FALSE  TRUE 
##    41     9

mj <- max(which(psort<(1:n)/n*0.05))
mj

## [1] 10

psort[1:mj]

##  [1] 0.001246334 0.001680415 0.002946009 0.003277343 0.003849424 0.005776099 0.006021007 0.006043941
##  [9] 0.006309793 0.008075164

p.bh <- p.adjust(p,method="BH")
which(p.bh<0.05)

##  [1]  1  2  3  4  5  6  7  8  9 10

p[which(p.bh<0.05)]

##  [1] 0.001680415 0.008075164 0.003849424 0.003277343 0.006021007 0.006043941 0.001246334 0.002946009
##  [9] 0.005776099 0.006309793

sort(p[which(p.bh<0.05)])

##  [1] 0.001246334 0.001680415 0.002946009 0.003277343 0.003849424 0.005776099 0.006021007 0.006043941
##  [9] 0.006309793 0.008075164

Statistics for Laboratory Scientists ( 140.615 )

Multiple Comparisons

Simulate 50 independent p-values from the null

Using 47 random and 3 very small p-values

Using 40 random and 10 moderately small p-values