Power Analysis by Data Simulation in R - Part I: Theoretical introduction to simulation

Why this blog?

In recent years, power-analysis has become a standard tool in the behavioral sciences. With an ongoing replication crisis, high-power research is the a key to improving replicability and to improving the reliability of findings. Especially with preregistration slowly becoming standard practice in psychology, power-analysis, the practice of estimating the required sample-size a-priori, is a more important step than ever to plan research projects accordingly.

Not so long ago, power-analysis was a rather complicated endeavor urging people to use algebraic solutions to calculate power of planned studies which can be demanding especially for non-trivial designs (i.e. basically anything that is not a correlation). Luckily, for many research designs power-analysis is nowadays readily available in software packages such as G*Power and even for relatively complex designs in specialized tools such as the great PANGEA tool for all kinds of generalized ANOVA designs and other tools by Jake Westfall.

However, while these tools are really great they also have (in my personal opinion) two drawbacks. First, they urge the user to familiarize themselves with a new piece of software with new user interfaces, that are not always intuitive. Second, and more importantly, these interfaces promote a statistical way of thinking that often leaves the user confused with what the ever-changing parameters (think of these \(\delta\) , \(d\), \(f\), \(f^2\) and vague “group-size” and “measurement-point” fields) that need to be filled in mean. Moreover, these parameters differ for most designs and give an impression that power-analysis is complicated business better left to statisticians or that it might not be worth the effort.

However, this impression changed dramatically for me, once I changed my ways and started doing power analysis by simulation. Moreover, learning how to simulate data can of course not only be of use for power analysis but is a useful skill to in every research project. When we simulate data, we can see whether what we think about the data-generating process will actually result in the patterns that we would expect. In other words, we can do theoretical experiments to see whether if everything goes as we would expect, we would also find the results that we would expect. Especially when our analyses become more complex, all the different parts in our model tend to interact and it is easy to get lost. By simulating data before conducting an experiment, we will be forced (and able) to specify more precisely how we think a theoretical model or prediction will be reflected in the data.

Thus, even though this tutorial will focus on simulating data for power analysis, you will learn a very useful skill on the side - simulating your own data and thereby conducting theoretical experiments before you even collect data. I know that there are already some excellent tutorials on data/power-simulation out there but they are often very brief and/or technical and assume a rather high level of prior knowledge about R and data-simulation in general. Therefore, I will spend some time on explaining theoretical concepts and slowly build up the simulation-code to hopefully enable the reader to understand the underlying principles and flexibly conduct power simulations themselves after reading this.

Throughout the tutorial I assume readers are familiar with R and some of it’s base functionality. This tutorial will consist of four different parts.

• In the first part (the one you are reading) I will give a short overview of how power-analysis by simulation works on a conceptual level and why I prefer it to available power-analysis software, even though there are clear drawbacks that I will also briefly mention. Moreover, I will introduce the concept of power as I want to bring us all on the same page and want to give a (hopefully) intuitive example about what we actually do in power calculation and how it relates to the simulation techniques that we will use for the rest of the tutorial.
• In the second part, I will discuss simulations for the simplest case of paired and two-sample t-tests.
• In the third part, we will explore different ANOVA and regression designs.
• In the fourth part we will move on to more complex mixed-effects and hierarchical models and even have a peak at Bayesian approaches to power analysis (well technically its not a power analysis but a true detection rate analysis). For this fourth part, I assume readers will be familiar with how to fit mixed-effect models in lmer and/or brms.

This part of the tutorial is by far the most wordy and longest part. Just as a little motivator to keep in mind during this sometimes lengthy tutorial: At the end of part IV of this tutorial, you will know how to do your own custom power analysis for mixed-effects models.

Power Analysis by Simulation (a sustainable alternative)

# first we write a function that calculates the probability for each number so we can call it in a loop
pbinom2 <- function(N, x, p){
  factorial(N)/(factorial(x)*factorial(N-x))*p^x*(1-p)^(N-x)
}

tosses_55to100 <- c(55:100)  # we define the amount of heads that we want to check for (all bigger than or equal to 55)
probs <- c() # we will make an empty collection that we will add the results for each number of heads to

for(i in tosses_55to100){
  probs <- append(probs, pbinom2(100,i,.5))
}
print(probs[1:10])

