Holm Error Control Simulation 2x2x2 ANOVA

ErrorControlANOVA.R

library(reshape2)
library(mvtnorm) 
library(ez) 
#Install multtest# try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("multtest")
library(mutoss) #load multiple testing library for Holm function

#2x2x2 within design
N<-50 #sample size per group
r=0.0 #correlation
sd=1 #true standard deviation
nSims<-250000
p3way<-matrix(, nrow = nSims, ncol = 8)
pb <- winProgressBar(title = "progress bar", min = 0, max = nSims, width = 300)
for(i in 1:nSims){ #for each simulated experiment
  setWinProgressBar(pb, i, title=paste( round(i/nSims*100, 2), "% done"))
  a <- as.data.frame(rmvnorm(n=N,mean=c(0,0,0,0),sigma=matrix(c(sd,r,r,r,r,sd,r,r,r,r,sd,r,r,r,r,sd), 4,4))) #change zeroes to introduce a true effect
  b <- as.data.frame(rmvnorm(n=N,mean=c(0,0,0,0),sigma=matrix(c(sd,r,r,r,r,sd,r,r,r,r,sd,r,r,r,r,sd), 4,4))) #change zeroes to introduce a true effect
  c<-cbind(a,b) #combine a and b into one big dataframe
  colnames(c) <- c(1,2,3,4,5,6,7,8) #add column names (within)
  c$subject<-c(1:N) #add subject number to dataframe
  c.long <- melt(c, id.vars = "subject", measure.vars = c(1,2,3,4,5,6,7,8), variable.name = "IV",value.name = "DV")
  c.long$group1<-as.factor(rep(1:2, each=N*8/2))#add groups (between)
  c.long$group2<-as.factor(rep(1:2, each=N*8/4, times = 2))#add groups (between)
  c.long$group3<-as.factor(rep(1:2, each=N*8/8, times = 4))#add groups (between)
  result<-ezANOVA(c.long, dv=.(DV), wid=.(subject), within=.(group1,group2,group3), detailed=TRUE, type="III")
  p3way[i,]<-result$ANOVA$p
}
close(pb)#close progress bar

p3way<-p3way[,-1] #delete first column (p-values for intercept)

#get percentage of p<0.05 (Type 1 errors or power) for the 7 tests
colSums(p3way<0.05)/nrow(p3way)

#apply Holm correction using mutoss package, save adjusted p-values
pholm<-t(apply(p3way, 1, function(x) holm(x, 0.05, silent=TRUE)$adjPValues))
#get percentage of p<0.05 (Type 1 errors or power) for the 7 tests after correction
colSums(pholm<0.05)/nrow(p3way)

#Create plot
library(ggplot2) 
library(extrafont)
font_import(pattern="OpenSans") #Download this font here: https://github.com/google/fonts/blob/d4f76456f96504b75d4707bc3aaa9852b9494b99/apache/opensans/OpenSans-Bold.ttf?raw=true
loadfonts(device="win")

x<-as.vector(p3way)

options(scipen=10000)
FontToUse<-"Open Sans"
SalientLineColor<-"#535353"
LineColor<-"#D0D0D0"
BackgroundColor<-"#F0F0F0"
plot1<-ggplot(as.data.frame(x), aes(x))  + 
  geom_histogram(colour="#535353", fill="#84D5F0", binwidth=0.05) +
  geom_hline(yintercept=37500, colour="#CC79A7", linetype="dashed", size = 1.5) +
  geom_vline(xintercept=0.05, colour="#CC79A7", linetype="dashed", size = 1.5) +
  xlab("\np-value") + ylab("Frequency of p-values\n")  + ggtitle("Holm adjusted p-value distribution") + 
  scale_x_continuous(breaks=c(seq(0,1,0.1))) +
  scale_y_continuous(breaks=c(0,100000,200000,300000,400000,500000,600000,700000, 800000, 900000, 100000)) +
  coord_cartesian(xlim = c(-0.1, 1.1),ylim = c(0, 1000000), expand=FALSE) +
  theme_bw() + 
  theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(),panel.grid.minor.y = element_blank()) + 
  theme(panel.background=element_rect(fill=BackgroundColor)) +
  theme(plot.background=element_rect(fill=BackgroundColor)) +
  theme(panel.border=element_rect(colour=BackgroundColor)) + 
  theme(panel.grid.major=element_line(colour=LineColor,size=.75)) + 
  theme(plot.title=element_text(face="bold",colour=SalientLineColor, vjust=0, size=24,family=FontToUse)) + 
  theme(axis.text.x=element_text(size=15,colour=SalientLineColor, face="bold", family=FontToUse)) +
  theme(axis.text.y=element_text(size=15,colour=SalientLineColor, face="bold", family=FontToUse)) +
  theme(axis.title.y=element_text(size=20,colour=SalientLineColor,face="bold", family=FontToUse)) +
  theme(axis.title.x=element_text(size=20,colour=SalientLineColor,face="bold", family=FontToUse)) + 
  theme(axis.ticks.x=element_line(colour=SalientLineColor, size=2)) +
  theme(axis.ticks.y=element_line(colour=BackgroundColor)) +
  theme(axis.line = element_line()) +
  theme(axis.line.x=element_line(size=1.2,colour=SalientLineColor)) +
  theme(axis.line.y=element_line(colour=BackgroundColor)) + 
  theme(plot.margin = unit(c(1,1,1,1), "cm"))

plot1