ベタなシミュレーション

実際の研究の手順を再現してるので、シミュレーションとしては無駄が多いです。

set.seed(42)

ある実験(t検定のp値)を表現。

f = function(m) t.test(rnorm(32, 0, 1), rnorm(32, m, 1))$p.value

alpha = 0.05で。

alpha = 0.05
# 最初の研究 (100個)
es = rnorm(100, 0, 0.5) # 母集団の平均値の差のこと
ps = sapply(es, f)
# p値
plot(ps, col=ifelse(ps<alpha, "red", "black")); abline(h=alpha)

# 追試(有意だった研究に対して、100回ずつ)
rps = sapply(es[ps<alpha], function(s) {r = sapply(1:100, function(x) f(s))})
# 追試が成功してる確率
plot(c(rps), col=ifelse(rps<alpha, "red", "black")); abline(h=alpha)