Introduction

We have recently, August 2020, added 5 new functions to Superpower to expand the package’s capabilities. In this vignette, we will introduce these new functions and how they can help with study planning.

These new functions include: power.ftest, optimal_alpha, ANOVA_compromise, alpha_standardized, and p_standardized

The overall goal is to provide ways for researchers to justify their alpha level when designing studies, and provide greater flexibility.

How are error rates minimized?

Within the new optimal_alpha function, there are two options within the error argument: minimal or balanced.

First, we have an equation for how a weighted combined error rate is calculated when error = minimal. (costT1T2 * x + priorH1H0 * y)/(priorH1H0 + costT1T2)

\(WCER = \frac{\frac{cost_\alpha}{cost_\beta} \cdot \alpha + \frac{Pr(H_1)}{Pr(H_0)}\cdot\beta}{\frac{Pr(H_1)}{Pr(H_0)}+\frac{cost_\alpha}{cost_\beta}}\)

Second, we have an equation for how a weighted combined error rate is calculated when error = balanced. max((costT1T2 * x - priorH1H0 * y)/(priorH1H0 + costT1T2), (priorH1H0 * y - costT1T2 * x)/(priorH1H0 + costT1T2))

\(WCER = max(\frac{\frac{cost_\alpha}{cost_\beta} \cdot \alpha - \frac{Pr(H_1)}{Pr(H_0)}\cdot\beta}{\frac{Pr(H_1)}{Pr(H_0)}+\frac{cost_\alpha}{cost_\beta}} | \frac{ \frac{Pr(H_1)}{Pr(H_0)}\cdot\beta - \frac{cost_\alpha}{cost_\beta} \cdot \alpha }{\frac{Pr(H_1)}{Pr(H_0)}+\frac{cost_\alpha}{cost_\beta}})\)

Minimizing Error Rates

Assume we plan to perform an independent samples t-test, where our smallest effect size of interest is d = 0.5, and we are planning to collect 64 participants in each condition. We would normally calculate power as follows:

power.t.test(delta = .5, sd = 1, n = 64, sig.level = 0.05, type = 'two.sample', alternative = 'two.sided')$power

This analysis tells us that we have 80% power with a 5% alpha level for our smallest effect size of interest, d = 0.5, when we collect 64 participants in each condition.

If we design 2000 studies like this, the number of Type 1 and Type 2 errors we make depend on how often the null hypothesis is true, and how often the alternative hypothesis is true. Let’s assume both are equally likely for now. This means that in 1000 studies the null hypothesis is true, and we will make 50 Type 1 errors. In 1000 studies the alternative hypothesis is true, and we will make 100-80 = 20% Type 2 errors, so in 200 studies we will not observe a significant result even if there is a true effect. Combining Type 1 and Type 2 errors, in the long run, we should expect 250 of our 2000 studies to yield an error.

The goal in Neyman-Pearson hypothesis testing is to control the number of errors we make, as we perform hypothesis tests. Researchers often rely on convention when setting error rates, and there is no special reason to set the Type 1 error rate at 5% and the Type 2 error rate at 20%, and there might be better choices when designing studies. For example, when collecting 64 participants per condition, it is also not optimally efficient to use these norms.

res <- optimal_alpha(power_function = "power.t.test(delta = .5, sd = 1, n=64, 
                     sig.level = x, type='two.sample', 
                     alternative='two.sided')$power")

print(res)

## 
## Optimal Alpha
##        Alpha      Beta Weighted.Combined.Error.Rate
## 1 0.09967232 0.1216627                    0.1106675
## 
## From power function:
## [1] "power.t.test(delta = .5, sd = 1, n=64, \n                     sig.level = x, type='two.sample', \n                     alternative='two.sided')$power"

If a researcher is interested in effects of d = 0.5, and plans to collect 64 participants in each condition, setting the Type 1 error rate to 10% will increase the power to 88%. If we would perform 2000 studies designed with these error rates, we would observe 100 Type 1 errors and 120 Type 2 errors. The combined error rate across 2000 studies is 220 instead of 250. In other words, by choosing a more optimal alpha level, we can design lines of research more efficiently, because we are less likely to make errors in our statistical inferences.

res2 <- optimal_alpha(power_function = "power.t.test(delta = .5, sd = 1, n=64, 
                      sig.level = x, type='two.sample', 
                      alternative='two.sided')$power", error = "balance")

print(res2)

## 
## Optimal Alpha
##       Alpha     Beta Weighted.Combined.Error.Rate
## 1 0.1110855 0.111082                 1.763508e-06
## 
## From power function:
## [1] "power.t.test(delta = .5, sd = 1, n=64, \n                      sig.level = x, type='two.sample', \n                      alternative='two.sided')$power"

Relative costs and prior probabilities

So far we have assumed a Type 1 error and Type 2 error are equally problematic. But you might believe Cohen (1988) was right, and Type 1 errors are exactly 4 times as bad as Type 2 errors. Or you might think they are twice as problematic, or 10 times as problematic. However you weigh them, as explained by Mudge et al., 2012, and Ulrich & Miller, 2019, you should incorporate those weights into your decisions.

The function has another optional argument, costT1T2, that allows you to specify the relative cost of Type1:Type2 errors. By default this is set to 1, but you can set it to 4 (or any other value) such that Type 1 errors are 4 times as costly as Type 2 errors. This will change the weight of Type 1 errors compared to Type 2 errors, and thus also the choice of the best alpha level.

