test.R

test.R

if(!require(ggplot2)){install.packages('ggplot2')}
library(ggplot2)
if(!require(pwr)){install.packages('pwr')}
library(pwr)

nSims <- 100000 #number of simulated experiments

M<-106 #Mean IQ score in the sample
n<-26 #set sample size
SD<-15 #SD of the simulated data
  
p <-numeric(nSims) #set up empty container for all simulated p-values
binwidth = 0.05
#Run simulation
for(i in 1:nSims){ #for each simulated experiment
  x<-rnorm(n = n, mean = M, sd = SD) #Simulate data with specified mean, standard deviation, and sample size
  z<-t.test(x, mu=100) #perform the t-test against mu (set to value you want to test against)
  p[i]<-z$p.value #get the p-value and store it
}

#Check power by summing significant p-values and dividing by number of simulations
(sum(p < 0.05)/nSims) #power
#Calculate power formally by power analysis
power<-pwr.t.test(d=(M-100)/SD, n=n,sig.level=0.05,type="one.sample",alternative="two.sided")$power #determines M when power > 0. When power = 0, will set  M = 100.

#plot pvalue distribution
#png(file="PvalueDistribution.png",width=4000,height=2000, res = 300)
ggplot(as.data.frame(p), aes(p))  + 
  geom_histogram(colour="black", fill="grey", binwidth = binwidth) +
  xlab("P-values") + ylab("number of p-values")  + ggtitle(paste("P-value Distribution with",round(power*100, digits=1),"% Power")) + theme_bw(base_size=20) + 
  scale_x_continuous(breaks=seq(0,1,0.1)) + coord_cartesian(xlim = c(0, 1), ylim = c(0, 100000)) +
  geom_hline(yintercept=nSims*binwidth, color="red", linetype="dashed", size=1.2) +
  theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank())
#dev.off()