R - multiple comparisons

along with chapter 16 of Draghici: Statistics and Data Analysis for Microarrays Using R and Bioconductor

• main function in R to obtain corrected p-vals p.adjust
  
• to get available methods type:
  

p.adjust.methods

[1] "holm"       "hochberg"   "hommel"     "bonferroni" "BH"        
[6] "BY"         "fdr"        "none"

• simulate 10 raw p-vals on gene level
  

raw.p.vals <- runif(10,0,0.5)
names(raw.p.vals) <- paste("gene",1:length(raw.p.vals),sep="")
raw.p.vals

     gene1      gene2      gene3      gene4      gene5      gene6      gene7 
0.30128864 0.45386933 0.11451989 0.17423164 0.36300090 0.32262681 0.24504471 
     gene8      gene9     gene10 
0.27808264 0.08773618 0.49969999

• get the corrected vals for example
  
  • Holm's stepwise correction
    
  
  

p.adjust(raw.p.vals,method="holm")

    gene1     gene2     gene3     gene4     gene5     gene6     gene7     gene8 
1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 
    gene9    gene10 
0.8773618 1.0000000

• False Discovery Rate
  

p.adjust(raw.p.vals,method="fdr")

    gene1     gene2     gene3     gene4     gene5     gene6     gene7     gene8 
0.4537511 0.4997000 0.4537511 0.4537511 0.4537511 0.4537511 0.4537511 0.4537511 
    gene9    gene10 
0.4537511 0.4997000

• the multtest package (available from bioconductor) contains also a function for correction (mt.rawp2adjp())
  

require(multtest)
mymethods <- c("Bonferroni", "Holm", "Hochberg", "SidakSS", "SidakSD","BH", "BY","ABH","TSBH")
mt.rawp2adjp(raw.p.vals,proc=mymethods)

$adjp
            rawp Bonferroni      Holm Hochberg   SidakSS   SidakSD        BH BY
 [1,] 0.08773618  0.8773618 0.8773618   0.4997 0.6007872 0.6007872 0.4537511  1
 [2,] 0.11451989  1.0000000 1.0000000   0.4997 0.7036615 0.6653358 0.4537511  1
 [3,] 0.17423164  1.0000000 1.0000000   0.4997 0.8525712 0.7837949 0.4537511  1
 [4,] 0.24504471  1.0000000 1.0000000   0.4997 0.9398532 0.8602188 0.4537511  1
 [5,] 0.27808264  1.0000000 1.0000000   0.4997 0.9615519 0.8602188 0.4537511  1
 [6,] 0.30128864  1.0000000 1.0000000   0.4997 0.9722682 0.8602188 0.4537511  1
 [7,] 0.32262681  1.0000000 1.0000000   0.4997 0.9796633 0.8602188 0.4537511  1
 [8,] 0.36300090  1.0000000 1.0000000   0.4997 0.9890001 0.8602188 0.4537511  1
 [9,] 0.45386933  1.0000000 1.0000000   0.4997 0.9976397 0.8602188 0.4997000  1
[10,] 0.49969999  1.0000000 1.0000000   0.4997 0.9990176 0.8602188 0.4997000  1
      ABH TSBH_0.05
 [1,]  NA 0.4537511
 [2,]  NA 0.4537511
 [3,]  NA 0.4537511
 [4,]  NA 0.4537511
 [5,]  NA 0.4537511
 [6,]  NA 0.4537511
 [7,]  NA 0.4537511
 [8,]  NA 0.4537511
 [9,]  NA 0.4997000
[10,]  NA 0.4997000

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

$h0.ABH
[1] NA

$h0.TSBH
h0.TSBH_0.05 
          10 

Warnmeldung:
In min(which(diff(h0.m, na.rm = TRUE) > 0), na.rm = TRUE) :
  kein nicht-fehlendes Argument für min; gebe Inf zurück