Post-hoc power analysis

posthoc_power_analysis.R

mu <- 0.1
sigma <- 1
nsamp <- 36
B <- 1000

real_power <- power.t.test(delta = mu, sd = sigma, n = nsamp, 
                           alternative = "one.sided",
                           type = "one.sample")$power

res <- replicate(B, {
    norm_samp <- rnorm(nsamp, mean = mu, sd = sigma)
    mu_hat <- mean(norm_samp)
    sigma_hat <- sd(norm_samp)
    
    power.t.test(delta = mu_hat, sd = sigma_hat, n = nsamp, 
                 alternative = "one.sided",
                 type = "one.sample")$power
})

hist(res, 50, main = 'Post-hoc power estimates', xlab = 'Power')
abline(v = real_power, lty = 2)