knitr::opts_chunk$set( # collapse = TRUE, comment = "#>" )

library("Superpower")

Recently, @shieh_ancova demonstrated that most software uses a slightly flawed approach to estimating power for ANCOVAs, and focused on one-way ANOVA designs to compare his method to the method first mentioned by @cohen_book.

As @shieh_ancova eloquently points out in their simulations, the method of @cohen_book overestimates power and the problem is exacerbated by having a high number of covariates with high proportion of the variance ($R^2$) explained by the covariates. The problem is worst (30% error between estimated and actual power) in the simulation when there are 10 covariates included in the model when explained variance is approximately 81% ($\rho = 0.9$) [@shieh_ancova see Table 3). While this may not seem like much of issue if you don't expect to encounter this scenario, it still demonstrates the @cohen_book method is inconsistent in producing appropriate estimates of power. I believe this is reason enough (when provided with a sustainable and implementable alternative) to abandon the old method of @cohen_book for a newer method that provides exact, rather than approximate, estimates of power for ANCOVA.

Thankfully, @shieh_ancova was diligent in his work and demonstrated, showing both the math and simulations, how a new exact method could be utilized. The direct method described in the paper is implemented in the `power_oneway_ancova`

function.

We can copy the example from @maxwell_delaney that Shieh also used. In this example there are 3 groups with means (`mu`

) of 400, 450, 500 respectively. The error variance is 10000 (`sd = 100`

). Rather than simulating dozens of examples, I will demonstrate one scenario below where there are 3 covariates, and the $R^2$ is equal to 0.5 (treatment effect excluded). This is demonstrated in @shieh_ancova, Table 2.

For `power_oneway_ancova`

we can demonstrate both the approximate and exact methods using the `type`

argument. We can leave the `n`

argument out in order to solve for the sample size required to reach 80% power. Please notice that `round_up`

is set to TRUE since we want have a whole number for sample sizes (rather than a fractional sample size).

power_oneway_ancova( mu = c(400,450,500), n_cov = 3, sd = 100, r2 = .25, alpha_level = .05, #n = c(17,17,17), beta_level = .2, round_up = TRUE, type = "approx" )

Notice that this method requires 3 more subjects in order to achieve a minimum of 80% power.

power_oneway_ancova( mu = c(400,450,500), n_cov = 3, sd = 100, r2 = .25, alpha_level = .05, #n = c(17,17,17), beta_level = .2, round_up = TRUE, type = "exact" )

Now, @shieh_ancova mentioned something very interesting at the end of section 3.

"Although the prescribed application of general linear hypothesis is discussed only from the perspective of a one-way ANCOVA design, the number of groups G may also represent the total number of combined factor levels of a multi-factor ANCOVA design. Hence, using a contrast matrix associated with a specific designated hypothesis, the same concept and process of assessing treatment effects can be readily extended to two-way and higher-order ANCOVA designs."

Therefore, all that is needed to extend the one-way ANOVA code provided by @shieh_ancova is to provide the appropriate contrast matrix for the main effect or interaction ANOVA-level effect that is desired. `Superpower`

accomplishes this with the `ANCOVA_analytic`

function which internally uses the `model.matrix`

function to form the appropriate contrast matrix.

This function operates similar to the `ANOVA_power`

and `ANOVA_exact`

functions. However, the `ANCOVA_analytic`

function doesn't require the use of `ANOVA_design`

first and relies upon the closed formulas from @shieh_ancova rather than a simulation to calculate statistical power. Please note that unlike the `power_oneway_ancova`

function there is no option to apply the approximation from @cohen_book for factorial designs.

We can extend the previous scenario with 3 groups to a factorial design with 2 groups across 3 conditions.

# Run function res1 = ANCOVA_analytic( design = "2b*3b", mu = c(400, 450, 500, 400, 500, 600), n_cov = 3, sd = 100, r2 = .25, alpha_level = .05, #n = 17, beta_level = .2, round_up = TRUE ) # Print main results res1

The results can also be printed as `power.htest`

objects by accessing the individual effects in the `res1`

object.

res1$aov_list$a res1$aov_list$b res1$aov_list$ab

We can also check the design by using the plot method.

```
plot(res1)
```

However, you may want to compare the power of ANOVA to an ANCOVA. In that case you can use the `ANOVA_design`

function and pass it onto the `ANCOVA_analytic`

function. But, this forces you to set the sample size.

des1 = ANOVA_design( design = "2b*3b", mu = c(400, 450, 500, 400, 500, 600), n = 17, sd = 100) res2 = ANCOVA_analytic( design_result = des1, n_cov = 3, r2 = .25, alpha_level = .05, round_up = TRUE ) res2

User specified contrasts can also be used for a power analysis. These can be provided in the `cmats`

argument of the `ANCOVA_analytic`

function or supplied directly to the `ANCOVA_contrast`

function as an independent test. The `ANCOVA_analytic`

function requires that the contrasts matrices are provided as matrices (i.e., `as.matrix`

). The contrasts can be accessed in the `con_list`

part of the results and can be named (in this case "test").

ANCOVA_analytic(design = "2b", mu = c(0,1), n = 15, cmats = list(test = matrix(c(-1,1), nrow = 1)), sd = 1, r2 = .2, n_cov = 1)$con_list$test # Same result ANCOVA_contrast(cmat = c(-1,1), n = 15, mu = c(0,1), sd = 1, r2 = .2, n_cov = 1)

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.