»  Home »  Resources & support »  FAQs »  How to calculate power using simulation

How can I use Monte Carlo simulations to estimate power in Stata?

How can I integrate a simulation program into the power command?

Title   How to calculate power using simulation
Authors Alan Feiveson, NASA
Chuck Huber, StataCorp
Mia (Dan) Lv, StataCorp

Power and sample size analysis is an important tool for planning your experiments. Stata's power command has several methods implemented that allow us to compute power or sample size for tests on means, proportions, variances, regression slopes, case-control analysis, and survival analysis, among others. For those complicated models that are not directly supported by the power suite of commands, we can consider using simulations to obtain the power. In Stata, we can use simulate to perform Monte Carlo simulations. We can also integrate our simulations into Stata's power commands so that we can easily create custom tables and graphs for a range of parameter values.

This FAQ is organized as follows:

  1. The general procedure
  2. How to use the program and simulate commands
  3. Verify your program using the null case
  4. How to integrate the simulation program into the power command

1. The general procedure

Statistical power is the probability of rejecting the null hypothesis when the null hypothesis is false. The general procedure for calculating power using Monte Carlo simulations is

  1. Use the underlying model to generate random data with (a) specified sample sizes, (b) parameter values based on the assumption that the alternative hypothesis is true (for example, mean = 75), and (c) nuisance parameters such as variances.
  2. Run the Stata estimation program (regress, glm, etc.) on these randomly generated data and obtain the p-value for the test we want.
  3. Save the p-value of the test and obtain the test result ("reject" or "fail to reject"). We reject the null hypothesis when the p-value is less than our significance level α.
  4. Repeat the steps 1–3 for N times (usually 1,000 or more).

The proportion of times that the null hypothesis is rejected (out of N) is our estimate of power.

Note: As a check, always run the simulation for the null case to make sure the estimate power is approximately α (to within the sampling error of the iteration number). More generally, in the null case, the distribution of the p-values should present a uniform distribution within the [0,1] interval.

2. How to use the program and simulate commands

Let's look at one example. Assume that we wish to test the effect of a drug dose on a Poisson-distributed response y in rats. For a sample of k rats, we have the Poisson regression model

log \( \mu_{i} = b_{1}x_{i} \)
\( y_{i} \tilde\ \mathrm{Poisson}(\mu_{i}) \)

for i=1, ..., k, where \( x_i \) is the dose in milligrams given to the \(i_{th} \) rat. Suppose we are trying to decide what sample size would give reasonable power for the test \( b_1=0 \) (the null hypothesis) if the true value of \( b_1 \) were 0.64 (the alternative hypothesis), under an experimental design with three levels of dosage (0.2, 0.5, and 1.0) repeated n times on \( k = 3n \) different rats.

Let's define a program to implement the steps 1–3 that we introduced above. For more instruction on how to use program to write a simple program in Stata, please read the documentation for [P] program.

capture program drop powersimu