##  [1] 0.048474297 0.038952560 0.030068643 0.022292270 0.015869073 0.010843867
##  [7] 0.007110732 0.004472880 0.002697928 0.001559739

n_heads <- 54:99  # all possible amount of heads that we want to try.

p_heads <- pbinom(q = n_heads, size = 100, prob = .5, lower.tail = F)  # get pbinom to show us the probability of so many heads for each of the values if a coin is fair

plot(n_heads, p_heads)
abline(h = .01)
abline(v = n_heads[p_heads < .01][1])

Back to Power

In this coin-toss example, we run a statistical test about the fairness of the coin. But what is the power of the test here (in case you forgot during my very long “brief” introduction, this post was about power)?

Actually, we already used three concepts from Frequentist statistical testing:

• The null-hypothesis: we assumed that the coin is fair
• The p-value : the probability that we calculated above when checking how likely 62 or more heads are (we got 0.006).
• The alpha level : (i.e. where people conventionally use .05 in psychological literature⁵) is the threshold that we picked for concluding it would surprise us enough to say that the coin must be unfair, i.e. the probability of .01. In other words, it is the chance that we conclude that a coin is unfair even if it is actually fair (the one percent of cases where 62 or more heads would happen even with a fair coin).

In a proper power analysis, all of these need to be specified in addition to an alternative hypothesis. Let’s give it a try and specify them to run a power-analysis

Specifying the desired power

Moreover, you also want to be sure that you would detect the unfairness of 55% if it is actually there. Let’s say you only want to have a 10 percent chance of not detecting it if it was there. 100 minus this chance is the power of our coin-toss study, i.e. 90 percent (or 1-.10 on the probability scale).

To summarize, our test has the following properties:

• alpha-level = .001
• alternative hypothesis (aka effect size / fairness-criterion) = 55% heads
• power = .90

Our job now is to figure out at which number of tosses we can be 90% sure to detect the unfairness of 55% with only a 0.1% chance of getting into a fight with Harvey by wrongly accusing him of unfairness.

To do this we can use the following r-code (I will explain below)

The first Power Calculation

power_at_n <- c(0) # initialize vector that stores power for each number of tosses
n_heads <- c() # save "critical" number of heads for that toss-amount that would result 
n_toss <- 2 # initialize the toss-counter
while(power_at_n[n_toss-1] < .90){ # continue as long as power is not 90%
  n_heads[n_toss] <- qbinom(.001, n_toss, .5, lower.tail = F) # retrieve critical value
  power_at_n[n_toss] <- pbinom(n_heads[n_toss], n_toss, .55, lower.tail = F) # calculate power (1-beta) for each coin-toss
  n_toss <- n_toss+1 # increase toss-number 
}

The above loop needs some explanation. It increases the toss-amount n_toss by 1 as long as it has not yet reached 90% power. To do this, we again use the qbinom function to find the number of heads (or more heads) that would only occur with a probability of .001. In other words, only 0.1% of the possible coin toss sequences would result in that specific amount of heads when tossing the coin n_toss times. This is the same thing we did earlier with a fixed amount of 100 tosses.

in the next line, we take this amount of heads that only occurs with a probabiltiy of .001 for the current toss amount n_toss and use pbinom to calculate the probability of getting at least this many heads with our hypothesized unfair coin that produces heads in 55% of the cases. This means, we calculate the percentage of coin toss sequences that contain at least this many heads. This probability is the power of the test. Why? Because if 90% of the sequences contain 55% or more heads, if we take the coin and toss it n_toss times,in 90% of the cases we will get one of those sequences that contain 55% or more heads.

Lets have a look at two of the values from this calculation to make this more clear. For instance, lets look at the values when the loop tried out 100 tosses:

