# Power analysis/sample size determination by simulation: Optimal algorithm for finding $N$

I want to calculate the needed sample size $$N$$ to achieve $$80\%$$ power for a complex longitudinal study design. Because of the complexity of the design and because I want to have control over every parameter, I'm using simulations to calculate the power at a particular sample size. I'm using R for the simulations. Unfortunately, even a simulation run with 1000 repetitions takes multiple minutes to finish. Because of that, I want to keep the number of simulation runs as low as possible.

If there was an analytical solution (a formula), I could use a root finding algorithm to find the $$N$$ that achieves $$80\%$$ power. Using simulations, however, the estimates of the power at a given sample size are noisy. The monte-carlo standard error of the simulation based power is simply $$\sqrt{\hat{p}(1 - \hat{p})/n}$$, where $$n$$ is the number of simulations and $$\hat{p}$$ the estimated power (the proportions of simulations were significance was achieved).

Question: What strategy/algorithm could I use to automate the process of finding the optimal $$N$$ with as little simulation runs as possible?

I don't have a specific algorithm for you, but here is some general advice.

First, you may want to look to see if you can speed up the simulation analysis process. There are little inefficiencies in a lot of existing code that generally costs less than a second when running once (so most people will not care), but that time can add up in simulations when running thousands of times. This could be something like the time it takes the program to convert categorical predictors into indicator variables, check for missing data, calculate transformations, calculate confidence intervals, tests on parameters that you are not interested in, etc. If you can move these steps outside of the analysis and just do them once up front (in the data generation phase, then skip the extra steps in the analysis) this could speed up your simulations.

Look at how you are doing your simulations, is there a way to be more efficient there? For example growing a data structure in a for loop in R is much slower than using functions like lapply, sapply, or replicate (or the functions in the purrr package). Using macro loops in SAS will take orders of magnitude longer than generating all the data in one data step and using by processing to analyze.

Most computers today have multiple processors, are you using parallel processing to speed things up? For short running simulations parallel processing can slow things down (the extra overhead), but as the individual simulation pieces take longer parallel processing can be a real time saver. Parallel processing can be as simple as having multiple copies of the program running on different sample sizes at the same time (this may depend on what software you are using).

Which part of the individual simulations is taking the most time? creating the data? or analyzing the data? If creating the data is taking the most time, you can look at multiple sample sizes in one simulation by generating the largest dataset, then using subsets of that for the other sample sizes, e.g. create data for 1,000 subjects over your time periods/events in the longitudinal design, then analyze using the first 200 subjects, the first 400 subjects, the first 600, the first 800, and finally the full 1,000. You only had to generate the data once, but have information on 5 different sample sizes.

As you point out, there is uncertainty in the estimated power from a simulation and a small number of simulations does not give very much information, but sometimes a small amount of information is enough. If you only run 100 simulations then the margin of error is about 10%, but if the sample size you are testing gives an estimate of 40% power then knowing the power could be anywhere between 30% and 50% still tells you that that sample size is too small to give you 80% power without the time needed to run 1,000 or more simulations. So run several sample sizes with only 100 simulations to get in the general area of 80% power, then run some with 1,000 simulations to further narrow it down, then possibly a final run with more simulations to verify.

Don't throw away the information from the simulations that do not give you what you want, you can fit a logistic regression using the the sample size as the predictor (using splines or polynomial terms to allow non-linearity) and the power (and number of runs) as the response. Plot the predicted curve along with its confidence interval to suggest what sample sizes are likely to give you your desired power, run a few more simulations with those sample sizes and update your logistic regression plot to suggest next sample sizes to test, etc.

• Thanks Greg, this is excellent advice! I especially like the application of logistic regression. The part that's taking the most time (by far) is the lmer function to fit a linear mixed model. I've already looked into the performance tips but it didn't make a lot of difference. I'm using R (will update my question accordingly). Commented Jan 4, 2021 at 16:52
• @COOLSerdash, If I remember correctly, lmer uses a Bayesian like McMC algorythm behind the scenes to do some of the inference (which is probably taking the most time). There may be some options to use a faster approximation at least for the narrowing down steps, the nlme package may be quicker than lme4 if it can fit your model. Doug Bates (who wrote the original lmer code) has moved a lot of his development work to the Julia language, I don't know if it would be faster to use an R/Julia interface package to run the Julia code. Commented Jan 4, 2021 at 17:05
• Hmmm, I had several tips I thought I could share, but you've already given them (plus others...). Commented Jan 4, 2021 at 17:27