POWER CALCULATIONS AND SIMULATIONS

John Tebes and Rohit Naimpally

Thank you to Chris Boyer and Erica Chuang for stress-testing this!

20 September 2016

Why are power calculations important?

Power influences many design aspects, including what research questions to pursue, how many treatment arms to employ, and even more fundamentally, whether or not to proceed with a potential research project. For example, it may be that a remedial education program boosts tests scores by 20 percent when comparing treatment and control groups, but due to limited power, the RCT is unable to detect this true effect (with 95% confidence). However, we can estimate whether a given design is likely to be able to detect a reasonable effect size ex ante, allowing us to properly manage partner expectations and make the most of limited research resources.

This exercise will cover two ways of calculating power for a RCT - (1) the conventional parametric method, and (2) a non-parametric "simulation" method.

This exercise is written in SMCL (pronounced "smickle"), Stata's output language. Your working directory needs to be set correctly using cd. Now, set your working directory to the directory that contains this SMCL file.
If you opened this file by double-clicking on it, you probably don't need to change your working directory.

The fact that we're using SMCL means that we can run commands from within the file. For example:

`"display "Hello world!""'

Commands that are one line long are blue and clickable.

Further, you will see links like this one that take you to help files, locations in a Stata manual, or even a webpage. These links work like those in an Internet browser,
so you can always press the Back button to return to your exact previous location. Further, if you right-click on a link, you have the option to open it in a new tab or window.

Questions to consider before running power calculations...

- What is the main specification (e.g. regression) we plan to run? [It doesn't have to be fully baked, but the more "baked" it is, the more precise we can make your power estimates.]

- What do we expect to be the mean of the outcome in the control group?

- How about the standard deviation (SD) of the outcome in control group? [If your outcome is binary, we can estimate this without baseline data by assuming this Bernoulli random variable takes a binomial distribution.]

- What sample sizes are feasible?

- What effect sizes could the intervention reasonably cause?

- What is the smallest, cost-effective effect size that we are interested in? [We often arrive at a reasonable answer to this question through discussions with partner organizations and literature reviews.]

Let's start with a simple parametric example.

EXAMPLE 1: Basic Parametric Example

+
+

| IMPORTANT NOTES! |

+
+

| To run the Stata sections, we will need Stata (version 13 preferably) and |

| Microsoft Excel. If you have not already, install the Stata commands that |

| are in the ado folder; we will need them for this exercise. |

| |

| This section creates globals that will be called upon in later sections, |

| so it's important to run all the Stata code in this section before |

| proceeding. |

| |

+
+

clear all

set more off

version 11

Let's load the Balsakhi dataset. We'll use this dataset to estimate the control
mean.

use "baroda_0102_1obs.dta",clear

For all parametric power calculations, we'll assume a conventional
95% confidence interval and 80% power.

What do we expect the mean of the outcome to be in the control group?

Note: Since power calculations are usually done prior to a study, we often use
baseline/pilot data on the study population, or government statistics for a comparable
population, to get an approximate for this outcome in the control group.

sum post_totnorm if bal==0

global sd = r(sd)

global control = r(mean)

Let's say, based on other studies, that we expect an effect size of a tenth of a standard deviation. Now let's calculate the sample size given that we know the likely effect size.

global effect = $sd/10

global treat = $control + $effect

sampsi $control $treat, power(.8) sd($sd)

Say, instead, we knew the sample size and wanted to calculate the Minimum
Detectable Effect Size (MDE). Although we cannot use the sampsi command to back out
the effect size, the command power, available in Stata versions
13 and higher, allows for this:

global n = 3142

power twomeans $control, power(.8) n($n) sd($sd) table

You can also manually back out the effect size using the formula on slide 41 here:

di (0.842+1.96)*sqrt(1/0.25)*sqrt(1/$n)*$sd

IMPORTANT NOTE: The effect size should be .1151 (some Windows computers have had glitches when running this code). Verify you get this result before proceeding. Let's compare this to the likely effect size that went into the first sample size calculation:

di $effect

As this shows, it doesn't matter which we start with - sample size or effect size.

Some other questions to answer before calculating power:

- Will this study be cluster-randomized?

- Will our main specification include controls (i.e. lagged dependent variables)?

- Do we expect only part of the treatment group to take-up the
intervention? If so, are we interested in estimating the local average treatment effect?

We'll address each of these questions one at a time to see how they affect power.

Table of contents

2. Building Intuition

3. Parametric Power Calculation with Controls

4. Parametric Power Calculation for Cluster-RCTs

5. Parametric Power Calculation with Partial Take-up
6. Non-Parametric Power Simulation

7. Answer Key

EXAMPLE 2: Building Intuition

Now, let's get a better intuition on how a larger or smaller sample size affects our power to pick up an effect.

Say our anticipated effect size is smaller than originally thought; how much larger would we need to make the sample in order to still pick up an effect?

Let's try an effect size that is half as large.

global new_effect = $effect/2

global new_treat = $control + $new_effect

sampsi $control $new_treat, power(.8) sd($sd)

Observation: The minimum sample required is four times as large.

Let's try an effect size that is a third of the original.

global new_effect2 = $effect/3

global new_treat2 = $control + $new_effect2

sampsi $control $new_treat2, power(.8) sd($sd)

Observation: The minimum sample required is now nine times as large.

Remember: If our MDE decreases by a factor of X, the required sample size increases by the square of X!

If you have Stata 13 or higher, you can verify this from the other side i.e. look at the impact on the MDE of increasing your sample size by a factor of X. Say X is 4:

global new_samp = $n*4

global new_effect3 = (0.842+1.96)*sqrt(1/0.25)*sqrt(1/$new_samp)*$sd

di $new_effect3/$effect

EXAMPLE 3: Parametric Power Calculation with Controls

Now, say we plan to control for baseline covariates in our main specification.
The inclusion of these controls will improve our power, since they explain some of the
variance in our outcome. For example, including prior test scores on the right-hand side as controls when our left-hand side variable is test scores during the study period, can
improve power if prior test scores predict future test scores. Economists will often control for a lagged dependent variable (recorded prior to randomization) if
it makes a first order improvement in power.

To see how potential controls affect power, we would ideally have access to a sample data set (e.g. historical or pilot data).
With these data, we would want to regress Y_i (the outcome) on X_i (the controls) to evaluate how much variance is
explained by the set of covariates we plan to include. With access to historical data, for example, this would involve regressing last year's test scores (Y_i = Y_t-1) on
test scores from the year before (X_i = Y_t-2).

From this regression, we are interested in the residual standard deviation of
the outcome variables, or the variance of the outcome that is NOT explained
by controls. This residual SD becomes the new SD we include in our parametric power
calculations.

Part 1. We have pilot/historical data.

Using Balsakhi data, this would be:

`"reg post_totnorm pre_totnorm if bal==0"'

