Simulate Z-scores v2

z-scores

#####EFFECT WITH PUBLICATION BIAS#########

nSims <- 10000 #number of simulated experiments
z1 <-numeric(nSims) #set up empty container for z-scores

for(i in 1:nSims){ #for each simulated experiment
  x<-rnorm(n = 23, mean = 100, sd = 20) #produce  simulated participants
  y<-rnorm(n = 23, mean = 114, sd = 20) #produce  simulated participants
  t<-t.test(x,y) #perform the t-test
  if(t$p.value<0.05) {
    z1[i]<-qnorm(1-(t$p.value/2)) #convert to the z-score and store it
  }  
}


z1<-z1[z1>0.0001] #remove empty cells of studies due to publication bias


#####EFFECT WITHOUT PUBLICATION BIAS#########

nSims <- 10000 #number of simulated experiments

z2 <-numeric(nSims) #set up empty container for z-scores

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


#####NO EFFECT WITHOUT PUBLICATION BIAS#########

nSims <- 10000 #number of simulated experiments
z3 <-numeric(nSims) #set up empty container for z-scores

for(i in 1:nSims){ #for each simulated experiment
  x<-rnorm(n = 23, mean = 100, sd = 40) #produce  simulated participants
  y<-rnorm(n = 23, mean = 100, sd = 40) #produce  simulated participants
  t<-t.test(x,y) #perform the t-test
  z3[i]<-qnorm(1-(t$p.value/2)) #convert to the z-score and store it
}

z<-c(z1,z2,z3)
#now plot the histogram NOTE if you get error some 'x' not counted; maybe 'breaks' do not span range of 'x' increase max scale end from 8 to higher value
hist(z, main="Histogram of z-scores", xlab=("Observed z-score"), seq(0,8,by=0.245))