How many simulations in my power analysis?

Today I was trying to figure out how to advise on the number of simulations to run when calculating power by simulation.

I tackled this question by running a simulation (of course).

library(ggplot2)
library(dplyr) # I love pipes
set.seed(8675309)

I wanted to figure out how close to the true power was the calculated power from a simulation where the number of replications ranges from 100 to 10K (in steps of 100) and power ranges from 0.5 to 1 in steps of .05 (the result is symmetric around 50%, so the figures below for 80% power also apply to 20% power)..

First, I made all possible combinations of replications and power.

x <- expand.grid(
  reps = seq(100, 1e4, 100),
  power = seq(0.5, 1, .05)
)

Then, for each combination, I calculated the proportion of significant analyses in 10K simulations. I assumed this would have a binomial distribution where size is the number of replications in each simulation and probability is the true power. I then calculated the absolute difference from the true value of power and reported the mean (I find it more intuitive than SD or variance).

x$diff <- mapply(function(size, prob) {
  sig <- rbinom(1e4, size, prob) / size
  diff <- abs(sig - prob)
  mean(diff)
}, size = x$reps, prob = x$power)

I plotted the results to see if they make sense. As the number of replications per simulation increases, the mean difference from the true power decreases. Accuracy is higher for larger values of power.

I also calculated the minimum number of replications to get a result that is, on average, less than 1% off from a power of 80%

filter(x, power == .8, diff < .01) %>% 
  pull(reps) %>%  min()

## [1] 1100

I also calculated the .95 quantile to see how many replications you need to run to get within 1% of the true value.

x$q95 <- mapply(function(size, prob) {
  sig <- rbinom(1e4, size, prob) / size
  diff <- abs(sig - prob)
  quantile(diff, .95)
}, size = x$reps, prob = x$power)

Turns out you need a lot more.

filter(x, power == .8, q95 < .01) %>% 
  pull(reps) %>%  min()

## [1] 6300