n_heads[100] = 65. This number is higher than the 62 heads we had above, as we are more strict now by specifying a stricter alpha level (.001 instead of .01). Next, this 65 was passed on to the pbinom function and we can look at the power, power_at_n[100] = 0.02. At this point we only have a power of .02 or 2%. This means that tossing an unfair coin that would give 55% heads 100 times, only 2% of the sequences would contain 55 or more heads, therfore making it unlikely that we would detect the unfairness. When we increase the number of tosses until the loop stops, we are at n_toss-1 = 1908 coin tosses. At this number of tosses, the amount of heads that would make us conclude that a fair coin is unfair (biased with 55% heads) with only 0.1% chance of being wrong is n_heads-[n_toss-1] = 1021. Thus in this case, with 1908 tosses, if we get 1021 heads or more, we conclude that the coin is unfair. What is the chance of getting at least that with our unfair coin? That’s what the pbinom function in the loop above tells us and it is power_at_n[n_toss-1] = 0.9, our specified 90%. We can also plot the power for each number of tosses that we tried in the loop. Figure 4 shows the increase in power with increasing sample-size.

plot(1:(n_toss-1), power_at_n, xlab = "Number of coin-tosses", ylab = "Power", axes = FALSE)
abline(h = .90, col = "red")
axis(side = 1, at = seq(0,(n_toss-1),by=100))
axis(side = 2, at = seq(0,1,by=0.1))

Thus, when tossing a coin 1908 times¹⁰, in 90% of the cases we will be able to tell it’s biased and can confidently confront Harvey like shown in Figure 5.

If you are not familiar with using these distribution functions (qbinom, pbinom etc.) in R, this might have been a lot of new information but this is basically what we do in power analysis:

1. We specify a null-hypothesis, an alternative hypothesis, an alpha-level and a desired power,
2. We try a small sample-size.
3. We retrieve the critical value, (the number of heads that would be more surprising than what we specified as our alpha-level).
4. We calculate the probability that the amount of heads that would make us reject the null-hypothesis would be observed with the unfair coin according to the alternative hypothesis.
5. We stop as soon as this probability is equal to the desired power.

Finally, an actual power simulation.

So far, we have not done any simulation but have merely analytically derived the power by making use of the binomial probability mass function¹¹. Thus, for easy toy-examples like this one we would not need to do a simulation. However, as soon as we deal with real examples, it is much more difficult to make use of this approach and if we have several predictors in our model, or if we deal with mixed-effect models or hierarchical models (as we will do in part IV of this tutorial) the above method is not feasible anymore.

What we can do however, for any model of any complexity and form, is to actually pretend we were repeatedly doing the experiment for each sample size and see how often we would be able to reject the null-hypothesis. For instance we could toss a coin 20 times and test whether we would reject the null-hypothesis. We could then repeat this process for 20 tosses very often, e.g. 1,000 times and see what the probability is that we would conclude that the null-hypothesis is false¹². This is what we do in power-simulation.

An obvious disadvantage is that instead of just calculating the power for each toss-amount (i.e. sample-size) only once, we need to try each toss-amount out 1,000 times. Thus, simulation takes much longer than a regular power-calculation, especially with more complex models and high sample-sizes. However, the advantage of the method is that we can just learn it once and adjust it for any situation that we will ever find ourselves in, not having to ever walk through tedious interfaces again, selecting arbitrary analyses and setting ever changing parameters to certain values. Another advantage (that will be discussed in detail later) is that we do not need to specify a precise alternative hypothesis and test it for that single value, but that we can actually remain more vague about what our alternative hypothesis (i.e. the effect size that we expect) will be. Oftentimes we do not know exactly what effect-size we can expect and we might like to tell the power-analysis about this uncertainty.

At last, let us do a power-simulation for the above example. Luckily we do not really have to toss a coin as R can do that for us by using the rbinom function, that will as often as we call it do a coin-tossing experiment for us with a specified sample-size. Lets first see how the rbinom function works if we would want to toss a coin 20 times.

set.seed(1) # make sure our simulation will give the same results if you try it
rbinom(n = 1, size = 20, prob = .50) # let r do 1 experiment with 20 coin tosses of a fair coin

## [1] 9

In the above code, R tossed a coin 20 times and it resulted in 9 heads. We could repeat this experiment again:

set.seed(2) # make sure our simulation will give other results than before
rbinom(n = 1, size = 20, prob = .5) # run the experiment again

## [1] 8

In this case, giving us 8 heads. By increasing the first argument to the rbinom function, we can tell R to repeat this experiment more often. Moreover, we can tell it to make use of an unfair coin directly, so we can directly put our alternative hypothesis in the simulation by changing the last number of the rbinom function from .50 to .55, to do the same test as above.