res3 <- optimal_alpha(power_function = "power.t.test(delta = .5, sd = 1, n=100, 
                      sig.level = x, type='two.sample', 
                      alternative='two.sided')$power",
                      error = "minimal", costT1T2 = 4)

print(res3)

## 
## Optimal Alpha
##        Alpha      Beta Weighted.Combined.Error.Rate
## 1 0.01918734 0.1211773                   0.03958533
## 
## From power function:
## [1] "power.t.test(delta = .5, sd = 1, n=100, \n                      sig.level = x, type='two.sample', \n                      alternative='two.sided')$power"

Now, the alpha level that minimized the weighted Type 1 and Type 2 error rates is 2%.

Similarly, you can take into account prior probabilities that either the null is true (and you will observe a Type 1 error), or that the alternative hypothesis is true (and you will observe a Type 2 error). By incorporating these expectations, you can minimize or balance error rates in the long run (assuming your priors are correct). Priors can be specified using the priorH1H0 argument, which by default is 1 (H1 and H0 are equally likely). Setting it to 4 means you think the alternative hypothesis (and hence, Type 2 errors) are 4 times more likely than that the null hypothesis is true (and hence, Type 1 errors).

res4 <- optimal_alpha(power_function = "power.t.test(delta = .5, sd = 1, n=100, 
                      sig.level = x, type='two.sample', 
                      alternative='two.sided')$power", 
                      error = "minimal", priorH1H0 = 2)

print(res4)

## 
## Optimal Alpha
##        Alpha       Beta Weighted.Combined.Error.Rate
## 1 0.07901659 0.03875692                   0.05217681
## 
## From power function:
## [1] "power.t.test(delta = .5, sd = 1, n=100, \n                      sig.level = x, type='two.sample', \n                      alternative='two.sided')$power"

If you think H1 is four times more likely to be true than H0, you need to worry less about Type 1 errors, and now the alpha that minimizes the weighted error rates is 8%. It is always difficult to decide upon priors (unless you are Omniscient Jones) but even if you ignore them, you are making the decision that H1 and H0 are equally plausible.

Extending to the Analysis of Variance

Now that the basics of finding an optimal alpha are out of the way we can move on to how these methods can be used for an ANOVA, and F-tests in general.

This is accomplished by the replacing the power.t.test function with the power.ftest function. As you can see below it operates very similar to the power.t.test; the main differences are that power.ftest requires arguments for the degrees of freedom rather than sample size. Further, rather than allow for input of the power as an argument, power_ftest allows the user to input the beta_level.

power.ftest(
  num_df = NULL,
  den_df = NULL,
  cohen_f = NULL,
  alpha_level = Superpower_options("alpha_level"),
  beta_level = NULL,
  liberal_lambda = Superpower_options("liberal_lambda")
)

The use of the function on its own is simple, but not very intuitive for power analysis on its own.

For example, let us suppose we have known degrees of freedom of 1 and 128 respectively, and an effect size of interest of cohen_f = .22. Now we simple want to find what the power would be given an alpha_level of .07.

power.ftest(
  num_df = 1,
  den_df = 128,
  cohen_f = .22,
  alpha_level = .045,
  beta_level = NULL,
  liberal_lambda = FALSE
)

## 
##      Power Calculation for F-test 
## 
##          num_df = 1
##          den_df = 128
##         cohen_f = 0.22
##     alpha_level = 0.045
##      beta_level = 0.3210373
##           power = 67.89627

The process is simple for using power.ftest; now we can see where it is actually useful.

design_result <- ANOVA_design(design = "2b*2w",
                   n = 40, 
                   mu = c(1.03, 1.41, 0.98, 1.01), 
                   sd = 1.03, 
                   r = 0.8, 
                   labelnames = c("voice", "human", "robot", 
                                  "emotion", "cheerful", "sad"),
                   plot = TRUE)

exact2_res = ANOVA_exact2(design_result, verbose = FALSE)

knitr::kable(exact2_res$main_results)

power.ftest(num_df = exact2_res$anova_table$num_df[3],
            den_df = exact2_res$anova_table$den_df[3],
            cohen_f = exact2_res$main_results$cohen_f[3],
            alpha_level = NULL,
            beta_level = .2)

## 
##      Power Calculation for F-test 
## 
##          num_df = 1
##          den_df = 78
##         cohen_f = 0.2707638
##     alpha_level = 0.125566
##      beta_level = 0.2
##           power = 80

res5 = optimal_alpha("power.ftest(num_df = exact2_res$anova_table$num_df[3],
            den_df = exact2_res$anova_table$den_df[3],
            cohen_f = exact2_res$main_results$cohen_f[3],
            alpha_level = x)$power/100",
            error = "minimal")

print(res5)

## 
## Optimal Alpha
##       Alpha      Beta Weighted.Combined.Error.Rate
## 1 0.1417106 0.1826088                    0.1621597
## 
## From power function:
## [1] "power.ftest(num_df = exact2_res$anova_table$num_df[3],\n            den_df = exact2_res$anova_table$den_df[3],\n            cohen_f = exact2_res$main_results$cohen_f[3],\n            alpha_level = x)$power/100"

comp_res = ANOVA_compromise(design_result,
                            emm = TRUE)

#ANOVA results
knitr::kable(comp_res$aov_comp)

#MANOVA results
knitr::kable(comp_res$manova_comp)

#emmeans results
knitr::kable(comp_res$emmeans_comp)

comp_res$aov_plotlist

## $voice

## 
## $emotion

## 
## $`voice:emotion`