####NULL hypotheses are all true
set.seed(1010093)
pValues <- rep(NA,1000)
for(i in 1:1000){
  y <- rnorm(20)
  x <- rnorm(20)
  pValues[i] <- summary(lm(y ~ x))$coeff[2,4]
}

# Controls type I error
sum(pValues < 0.05)


# Controls FWER 
sum(p.adjust(pValues,method="bonferroni") < 0.05)


###we could do more times

ben.rejct = c()
for(j in 1:100){
pValues <- rep(NA,1000)
	for(i in 1:1000){
	  y <- rnorm(20)
	  x <- rnorm(20)
	  pValues[i] <- summary(lm(y ~ x))$coeff[2,4]
	}
	tmp = sum(p.adjust(pValues,method="bonferroni") < 0.05)
	ben.rejct = c(ben.rejct,tmp)
	}

##control FWER by hochberg
sum(p.adjust(pValues,method="hochberg") < 0.05)

# Controls FDR 
sum(p.adjust(pValues,method="BH") < 0.05)



####50% of NULL hypothesis are true
set.seed(1010093)
pValues <- rep(NA,1000)
for(i in 1:1000){
  x <- rnorm(20)
  # First 500 beta=0, last 500 beta=2
  if(i <= 500){y <- rnorm(20)}else{ y <- rnorm(20,mean=2*x)}
  pValues[i] <- summary(lm(y ~ x))$coeff[2,4]
}
trueStatus <- rep(c("zero","notZero"),each=500)
tmp = table(pValues < 0.05, trueStatus)
rownames(tmp) = c("retain","reject")
tmp


# Controls FWER 
tmp = table(p.adjust(pValues,method="bonferroni") < 0.05,trueStatus)
rownames(tmp) = c("retain","reject")
tmp


tmp = table(p.adjust(pValues,method="hommel") < 0.05,trueStatus)
rownames(tmp) = c("retain","reject")
tmp


# Controls FDR 
tmp = table(p.adjust(pValues,method="BH") < 0.05,trueStatus)
rownames(tmp) = c("retain","reject")
tmp

# Controls FDR
tmp = table(p.adjust(pValues,method="BY") < 0.05,trueStatus)
rownames(tmp) = c("retain","reject")
tmp

par(mfrow=c(1,2))
plot(pValues,p.adjust(pValues,method="bonferroni"),pch=19)
plot(pValues,p.adjust(pValues,method="BH"),pch=19)



par(mfrow=c(1,2))
plot(pValues,p.adjust(pValues,method="BH"),pch=19)
plot(pValues,p.adjust(pValues,method="BY"),pch=19)

par(mfrow=c(1,2))
plot(pValues,p.adjust(pValues,method="bonferroni"),pch=19)
plot(pValues,p.adjust(pValues,method="hommel"),pch=19)



p1 = p.adjust(pValues,method="bonferroni") 
p2 = p.adjust(pValues,method="hommel")

sum(p2<p1)