Simulate Pdistribution ML Fig4

gistfile1.r

nSims <- 11000 #number of simulated experiments
psig <-numeric(nSims) #set up empty container for all simulated p-values
pnonsig <-numeric(nSims) #set up empty container for all simulated p-values
z <-numeric(nSims) #set up empty container for z-scores

for(i in 1:nSims){ #for each simulated experiment
  x<-rnorm(n = 38, mean = 100, sd = 20) #produce  simulated participants
  y<-rnorm(n = 38, mean = 108, sd = 20) #produce  simulated participants
  t<-t.test(x,y) #perform the t-test
  if(t$p.value<=0.05&t$p.value>=0.01) {
    psig[i]<-t$p.value #get the p-value and store it
  }
  if(t$p.value>0.05&t$p.value<=0.1) {
    pnonsig[i]<-t$p.value #get the p-value and store it
  }
}

psig<-psig[psig>0.0000000001] #remove empty cells of studies due to publication bias
pnonsig<-pnonsig[pnonsig>0.0000000001] #remove empty cells of studies due to publication bias
length(psig)
length(pnonsig)
#mimick publication bias by removing 50% of non-significant values
pnonsig<-head(pnonsig, length(pnonsig)/2)
length(pnonsig)


#Simulate the Type 1 errors
nSims <- 4500 #number of simulated experiments
p <-numeric(nSims) #set up empty container for all simulated p-values

for(i in 1:nSims){ #for each simulated experiment
  x<-rnorm(n = 23, mean = 100, sd = 20) #produce  simulated participants
  y1<-rnorm(n = 23, mean = 100, sd = 20) #produce  simulated participants
  y2<-rnorm(n = 23, mean = 100, sd = 20) #produce  simulated participants
  y3<-rnorm(n = 23, mean = 100, sd = 20) #produce  simulated participants
  y4<-rnorm(n = 23, mean = 100, sd = 20) #produce  simulated participants
  y5<-rnorm(n = 23, mean = 100, sd = 20) #produce  simulated participants
  t<-t.test(x,y1) #perform the t-test
  if(t$p.value<0.05) {
    p[i]<-t$p.value #get the p-value and store it
  }
  if(t$p.value>0.05) {
    t<-t.test(x,y2) #perform the second t-test
    if(t$p.value<0.05) {
      p[i]<-t$p.value #get the p-value and store it
    }
  }
  if(t$p.value>0.05) {
    t<-t.test(x,y3) #perform the third t-test
    if(t$p.value<0.05) {
      p[i]<-t$p.value #get the p-value and store it
    }
  }
  if(t$p.value>0.05) {
    t<-t.test(x,y4) #perform the third t-test
    if(t$p.value<0.05) {
      p[i]<-t$p.value #get the p-value and store it
    }
  }
  if(t$p.value>0.05) {
    t<-t.test(x,y5) #perform the fifth t-test
    p[i]<-t$p.value #get the p-value (regardless of whether significant and store it
  }
}
pmult<-p[p<0.1]
pmult<-pmult[pmult>0.01]
length(pmult)
hist(pmult, main="Testing exp. condition against 1/5 control conditions", xlab=("Observed p-value"), breaks =10)




ptot<-c(psig,pnonsig,pmult)

hist(ptot, xaxt="n", main="Histogram of p-values", xlab=("Observed p-value"), breaks =20)
axis(side = 1, at = c(0.01, 0.015, 0.02, 0.025, 0.03, 0.035, 0.04, 0.045, 0.05,0.055,0.06,0.065,0.07,0.075,0.08,0.085,0.09,0.095,0.1))
axis(side = 2, at = c(1000,2000,3000,4000,5000,6000))
abline(v=0.05, lwd=2)

length(ptot)