Compromised power analysis estimates

In compromised power analyses, \(\alpha\) and \(1-\beta\) are computed as functions of the effect size, sample size, model, and most importantly the error probability ratio \(q=\frac{\beta}{\alpha}\) to indicate a type of trade-off between Type I and II errors. Setting \(\alpha=\beta\) results in a \(q=1\), indicating an equal trade-off between Type I and Type II errors, while a \(q=4\) indicates that Type II errors are four-times as costly to make than a Type I error (hence, rejecting the null when it is in fact true is worse than retaining the null hypothesis when it is in fact false). This is particularly useful when the sample size required from prior power analysis suggest a much larger size than a researcher can afford.

For instance, suppose that only \(N=100\) per group are possible for an independent samples \(t\)-test analysis with a Welch correction, and the true effect size is thought to be \(d=.3\). In such a simulation, and assuming for the moment that \(\alpha=.05\) (will be changed momentarily thanks to the store_results = TRUE flag), then the following indicates the \(q=\frac{\beta}{\alpha}\) error ratio estimate for this fixed \(\alpha\) level.

library(SimDesign)

Design <- createDesign(N = 100,
                       d = .3, 
                       alpha = .05)
Design    

# A tibble: 1 × 3
      N     d alpha
  <dbl> <dbl> <dbl>
1   100   0.3  0.05

# (optional) print RStudio special disabling flags and add to script
##   (fixes missing variables when Attach used)
Attach(Design, RStudio_flags = TRUE)

# !diagnostics suppress=N,d,alpha

# !diagnostics suppress=N,d,alpha

#~~~~~~~~~~~~~~~~~~~~~~~~
#### Step 2 --- Define generate, analyse, and summarise functions

Generate <- function(condition, fixed_objects = NULL) {
    Attach(condition)
    group1 <- rnorm(N)
    group2 <- rnorm(N, mean=d)
    dat <- data.frame(DV = c(group1, group2),
                      group = rep(c('G1', 'G2'), each=N))
    dat
}

Analyse <- function(condition, dat, fixed_objects = NULL) {
    p <- t.test(DV ~ group, data=dat)$p.value
    p
}

Summarise <- function(condition, results, fixed_objects = NULL) {
    rate <- EDR(results, alpha=condition$alpha, unname = TRUE)
    ret <- c(q = (1-rate) / condition$alpha)
    ret
}

#~~~~~~~~~~~~~~~~~~~~~~~~
#### Step 3 --- Compute q ratio

sim <- runSimulation(design=Design, replications=100000, 
                     generate=Generate, analyse=Analyse,
                     summarise=Summarise)

Replications: 100000;   RAM Used: 53.7 Mb;   
 Conditions: N=100, d=0.3, alpha=0.05

sim

# A tibble: 1 × 9
      N     d alpha     q REPLICATIONS SIM_TIME   RAM_USED       SEED COMPLETED 
  <dbl> <dbl> <dbl> <dbl>        <dbl> <chr>      <chr>         <int> <chr>     
1   100   0.3  0.05 8.789       100000 02m 10.13s 55.8 Mb  1930362837 Thu Apr 1…

However, because store_results=TRUE is used (the default) the results can be reSummarise()ed using a different \(\alpha\) cut-off, where it is instead possible to obtain some target \(q\) given the stored stimulation results. For example,

# compute beta/alpha ratio given different alpha
compromise <- function(alpha, sim, Design){
    Design$alpha <- alpha
    out <- reSummarise(Summarise, results=sim, Design=Design)    
    out$q
}

# more liberal Type I error (but lower Type II error)
compromise(.3, sim=sim, Design=Design)  

[1] 0.463

# more conservative Type I error (but higher Type II error)
compromise(.01, sim=sim, Design=Design) 

[1] 68.2

which indicates different \(q\) ratios. If a specific ratio is desired, then root-solving methods can be used to obtain the desired \(f(\alpha) = q\).

# define root function f(alpha) - q = 0
compromise_root <- function(alpha, target.q, ...)
    compromise(alpha, ...) - target.q

# solve alpha given equal beta/alpha trade-off  
root1 <- uniroot(compromise_root, c(.01, .3), target.q=1, 
        sim=sim, Design=Design)
root1

$root
[1] 0.2

$f.root
[1] -6.58e-05

$iter
[1] 8

$init.it
[1] NA

$estim.prec
[1] 6.1e-05

# solve alpha beta/alpha trade-off 4 times worse (beta = 4 * alpha)
root4 <- uniroot(compromise_root, c(.01, .3), target.q=4, 
        sim=sim, Design=Design)
root4

$root
[1] 0.0863

$f.root
[1] -0.00091

$iter
[1] 7

$init.it
[1] NA