`"predict res, res"'

`"sum res"'

`"global res_sd = r(sd)"'

If we knew the effect size and wanted to know the sample size needed...

`"sampsi $control $treat, power(.8) sd($res_sd)"'

If we knew the sample size was 3,142 and wanted to calculate the effect size (in Stata 13 or higher)...

`"power twomeans $control, power(.8) n($n) sd($res_sd) table"'

Questions:

1. What percent of the variance of study period test scores is explained by test scores at baseline? (Hint: Look at the R^2 statistic of the regression.)

2. How does this affect our sample size (compared to not including controls)?

3. How about our MDE?

Part 2. We do not have pilot/historical data for our study population, but plan to include controls.

We must first guess the percentage of variance in outcome variables that we
expect to be explained by controls. (This can be done by looking at other datasets
with similar outcomes, etc.)

Say our controls explain 49% of the outcome. We can then calculate the residual standard
deviation as follows.

`"global res_sd2 = sqrt(.51*$sd^2)"'

Then, we run our power calculations with this residual SD.

`"sampsi $control $treat, power(.8) sd($res_sd2)"'

We can also trace out an MDE "range" using various assumptions for how much of the variance of the outcome is explained by covariates:

`"do controlling_for_covariates.do"'

EXAMPLE 4: Parametric Power Calculation for Cluster-RCTs