set.seed(1)
n_heads <- rbinom(n = 1000, size = 20, prob = .55) # run 1,000 experiments, of 20 coin tosses each, at once
str(n_heads) # show structure of vector

##  int [1:1000] 12 12 11 8 13 8 7 10 10 14 ...

Now, R repeated the 20 coin-toss experiment 1,000 times with an unfair coin of 55% chance of resulting in heads giving us 1,000 times the amount of heads that it got. Let us again now test how big our power was in this case, again with an alpha-level of .001.

p_heads <- pbinom(n_heads, 20, .50, lower.tail = F) # calculate the probability of observing this many heads if the coin would be fair (which it is not cause we simulated with 55% heads)

exp_power <- mean(p_heads < .001) # check where this chance drops below our alpha level

The first line in the above code does exactly what we did earlier, just with a little change. We take the amounts of heads that we got from an unfair coin, and check in how many cases we would conclude - assuming that the coin would actually be fair, thus using .50 as the probability in pbinom - that the observed amounts of heads is too unlikely for us to believe that the coin was fair. We save these probabilities to a vector. The second line calculates the observed power of our experiment. To understand what it does, let us have a look at the vector p_heads it looks the following (here only the first 10 out of 1,000 values):

p_heads[1:10]

##  [1] 0.13158798 0.13158798 0.25172234 0.74827766 0.05765915 0.74827766
##  [7] 0.86841202 0.41190147 0.41190147 0.02069473

As shown above, the p_heads vector contains the probability of observing each amount of heads n_heads from our 1,000 experiments assuming the experiment was done with a fair coin (which it was not). Now we would like to check how many of these probabilities are at least as small as our alpha-level, i.e. surprising enough to conclude the coin was not fair. We do this by checking for each value whether it was .001 or smaller (p_heads < .001). This will result in another vector of 0 when the condition is false and 1 when the condition is true¹³. Taking the mean of this vector will give us the probability of rejecting the null-hypothesis while we actually know it is incorrect (as we put the bias in the coin ourselves). Again, this will give us the power. In the present case with 20 coin-tosses the power is exp_power = 0.004 or 0.4 percent. This is obviously very low and not surprising given that we already know from the calculation above that we need a lot more coin-tosses than 20 to get the desired power of 90%. To get to our desired power in this example, we need to change the code above so it will try different sample-sizes again.

set.seed(1)
exp_power_at_n <- c(0) # create a vector where we can store the power for each sample-size
n_toss_start <- 19 # start at 21 tosses
n_toss_loop <- 2 # additional number of tosses tried (above 20)
while(exp_power_at_n[n_toss_loop-1] < .90){ # continue increasing the sample-size until power = 90%
  n_toss <- n_toss_start+n_toss_loop # calculate the current number of tosses
  n_heads <- rbinom(1000, n_toss, .55) # run 1000 experiments for any given number of tosses and store number of heads
  p_heads <- pbinom(n_heads, n_toss, .50, lower.tail = F) # calculate the probability of getting at least that many heads if the coin would be fair 
  exp_power_at_n[n_toss_loop] <- mean(p_heads < .001) # calculate power by checking what proportion of the probabilities is smaller than or equal to our alpha-level
  n_toss_loop = n_toss_loop+1
}

exp_power_at_n <- exp_power_at_n[-1] # remove the first 0 that we used to populate the vector for the first iteration of the loop

The above code is similar to what we have done earlier when we tried only 1 sample-size. This time, we iterate over different sample-sizes in a loop and store the power for each in the vector exp_power_at_n. As in the earlier calculation, we can now use this to see how many tosses we would need by having a look at where the loop stopped, i.e. when it reached 90% power, which is at n_toss-1 = 1795 tosses at which the power was exp_power_at_n[length(exp_power_at_n)] = 0.903. We can also plot all these values again as done in Figure 6.

plot(21:n_toss, exp_power_at_n, xlab = "Number of coin-tosses", ylab = "Power", ylim = c(0,1), axes = FALSE)
abline(h = .90, col = "red")
axis(side = 1, at = seq(0,(n_toss),by=100))
axis(side = 2, at = seq(0,1,by=0.1))