$estim.prec
[1] 6.1e-05

Hence, when equal \(\beta\) and \(\alpha\) errors are desirable then the \(\alpha\) to utilize is 0.2, while if the Type I errors are 4 times more costly than the Type II errors then \(\alpha\) should be selected to be approximately 0.086.

Compromise analysis with empirical \(\alpha\) estimate

In situations where the Type I error rate controlled by a select \(\alpha\), but the true/empirical Type I error rate associated with this \(\alpha\) is in fact sub-optimal (e.g., small samples that utilize maximum-likelihood estimators), then it is possible to define the compromise ratio \(q = \frac{\beta}{\alpha}\) in terms of the empirical Type I error rate rather than the assumed nominal \(\alpha\). To do so requires more computation as the null model must also be generated and analysed, however the resulting compromise ratio and root-solved cut-offs should perform more honestly in practice. Below is an example which utilizes the empirical Type I error estimate rather than the assume Type I error = \(\alpha\).

#####################
# Same as above, however if Type I error not nominal then may wish to use 
# empirical Type I error estimate instead
library(SimDesign)

Design <- createDesign(N = 100,
                       d = .3, 
                       alpha = .05)
Design    

# A tibble: 1 × 3
      N     d alpha
  <dbl> <dbl> <dbl>
1   100   0.3  0.05

#~~~~~~~~~~~~~~~~~~~~~~~~
#### Step 2 --- Define generate, analyse, and summarise functions

Generate <- function(condition, fixed_objects = NULL) {
    Attach(condition)
    group1 <- rnorm(N)
    group2 <- rnorm(N, mean=d)
    group3 <- rnorm(N)   # For H0 tests
    dat <- data.frame(DV = c(group1, group2),
                      DV.null = c(group1, group3),
                      group = rep(c('G1', 'G2'), each=N))
    dat
}

Analyse <- function(condition, dat, fixed_objects = NULL) {
    p.power <- t.test(DV ~ group, data=dat)$p.value       # Power: (1-beta)
    p.null <- t.test(DV.null ~ group, data=dat)$p.value   # Type I error: alpha
    nc(p.null, p.power)
}

Summarise <- function(condition, results, fixed_objects = NULL) {
    rate <- EDR(results, alpha=condition$alpha)
    ret <- c(q = unname((1-rate["p.power"]) / rate["p.null"]))
    ret
}

#~~~~~~~~~~~~~~~~~~~~~~~~
#### Step 3 

sim <- runSimulation(design=Design, replications=100000, 
                     generate=Generate, analyse=Analyse,
                     summarise=Summarise, store_results=TRUE)

Replications: 100000;   RAM Used: 59.6 Mb;   
 Conditions: N=100, d=0.3, alpha=0.05

sim

# A tibble: 1 × 9
      N     d alpha      q REPLICATIONS SIM_TIME   RAM_USED      SEED COMPLETED 
  <dbl> <dbl> <dbl>  <dbl>        <dbl> <chr>      <chr>        <int> <chr>     
1   100   0.3  0.05 8.9691       100000 03m 39.79s 61.1 Mb  447974663 Thu Apr 1…

# compute beta/alpha ratio given different alpha
compromise <- function(alpha, sim, Design){
    Design$alpha <- alpha
    out <- reSummarise(Summarise, results=sim, Design=Design)    
    out$q
}

compromise(.3, sim=sim, Design=Design)

[1] 0.46

compromise(.01, sim=sim, Design=Design)

[1] 71.1

compromise_root <- function(alpha, target.q, ...)
    compromise(alpha, ...) - target.q


# equal beta/alpha trade-off  
root1 <- uniroot(compromise_root, c(.01, .3), target.q=1, 
        sim=sim, Design=Design)
root1

$root
[1] 0.2

$f.root
[1] -0.00015

$iter
[1] 8

$init.it
[1] NA

$estim.prec
[1] 6.1e-05

# beta/alpha trade-off 4 times worse for alpha (beta = 4 * alpha)
root4 <- uniroot(compromise_root, c(.01, .3), target.q=4, 
        sim=sim, Design=Design)
root4

$root
[1] 0.0868

$f.root
[1] -0.00139

$iter
[1] 7

$init.it
[1] NA

$estim.prec
[1] 6.1e-05

Based on the empirical \(\hat{\alpha}\) and \(\hat{\beta}\) estimates, when equal \(\beta\) and \(\alpha\) errors are desirable then the \(\alpha\) to utilize is 0.2, while if the Type I errors are 4 times more costly than the Type II errors then \(\alpha\) should be selected to be approximately 0.087. These agree with the results above for this particular analysis because the empirical \(\alpha\) closely matches the theoretical \(\alpha\), however in situations where these quantities do not match the solutions can and will differ.