Many designs randomize at the group level instead of at the individual level. For such designs, we need to adjust our power calculations so that they incorporate
the fact that individuals within the same group may be subject to similar shocks, and thereby have correlated outcomes. Duflo et al. presents a modified
parametric approach, which takes into account the intra-cluster correlation (ICC) that arises from randomization at the group level.

We can think of cluster-RCTs as follows:

- When ICC = 0, then our N is effectively the number of individuals in the study.

- When ICC = 1, then our N is effectively just the number of clusters.

- Usually the ICC lies somewhere between 0 and 1, requiring that we adjust our power calculations to account for this.

Below we adjust Stata's power estimates based on Duflo et al.'s model.

Note: This model assumes that all clusters are of the same size and have the same number of individuals.
It's usually okay if this is violated in reality, but you would not want to use these adjustments if groups are
dramatically different in size (e.g. group one has 10 individuals, group two has 1,000 individuals.) More on this
model is explained in Duflo et al.'s article "Using Randomization in Development Economics Research: A Toolkit."

Part 1. Calculating MDE

First, let's calculate the intra-cluster correlation (ICC) which measures how
correlated the error terms of individuals in the same cluster are.

`"loneway post_totnorm divid if bal==0"'

`"global rho = r(rho)"'

Now, let's specify the number of individuals in each cluster...

`"global m = 53"'

... and the number of clusters (as documented in the Balsaki experiment).

`"global j = 193"'

Let's assume 95% confidence intervals and 80% power...

`"global t = 1.96 + .842"'

... and 50% of the study population is assigned to treatment and 50% to control.

`"global P = .5"'

Now, we have Duflo et al.'s power adjustment:

`"global mde = $t*sqrt(1/($P*(1-$P)*$j))*sqrt($rho + (1-$rho)/$m)*$sd"'

`"global mde_ldv = $t*sqrt(1/($P*(1-$P)*$j))*sqrt($rho + (1-$rho)/$m)*$res_sd"'

And lastly, the total N of our study and the number who are treated, respectively.

`"global n = $j*$m"'

`"global treated = $n*$P"'

`"matrix mde_clus = (.05,0.8,$rho,$j,$m,$n,$treated, $control, $sd, $mde) \ (.05,0.8,$rho,$j,$m,$n,$treated, $control, $sd, $mde_ldv)"'

`"matrix colnames mde_clus = "Signif" "Power" "ICC" "Clusters" "Cluster_size" "N" "Treated" "Cntrl_mn" "Cntrl_SD" "MDE""'

`"matrix rownames mde_clus = "No_controls" "Control_LDV""'

`"matrix list mde_clus"'

Part 2. Calculating Sample Size

Stata's sampclus command can be used in combination with the sampsi command used above to incorporate clustering
into our sample size calculations. A more parsimonious alternative is the clustersampsi command. However, note that clustersampsi rounds
the group means to the tenth place after the decimal, making for a less precise calculation. We demonstrate the use of both commands below:

`"global treat2 = $control + $mde"'

`"clustersampsi, samplesize mu1($control) mu2($treat2) m(53) rho($rho) beta(0.8)"'

To use sampclus, we first run the sampsi command.

`"sampsi $control $treat2, alpha(.05) power(.8) sd($sd)"'

Then, we adjust for the number of individuals per cluster.

`"sampclus, obsclus(53) rho($rho)"'

If we include the lagged dependent variable as a control we have...