program powersimu, rclass
version 17.0

   syntax, n(integer)	      /// sample size for each dosage group
	b1(real)             ///  b1 under the alternative
	[ alpha(real 0.05)   ///  alpha level
	d1(real 1)           ///  fixed dosage level 1
	d2(real 2)           ///  fixed dosage level 2
	d3(real 3)	     ///  fixed dosage level 3

   set obs 3
   generate x=`d1' in 1
   replace x=`d2' in 2
   replace x=`d3' in 3

   expand `n'

   generate mu=exp(`b1'*x)

// generate random numbers as dependent variable in the Poisson regression
   gen y = rpoisson(mu)        

// fit the Poisson regression 
   poisson y x, noconstant 

   mat a=r(table)
   local p1=el(a,rownumb(a,"pvalue"),colnumb(a,"x"))
   return scalar pvalue = `p1'
   return scalar reject = (`p1'<`alpha') 

You could save the above program in an ado-file called powersimu.ado. Now we can run the program powersimu with the input parameters, and it will give you the result of one iteration. Please note that the above program calls the random-number generator so that we use set seed to get reproducible results.

. set seed 1234
. powersimu, n(10) b1(0.64) d1(.2) d2(.5) d3(1.0) alpha(0.05)
. return list

The output is

             r(reject) =  1
             r(pvalue) =  .0000330350943392

In this iteration, the null hypothesis \( (b_1 = 0) \) was rejected because the p-value is smaller than α.

In order to estimate the power by repeating the above procedure, we can run this program with the simulate prefix.

. simulate reject=r(reject) pvalue=r(pvalue), reps(1000) seed(1234): 
> powersimu, n(10) b1(0.64) d1(.2) d2(.5) d3(1.0) alpha(0.05)

In the above line, we specify 1000 as the number of replications in the simulation. We also set the random-number seed for the simulation by specifying the seed() option in order to get reproducible results. After we run the simulation, the test results and p-values are saved in the new variables reject and pvalue, respectively. You can use summarize to calculate the mean of reject, which equals the proportion of times out of 1000 iterations that the null hypothesis was rejected. That proportion is your estimate of statistical power given the input parameters.

. summarize reject

Variable Obs Mean Std. dev. Min Max
reject 1,000 .788 .4089294 0 1

In this example, the proportion equals 0.788, which means that you can expect 78.8% power for this test when the coefficient of drug dose equals 0.64, given a sample size of 30. You could increase the number of replications in the simulation to get a more accurate power estimate.

3. Verify your program using the null case

Now, we would like to run the program with \( b1 = 0 \) with a relatively large value of N (replication number) to see if the proportion of rejections is close to α. If not, then either the normal approximation used by the command is not adequate for this case or we may have programmed the procedure incorrectly. Let's run the simulation with N = 1500.

. simulate reject=r(reject) pvalue=r(pvalue), reps(1500) seed(1234):
> powersimu, n(10) b1(0) d1(.2) d2(.5) d3(1.0) alpha(0.05)

Then, we can estimate the proportion of observations with reject = 1 using ci proportion. The point estimate of proportion is .0546667, and a 95% confidence interval is (0.0437109, 0.0674042), given α = 0.05.

. ci proportion reject

Binomial exact
Variable Obs Proportion Std. err. [95% conf. interval]
reject 1,500 .0546667 .0058696 .0437109 .0674042

Now, let's check the entire distribution of simulated p-values for a uniform distribution using histogram. We can see that it is approximately uniform distributed.

. histogram pvalue

These results suggest that the normal approximation z test is quite good for 30 observations (n = 10).

4. How to integrate the simulation program into the power command

With the power suite of commands, we can specify multiple values for each parameter such as the coefficient, the standard deviation, the means, or the α level. And power will create custom tables and graphs that present the results under multiple parameter values. In addition, with power we can add our own methods to the suite of commands. We will show you how to add the Poisson regression simulation program we just created to power. First, we need to define a program named power_cmd_mymethod. In this program, we need to save the power, sample size, and other parameters as the stored results. Because we already know how to obtain power using the program powersimu and simulate, the new program will simply call powersimu to obtain the power.

capture program drop power_cmd_powersimu
program power_cmd_powersimu, rclass
version 17.0         

   syntax, n(integer)	     /// sample size for each dosage group
	b1(real)             /// b1 under the alternative
	[ alpha(real 0.05)   /// Alpha level
	d1(real 1)           /// fixed dosage level 1
	d2(real 2)           /// fixed dosage level 2
	d3(real 3)	     /// fixed dosage level 3
	reps(integer 100)    /// Number of repetitions
// call the powersimu program
   quietly simulate reject=r(reject), reps(`reps') seed(12345):    ///
   powersimu, n(`n') b1(`b1') d1(`d1') d2(`d2') d3(`d3') alpha(`alpha')
   quietly summarize reject
   return scalar power = r(mean)
   return scalar N = 3*`n'
   return scalar alpha = `alpha'
   return scalar b1 = `b1'

You can save the above program in an ado-file called power_cmd_powersimu.ado.

We've almost finished the work now. But you will need to write one more small initializer program if you wish to enter a range of values for the parameters other than n() or alpha(). The program must be named power_cmd_mymethod_init, so we will name our program power_cmd_powersimu_init.

The initializer program is defined as a sclass program. In the program, the line sreturn local pss_colnames initializes columns in the output table for the parameters listed in double quotes. The line sreturn local pss_numopts allows you to specify numlists for the parameters placed in double quotes. For more detailed information regarding writing the initializer program for user-defined power commands, please read the entry to [PSS] power usermethod.

You can save the following program in an ado-file called power_cmd_powersimu_init.ado.

capture program drop power_cmd_powersimu_init
program power_cmd_powersimu_init, sclass
   sreturn clear
   sreturn local pss_numopts "b1"
   sreturn local pss_colnames "b1"

Now, we can use power powersimu to calculate power for a range of coefficient values assuming different alternative hypotheses and a range of sample sizes. Please note that power powersimu calls the program power_cmd_powersimu rather than powersimu. In the generated table, we can easily see the estimated power values for different input parameters. And we can plot the power analysis results by specifying the graph option. You can also change the x-dimension variable in your plot by specifying xdimension in the graph option. Please note that the sample size equals 3*n (repeat time in each dosage level) in this power analysis because we have three levels of drug dosage.

. program drop _all
. power powersimu, reps(1000) n(10(10)30) b1(0.5(0.05)0.85) d1(.2) d2(.5) d3(1.0)
> table graph(xdimension(b1)) alpha(0.05)
Estimated power
Two-sided test

alpha power N b1
.05 .58 30 .5
.05 .673 30 .55
.05 .745 30 .6
.05 .813 30 .65
.05 .849 30 .7
.05 .912 30 .75
.05 .945 30 .8
.05 .974 30 .85
.05 .842 60 .5
.05 .896 60 .55
.05 .951 60 .6
.05 .977 60 .65
.05 .985 60 .7
.05 .994 60 .75
.05 1 60 .8
.05 1 60 .85
.05 .946 90 .5
.05 .979 90 .55
.05 .992 90 .6
.05 .997 90 .65
.05 1 90 .7
.05 1 90 .75
.05 1 90 .8
.05 1 90 .85

As we can see in the graph, a larger true value of \( b_1 \) will lead to a larger power because there is a bigger chance to reject the null hypothesis \( b_1=0 \) when the true value of \( b_1 \) is larger. The graph also suggests that the increase in power is substantial between N = 30 and N = 60 and becomes less beneficial when N is beyond 60.

On the other hand, if we do not specify the xdimension of graph, the sample size N will be the xdimension by default. If we run the following power command with the notable option, the table will be suppressed and we will see only the graph.

. power powersimu, reps(1000) n(10(5)30) b1(0.5(0.1)0.7) d1(.2) d2(.5) d3(1.0) 
> table graph alpha(0.05) notable 

As expected, as sample size increases, power increases. As power gets closer to 1, the effect of increasing the sample size becomes smaller.





The Stata Blog: Not Elsewhere Classified Find us on Facebook Follow us on Twitter LinkedIn YouTube Instagram
© Copyright 1996–2023 StataCorp LLC   •   Terms of use   •   Privacy   •   Contact us