The shape of the line looks very similar to the earlier calculation, however the line appears to be thicker. Moreover, it might be surprising that we did not get the same amount of tosses that we got from the calculation. Actually they differ quite a lot (1908 vs. 1795). This is due to the fact that even when running 1,000 experiments for each sample-size there is still imprecision in the simulation. Each coin-flip is random and even if we repeat a experiment 1,000 times this randomness is still in there (think of how many possible sequences there would be!). This randomness is why we cat a different number in the simulation and why the line in Figure 6 is thicker than in Figure 4.

If we want to approach the results of the calculation more closely in our simulation (i.e. get a more precise power-estimate) we can increase the number of experiments that R will run per sample-size (i.e. the number of simulations). For example, we could repeat the simulation with 100,000 experiments per sample-size. You have to be patient here, this already takes a few minutes maybe.

set.seed(1)
exp_power_at_n <- c(0) # create a vector where we can store the power for each sample-size
n_toss_start <- 19 # start at 21 tosses
n_toss_loop <- 2 # additional number of tosses tried (above 20)
while(exp_power_at_n[n_toss_loop-1] < .90){ # continue increasing the sample-size until power = 90%
  n_toss <- n_toss_start+n_toss_loop # calculate the current number of tosses
  n_heads <- rbinom(100000, n_toss, .55) # run 1000 experiments for any given number of tosses and store number of heads
  p_heads <- pbinom(n_heads, n_toss, .50, lower.tail = F) # calculate the probability of getting at least that many heads if the coin would be fair 
  exp_power_at_n[n_toss_loop] <- mean(p_heads < .001) # calculate power by checking what proportion of the probabilities is smaller than or equal to our alpha-level
  n_toss_loop = n_toss_loop+1
}

exp_power_at_n <- exp_power_at_n[-1] # remove the first 0 that we used to populate the vector for the first iteration of the loop

If you run this code, you will see that it took substantially longer than the previous simulation with only 1,000 repetitions. If we now look at the outcome again we find that the number of tosses n_toss-1 = 1870 is already slightly closer to the calculated value but still not exactly the same. In my personal opinion, this imprecision is both advantage of simulation as well as disadvantage. It is a disadvantage in that it is less precise than the calculation above in theory. However, it is an advantage as it adds some noise to the power-estimation process that is actually also present in real life. This is, even if there is an effect in the population, each new sample will always be different which is also the case in the simulation. If we look at Figure 7, we also see that now the power-increase follows a more straight line again that is very close to the one from the exact calculation and not as thick anymore as in Figure 6.

plot(21:n_toss, exp_power_at_n, xlab = "Number of coin-tosses", ylab = "Power", ylim = c(0,1), axes = FALSE)
abline(h = .90, col = "red")
axis(side = 1, at = seq(0,(n_toss),by=100))
axis(side = 2, at = seq(0,1,by=0.1))

Summary

In this part of the tutorial I tried to bring us all on the same page about why we would do data or power simulation and how it is related to exact power calculation.

• In principle, in power calculation we count the amount of potential data sequences (e.g. possible sequences of HEADS-TAILS in coin-tosses) that could happen in an experiment of a specific size (e.g. amount of coin-tosses).
• Afterwards, we calculate the critical value (e.g. number of heads) that would make us reject the null-hypothesis (e.g. a coin is fair) as only a very small percentage of sequences (“small” is defined by the alpha-level here) would surpass the critical value.
• Finally, we calculate the percentage of sequences according to the alternative hypothesis (e.g. a coin is unfair with a bias of 55%) that would surpass the critical value. This is the power of our test.

In power-simulation we do something very similar, but instead of calculating how many possible data sequences there are and how many of them result in a critical value, we just try this out a lot of times until we get an approximation of how often we would observe data that are inconsistent with the null-hypothesis given our specified alpha-level.

I hope you forgive my rather lengthy introduction on power-analysis before actually doing the simulation, as I thought it would be a good foundation before we move on to more complex and realistic situations.

• In part II of this tutorial, we will move on to these more realistic situations that we might actually be interested in psychology by looking at how we do simulations for t-tests.
• In part III we will continue with ANOVA and regression designs.
• In part IV we will look at mixed-effects/multilevel models.

Footnotes