Simulate CI

gistfile1.r

#Adapted from on EV Nordheim, MK Clayton & BS Yandell https://www.stat.wisc.edu/~yandell/st571/R/append7.pdf

n.draw = 20 #number of tests you draw
mu = 0.3 #true difference
SD = 1

#determine sample sizes for three tests
library(pwr)
pwr.t.test(d=0.3, power=0.9, sig.level=0.05, type = "one.sample", alternative="two.sided")
pwr.t.test(d=0.3, power=0.9, sig.level=0.01, type = "one.sample", alternative="two.sided") 
pwr.t.test(d=0.3, power=0.9, sig.level=0.001, type = "one.sample", alternative="two.sided")

#draw p=0.05 tests
draws1 = matrix(rnorm(n.draw * 118, mu, SD), 118)
get.conf.int = function(x) t.test(x,conf.level=0.95)$conf.int
conf.int1 = apply(draws1, 2, get.conf.int)
sum(conf.int1[1, ] <= mu & conf.int1[2, ] >= mu)
sum(conf.int1[1, ] > 0)

#draw p=0.01 tests
draws2 = matrix(rnorm(n.draw * 167, mu, SD), 167)
get.conf.int = function(x) t.test(x,conf.level=0.95)$conf.int
conf.int2 = apply(draws2, 2, get.conf.int)
sum(conf.int2[1, ] <= mu & conf.int2[2, ] >= mu)
sum(conf.int2[1, ] > 0)

#draw p=0.001 tests
draws3 = matrix(rnorm(n.draw * 238, mu, SD), 238)
get.conf.int = function(x) t.test(x,conf.level=0.95)$conf.int
conf.int3 = apply(draws3, 2, get.conf.int)
sum(conf.int3[1, ] <= mu & conf.int3[2, ] >= mu)
sum(conf.int3[1, ] > 0)

#combine draws
conf.int.tot<-cbind(conf.int1,conf.int2,conf.int3)

#plot CI
plot(range(0, conf.int.tot), c(0, 1 + n.draw*3), type = "n", xlab = "mean", ylab = "sample run")
for (i in 1:20) lines(conf.int.tot[, i], rep(i, 2), lwd = 2, col = "gray48")
for (i in 21:40) lines(conf.int.tot[, i], rep(i, 2), lwd = 2, col = "gray35")
for (i in 41:60) lines(conf.int.tot[, i], rep(i, 2), lwd = 2, col = "gray8")
abline(v = mu, lwd = 2, lty = 2)
abline(v = 0, lwd = 2, col = "lightgray", lty = 3)