One way to assess the power of a factorial anova design is through the use of Monte-Carlo simulation. With this approach, one generate hundreds (or thousands) of randomly generated datasets. Analyze each one and retain the p-values from the analyses. Then it is a simple matter to count up the proportion of p-values that are less than or equal to your nominal alpha level.

The user written program **factorialsim** (**search factorialsim**) will perform Monte-Carlo
power analyses for two-way factorial anova designs. To use **factorialsim** you need to create
a data file named **simdat.dta**. This data file will contain one row for each cell in
your design. Here is an example:

+----------------------+ | a b n m s | |----------------------| 1. | 1 1 8 20 9 | 2. | 1 2 8 27 9 | 3. | 1 3 8 19 9 | 4. | 2 1 12 20 10 | 5. | 2 2 12 25 10 | 6. | 2 3 12 30 10 | +----------------------+`list, sep(0)`

Variable **a** is the level of the first factor variable; variable **b** is the level
of the second factor variable; **n** is the number of observations in the cell; **m**
is the cell mean; and **s** is the cell standard deviation. The variables in the data file
must have the exact same names as shown above.

For each replication, **factorialsim** expands the **simdat** dataset on variable **n**.
Then, the program generates a random response variable with mean equal **m** and standard deviation
equal **s**. The program makes use of Stata’s
**simulate** command to collect and retain the Monte-Carlo results before
displaying the observed proportion of each of the p-values.

The examples below show how to use **factorialsim** for power and robustness analyses.

## Example 1 – Prospective Power Analysis

Let’s say that a researcher has decided that a 2×3 factorial design meets the need of his research project. Through a literature review and a pilot-study he has, what he thinks, are reasonable estimates for cell means and standard deviations. Further, he figures that it will be fairly easy to obtain 60 subjects for his studies. Of course, he will need to find more if the power analysis reveals low power especially for the important interaction hypothesis.

The researcher has also decided to put more subjects into the experimental groups than into the
control groups. Below is the code the researcher used to create the **simdat** data file.

`clear input a b n m s`

1 1 8 20 9 1 2 8 27 9 1 3 8 19 9 2 1 12 20 10 2 2 12 25 10 2 3 12 30 10`end save simdat, replace`

Now that the **simdat** data file has been created we will run the most basic **factorialsim**
command.

command: factorial_sim, obs(0) pab: r(pab) pa: r(pa) pb: r(pb) Simulations (100) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .................................................. 50 .................................................. 100 alpha = .05 replications = 100 simulated power for a*b = .54 simulated power for a = .25 simulated power for b = .48`factorialsim`

We won’t discuss these results because too few replications (only 100) were used.
We will rerun the command setting the **reps** option to 1,000. We will also
use the **nodots** option to suppress the display of the dots for each simulated
dataset.

command: factorial_sim, obs(0) pab: r(pab) pa: r(pa) pb: r(pb) alpha = .05 replications = 1000 simulated power for a*b = .486 simulated power for a = .229 simulated power for b = .415`factorialsim, reps(1000) nodots`

This proposed design seems to be under powered. We would like to get the power up to at
least 0.8. We can see what would happen if we were to increase each of the cell sizes to
20 observations using the **obs** command.

command: factorial_sim, obs(20) pab: r(pab) pa: r(pa) pb: r(pb) alpha = .05 replications = 1000 simulated power for a*b = .834 simulated power for a = .41 simulated power for b = .75`factorialsim, reps(1000) obs(20) nodots`

Now the power looks good for the interaction and is much better for the **b** main effect.
However, the **a** main effect still lacks power. Since the interaction has good power the
researcher must decide whether to stick with this sample size or increase it until the
main effect for **a** has a simulated power of 0.8.

It is also possible to run the analysis with the **robust** option. The **robust**
option causes **factorialsim** to run the model using **regress** with Huber-White
robust standard errors. Below we will rerun the analysis with **obs(20)** as in the
previous example.

command: factorial_simr, obs(20) pab: r(pab) pa: r(pa) pb: r(pb) alpha = .05 replications = 1000 simulated power for a*b = .838 simulated power for a = .056 simulated power for b = .748`factorialsim, reps(1000) nodots robust obs(20)`

The power for **a*b** and **b** held up very well. However, the power for **a**
has dropped dramatically due to heteroscedasticity between the two levels of **a**.

## Example 2 – Observed Power Analysis

Observed (post-hoc) power analysis is one that is done after the data have been collected and analyzed.

We begin by loading a dataset (**hsbdemo**) and running an anova.

Number of obs = 200 R-squared = 0.2590 Root MSE = 8.26386 Adj R-squared = 0.2399 Source | Partial SS df MS F Prob > F ------------+---------------------------------------------------- Model | 4630.36091 5 926.072182 13.56 0.0000 | female | 1261.85329 1 1261.85329 18.48 0.0000 prog | 3274.35082 2 1637.17541 23.97 0.0000 female#prog | 325.958189 2 162.979094 2.39 0.0946 | Residual | 13248.5141 194 68.2913097 ------------+---------------------------------------------------- Total | 17878.875 199 89.843593`use https://stats.idre.ucla.edu/stat/data/hsbdemo, clear`

