library("multtest", lib.loc = "C:/Rpackages/multtest_2.12.0")

?mt.rawp2adjp

pvals <- c(0.01, 0.02, 0.005, 0.13, 0.045) 

## No correction: Reject 4 Nullhyps. at 5% level

## FWER: Bonferroni 
## (control at 5% level; i.e., prob. of having one or more false pos. is at most 5%)
resB <- mt.rawp2adjp(pvals, proc = "Bonferroni")
resB
rejectB <- resB$index[resB$adjp[,2] < 0.05]
rejectB ## reject only one Nullhyp.

## FWER: Holm-Bonferroni
## (control at 5% level; i.e., prob. of having one or more false pos. is at most 5%)
resH <- mt.rawp2adjp(pvals, proc = "Holm")
resH
rejectH <- resH$index[resH$adjp[,2] < 0.05]
rejectH ## reject two Nullhyps.: Holm is more powerful than Bonferroni

## FDR: Benjamini-Hochberg 
## (control at 10% level; i.e., about 10% of the found are actually not there)
resBH <- mt.rawp2adjp(pvals, proc = "BH")
resBH
rejectBH <- resBH$index[resBH$adjp[,2] < 0.05]
rejectBH ## reject three Nullhyps.: Holm is more powerful than Bonferroni


#######################
## Case study: Which genes are associated with different types of Leukemia?
#######################
data(golub)
?golub

head(golub)
golub.cl

teststat <- mt.teststat(golub, golub.cl)
str(teststat)

rawp0 <- 2 * (1 - pnorm(abs(teststat))) ## use normals distr. instead of t as simplification
image(golub)
## no correction
sum(rawp0 < 0.05)
## FDR < 10%
resBH <- mt.rawp2adjp(rawp0, proc = "BH")
sum(resBH$adjp[,"BH"] < 0.1)
## FWER < 5% using Holm
resBH <- mt.rawp2adjp(rawp0, proc = "Holm")
sum(resBH$adjp[,"Holm"] < 0.05)