####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