`anova write female##prog`

Next, we create a dataset that contains the frequency, mean and standard deviation for each
cell using the **collapse** command.

a b m s n 0 1 49.14286 10.36478 21 0 2 54.61702 8.656622 47 0 3 41.82609 8.003705 23 1 1 53.25 8.205248 24 1 2 57.58621 7.115672 58 1 3 50.96296 8.341193 27`collapse collapse (mean) m=write (sd) s=write (count) n=write, by(female prog)`

`rename female a rename prog b`

list, clean nolabel noobs`save simdat, replace`

Now we can run **factorialsim** to get the Monte-Carlo post-hoc power estimates.

command: factorial_sim, obs(0) pab: r(pab) pa: r(pa) pb: r(pb) alpha = .05 replications = 1000 simulated power for a*b = .49 simulated power for a = .98 simulated power for b = 1`factorialsim, reps(1000) nodots`

The power for the **a** and **b** main effects are more than adequate while the power for the
interaction is a bit on the low side. If we play around with the **obs** option
we will see that it requires around 65 observation per cell to obtain a power of
.8 for the interaction this model.

command: factorial_sim, obs(65) pab: r(pab) pa: r(pa) pb: r(pb) alpha = .05 replications = 1000 simulated power for a*b = .814 simulated power for a = 1 simulated power for b = 1`factorialsim, reps(1000) obs(65) nodots`

## Example 3 – Robustness Analysis

To assess robustness of a model, you just set all the mean values to be equal. **factorialsim**
accomplishes this by setting all of the means to zero. Next you run the Monte-Carlo simulation.
If the portion of p-values is close to the nominal alpha level then the model displays good
robustness. Depending on cell sizes and variances the observed p-values can be much smaller
then alpha or much larger.

Just to show that robustness works well under optimal conditions, we will begin with an example in which all the cells are the same size and have the same standard deviations.

command: factorial_sim, obs(0) zero pab: r(pab) pa: r(pa) pb: r(pb) alpha = .05 replications = 1000 simulated robustness for a*b = .047 simulated robustness for a = .05 simulated robustness for b = .056`clear input a b n m s`

1 1 20 20 10 1 2 20 27 10 1 3 20 19 10 2 1 20 20 10 2 2 20 25 10 2 3 20 30 10`end save simdat, replace`

`factorialsim, reps(1000) nodots zero`

Note that all of the Monte-Carlo simulated alpha values are very close to the nominal 0.05 alpha levels.

Next, we will run an example in which cell size and variability differ. In this case the smaller cells in the design have the larger variability.

command: factorial_sim, obs(0) zero pab: r(pab) pa: r(pa) pb: r(pb) alpha = .05 replications = 1000 simulated robustness for a*b = .214 simulated robustness for a = .148 simulated robustness for b = .207`clear input a b n m s`

1 1 10 20 20 1 2 10 27 20 1 3 10 19 20 2 1 30 20 10 2 2 30 25 10 2 3 30 30 10 end save simdat, replace`factorialsim, reps(1000) nodots zero`

The observed proportions across our 1,000 samples is much larger that the nominal 0.05 level that is assumed.

Here’s what this analysis with the **robust** option looks like.

command: factorial_simr, obs(0) zero pab: r(pab) pa: r(pa) pb: r(pb) alpha = .05 replications = 1000 simulated robustness for a*b = .088 simulated robustness for a = .099 simulated robustness for b = .087`factorialsim, reps(1000) nodots zero robust`

The values with **robust** option are much better.

Now let’s run it again, this time with the smaller cells in the design having the lower variability.

command: factorial_sim, obs(0) zero pab: r(pab) pa: r(pa) pb: r(pb) alpha = .05 replications = 1000 simulated robustness for a*b = .003 simulated robustness for a = .005 simulated robustness for b = .002`clear input a b n m s`

1 1 10 20 10 1 2 10 27 10 1 3 10 19 10 2 1 30 20 20 2 2 30 25 20 2 3 30 30 20`end save simdat, replace`

`factorialsim, reps(1000) nodots zero`

In this case the observed proportions are all much lower than the nominal alpha level of 0.05.

Once again, we will rerun the analysis with the **robust** option.

command: factorial_simr, obs(0) zero pab: r(pab) pa: r(pa) pb: r(pb) alpha = .05 replications = 1000 simulated robustness for a*b = .055 simulated robustness for a = .064 simulated robustness for b = .088`factorialsim, reps(1000) nodots zero robust`

The p-values, once again, are much closer to the nominal 0.05 levels than without the
**robust** option.

Thus, **factorialsim** can be a useful tool in exploring robustness in two-way factorial
designs.