Power analysis

draw_sample <- function(x, n) {
  ii <- sample(x, size = n, replace = TRUE)
  x[ii]
}


calc_p_value <- function(effect_size, n) {
  dd <- tibble(x = rnorm(n), 
               y = rnorm(n, effect_size)) |> 
    pivot_longer(everything())
  t.test(value ~ name, data = dd)$p.value
}

calc_power <- Vectorize(
  function(effect_size, n, nr_of_trials) {
    p_values <- replicate(nr_of_trials, calc_p_value(effect_size, n))
    return(mean(p_values < 0.05))
  }
)

calc_power(0.25, 4, 1e3)

[1] 0.053

get_n_needed <- Vectorize(
  function(my_delta, my_power) {
    power.t.test(delta = my_delta, sd = 1, sig.level = 0.05,
                 power = my_power)$n
  }
)

dat <- 
  expand_grid(power = plogis(-3:3),
              effect_size = c(2^(-2:2))) |> 
  # slice(1:5) |> 
  mutate(n_needed_exact = get_n_needed(effect_size, power),
         n_needed = ceiling(n_needed_exact),
         power_simulated = calc_power(effect_size, n_needed, 1e3)) |> 
  pivot_longer(cols = c(power, power_simulated))

dat |> 
  ggplot(aes(n_needed_exact, value, colour = name)) + 
  facet_grid(cols = vars(effect_size), 
             scales = 'free', 
             labeller = label_both) + 
  geom_hline(yintercept = 0.8, lty = 'dashed') + 
  geom_line() + 
  scale_x_continuous(trans = 'log10') + 
  ylim(0, 1) +
  labs(colour = 'Effect size')