`"global treat_ldv =$control + $mde_ldv"'

`"sampsi $control $treat_ldv, alpha(.05) power(.8) sd($res_sd)"'

`"sampclus, obsclus(53) rho($rho)"'

Note that again it doesn't matter if we start with the MDE, or with the sample size.

Questions:

1. Why do we have to adjust power for clustering when running a cluster-RCT?

2. Assuming ICC>0, does adding a new cluster of 5 individuals or adding 5
individuals to already-existing clusters give us more power to detect effects?

EXAMPLE 5: Parametric Power Calculation with Partial Take-up

In randomized designs, it is common that there is partial take-up of the intervention.
For example, in the Oregon Health Insurance Experiment, the offer to apply for
health insurance was associated with only a 25 percentage point increase in
take-up of health insurance. When take-up is not 100 percent, researchers are often interested in
the answers to second stage questions, such as what the average effect of
becoming insured is on health care utilization.

Below, we provide code that adjusts Stata's power to take into account that only
some individuals in the treatment group take up the intervention.

The measure of take-up that we care about is "effective take-up", or the
percentage of individuals in the treatment group that takes up the intervention
MINUS the percentage of individuals in the control group that takes up.

Let's say 90% take up in the treatment group, and 10% do so in the control group. We then have an effective take-up rate of 80%:

`"global eff_tu = .9 - .1"'

Now we calculate the adjusted MDE by multiplying the unadjusted effect size times the effective take-up rate...

`"global adjusted_mde=((0.842+1.96)*sqrt(1/0.25)*sqrt(1/$n)*$sd)*$eff_tu"'

`"di "Adjusted MDE = $adjusted_mde""'

If we want to estimate sample size, we follow a similar procedure:

`"global treat_adjusted = $control + $adjusted_mde"'

`"sampsi $control $treat_adjusted, power(.8) sd($sd)"'

Note that unlike with the MDE, partial take-up affects sample size quadratically; we divide estimates of sample size by the square of the effective take-up rate!

Note: For more on partial take-up and how it affects estimation of the power to
detect local average treatment effects, please see the aforementioned article by Duflo et al.

EXAMPLE 6: Non-parametric Power Simulations

Non-parametric power simulations do better than parametric power calculations when
we have access to good data (historical, baseline, or pilot) on our study
population. From these data, we can simulate a fake dataset that assumes the
treatment has no effect and then see what effects we are powered to detect, by
looking at the simulated 95% confidence interval around our null effect. We should expect that any effect greater than this confidence interval would be detected by our study.

In particular, power simulations do not require the assumption that the
sampling distribution of your Beta coefficient(s) of interest takes a normal
distribution in your (finite) sample. You may be particularly worried about this
assumption (of parametric power calculations) if your sample is very small.

To do non-parametric power simulations we need to create a (reasonable) "fake" or simulated dataset. For example, if you have baseline data for the 3 months prior to a 12 month trial, then
a reasonable way to expand this dataset would be to simply randomly draw days with replacement until you have 365 days in your dataset. Similarly, you could use historical data
from the two years prior to the study to estimate the confidence interval around a null effect in the past year, with data on your outcome variable from two years ago serving as controls.

Let's do an example using Balsakhi data to make this procedure more clear.

Step 1. Upload pre-period data

For this example, we expand study data for the control group of the Balsakhi
experiment by 2, since this is similar to what would be available from historical data.

`"program create_data"'

`"use baroda_0102_1obs,clear"'

`"keep if bal ==0"'

`"expand 2"'

`"bys studentid: replace studentid = studentid*-1 if _n==2"'

`"keep divid studentid *_tot*"'

`"save bal_power_data, replace"'

`"end"'

`"create_data"'

Step 2. Write randomization code as you plan to randomize.

Your randomization code should take into account how you plan to randomize,
including your plan to stratify, etc. For simplicity, we just randomize by
grade-year in this example and do not stratify (although the actual RCT did stratify).

`"program randomization"'

