probabilityofp-values.R

nSims <- 100000 #number of simulated experiments (the more, the more accurate the numbers you get, but the longer it takes. I used 1000000 simulations for my blog)
N<-32 #number of participants
lowp<-0.04
highp<-0.05

#set up some variables
p<-numeric(nSims) 
p2<-numeric(nSims) 
obs_pwr<-numeric(nSims) 
t<-numeric(nSims) 
d<-numeric(nSims) 

for(i in 1:nSims){ #for each simulated experiment
  x<-rnorm(n = N, mean = 100, sd = 20) #produce 100 simulated participants
  y<-rnorm(n = N, mean = 110, sd = 20) #produce 100 simulated participants
  z<-t.test(x,y) #perform the t-test
  p[i]<-z$p.value #get the p-value and store it
}

#Calculate power in simulation
cat ("The power is",(sum(p < 0.05)/nSims*100),"%")

#Count probability of finding p between boundaries specified above

p2<-p[p>lowp & p<=highp]
cat ("The probability of finding a p-value between ",lowp," and ", highp," is ",(length(p2)/nSims*100),"%, which makes it ",((length(p2)/nSims*100)/(((highp-lowp)*100)))," times more probable under the alternative hypothesis than the null-hypothesis (numbers below 1 mean the observed p-value is more likely under the null-hypothesis than under the alternative hypothesis")

#now plot histograms of p-values (the most left bar contains all p-values between 0.00 and 0.05)
hist(p, main="Histogram of p-values", xlab=("Observed p-value"), breaks = 20)