retrospective_powerR

library("pwr")
library("ggplot2")

n_sims <- 1000L
n <- 10L
mu <- 0.5
sigma <- 1

true_power <- power.t.test(
    n = n,
    delta = mu,
    power = NULL,
    sd = sigma,
    type = "one.sample",
    strict = TRUE
)$power

powers <- rep(0.0, n_sims)
p_values <- rep(0.0, n_sims)

for (s in 1L:n_sims) {
    x <- rnorm(n, mu, sigma)

    p_values[s] <- t.test(x)$p.value

    powers[s] <- power.t.test(
        n = n,
        delta = mean(x),
        power = NULL,
        sd = sd(x),
        type = "one.sample",
        strict = TRUE
    )$power
}

df <- data.frame(power = powers, p_value = p_values)

ggplot(df, aes(x = power, y = p_value)) +
    geom_point() +
    geom_vline(xintercept = true_power)

ggplot(df, aes(x = power)) +
    geom_density() +
    geom_vline(xintercept = true_power)

mean(df$power) / true_power
mean(df$power[df$p_value < 0.05]) / true_power