Skip to content

Instantly share code, notes, and snippets.

@Lakens
Last active August 29, 2015 14:17
Show Gist options
  • Save Lakens/799cef8d620413623319 to your computer and use it in GitHub Desktop.
Save Lakens/799cef8d620413623319 to your computer and use it in GitHub Desktop.
probabilities of p-values
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment