power_analysis.R

library(ggplot2)

# ------------------------------

gen_dat = function(baseline, change, sd) {
  data.frame(
    y = c(rnorm(n=320, mean=baseline, sd=sd),
      rnorm(n=90, mean=baseline+change, sd=sd)),
    x = c(rep("out", 320), rep("in", 90))
  )
}

get_pvalue = function(data) {
  summary(lm(y ~ x, data))$coef[2, 4]
}

# ------------------------------

set.seed(18959871)
baseline = 4
change = seq(0.1, 1, 0.1)
sd = c(2.1, 2.7)

results = expand.grid(baseline = baseline, change = change, sd = sd)
results$power = NA

for (i in seq_len(nrow(results))) {
  ps = sapply(1:1000, function(x) {
    get_pvalue(gen_dat(results$baseline[i], results$change[i], sd = results$sd[i]))
  })

  results$power[i] <- mean(ps < 0.05)
  print(i)

}

results$sd = factor(results$sd)

ggplot(results, aes(x = change, y = power, color = sd)) +
  geom_line() + theme_bw() +
  geom_hline(yintercept = 0.8, color = "red") +
  labs(title = "Statistical power by effect size and SD",
    subtitle = "Assuming normal distribution, mean = 4",
    y = "Power", x = "Effect size")