`"preserve"'

`"bys divid: keep if _n==1"'

`"gen treatment = (runiform()>=.5)"'

`"keep divid treatment"'

`"save rand, replace"'

`"restore"'

`"end"'

Step 3. Write simulation code.

This power simulation program does the following:

(a) specifies main regression equation,

(b) simulates treatment assignment and runs main regressions 1000 times,

(c) summarizes the results from these regressions into an easy-to-read matrix.

To reduce the run-time, we have set this program to run over 100 iterations, but if
you are truly calculating power via simulation, you should set this at 1000.

`"do simulation_program.do"'

Part 1. Basic MDE (clustering on divid because data is from C-RCT)

`"simulation_program, reps(100) name(mde_cluster) outcome(post_totnorm) regression("reg post_totnorm treatment, cluster(divid)")"'

Check out the results recorded in the excel file named "JPAL Power Exercises". We can interpret the high (or low) side of the confidence interval as the effect size we are powered to detect.
In this example, you should find that we can detect an ~.20 standard deviation increase (or decrease) in test scores.

Note: A check on whether the program is working properly, is to see if the
p-value is ~.05, since we assume a 95% confidence interval and no effect.

Part 2. Controlling for LDV

`"simulation_program, reps(100) name(mde_cluster_LDV) outcome(post_totnorm) regression("reg post_totnorm treatment pre_totnorm, cluster(divid)")"'

Since the Balsakhi dataset has many observations, we find that our simulated power estimates of MDE are very similar to those found through parametric estimates. To see this, compare Example 4 (parts 1 & 2) to Example 6 (parts 1 & 2).

As you can see, it is much easier to cluster or control for covariates using
this method; you simply run the regression specification that you intend to use at the analysis stage.

Similarly, if you wanted to estimate the (local) average treatment effect in the presence of partial take-up of the intervention, then you can simply make basic
assumptions on the effective take-up rate, and run your two-stage least squares regression in Stata.

Generally speaking, parametric power calculations are great for back-of-the-envelope
calculations of power, while non-parametric simulations provide more precise estimates
of power when you have good baseline data.

Questions:

1. How is "non-parametric" power different from "parametric" power?

2. How do our parametric and non-parametric estimates of power compare? (Hint: Compare Example 4 to Example 6.)

3. When would you want to run parametric power calculations?

4. Non-parametric simulations?

7. Answer Key

Ex. 3, part 1:

1. 49%

2. Including controls cuts sample size almost in half.

3. Reduces MDE to a lesser extent (~28%).

Ex. 4, part 2:

1. Assignment is only random at the cluster level; thus we must cluster our standard errors in our main specification. To this end, our power calculations must also take this into account,
since, by clustering in our main specification, we will lose all precision gained from intra-cluster correlation in outcomes.

2. Adding a new cluster of 5 individuals, since in the limiting case of ICC=1, adding 5 individuals to previous clusters would have no effect on power, while adding 5 individuals in
a new cluster would increase our effective N by 1.

Ex. 6:

1. Non-parametric power simulations do not require the assumption that beta coefficients take a normal distribution (which follows from the central limit theorem), while parametric power
calculations do.

2. Our power simulations in Example 6 find that we are powered to detect an ~.20 SD increase in test scores without controls and a ~.14 SD increase in test scores when including a lagged
dependent variable. Using parametric methods, Example 4 illustrates these are .19 SD and .13 SD, respectively. Thus, we find very similar results using each method.

3. We would want to run parametric power calculations when we want a quick, rough estimate of power, or when we do not have access to high quality baseline data.

4. We would want to run non-parametric simulations when we have access to high-quality baseline data and have reason to believe parametric assumptions may be violated. For example, we would
prefer to run simulations if our randomization is complex (i.e. multiple randomizations), have a small sample size, or are running complex specifications.