R code for p curves

I have finally got around to posting the R code for my p curve simulation. Those familiar with R will realize how crude it is (I've been caught up with other urgent stuff and had no time to explore further).

You are welcome to play with (and improve!) the code. Changing delta will alter the (at present) fixed effect size. It would be more realistic to vary this (and the sample sizes). A good starting point for the effect size distribution (in the population) might be a normal distribution with say a mean of zero and a variance of 1 (see Gillett, 1994).

delta <- 0.5
m1 <- 10
sd <- 2
m2 <- m1 + sd*delta
n1 <- n2 <-25
 
n.sims <- 500
p.data <- replicate(n.sims, t.test(rnorm(n1, m1,sd), rnorm(n2, m2,sd))$p.val, simplify=T)
 
par(mfrow=c(5,3))
for (i in 1:15) {
 p.data <- replicate(n.sims, t.test(rnorm(n1, m1,sd), rnorm(n2, m2,sd))$p.val, simplify=T)
 hist(p.data, xlim=c(0,0.1), breaks = 99, col = 'gray')
 }

R code html script courtesy of Pretty R at inside-R.org

References

Gillett, R. (1994). Post Hoc Power Analysis. Journal of Applied Psychology, 79, 783-785.