Power calculation with a binary outcome (using a R package & manual simulation)

Question

I am trying to compare three different methods for power calculation under the same setting. Although I understand that the power estimates from each method cannot be exactly identical due to randomness, I expected them to be very similar. However, I observed somewhat different power estimates.

The basic setting is as follows:

1. Outcome: Binary (success 1; failure 0)
2. Two groups to compare: G1 vs G2
3. G1: # of samples - 150; a rate of success 0.2
4. G2: # of samples - 30; a rate of success 0.4
5. Significance level of 0.2 (not the usual 0.05)

That is,

p1 <- 0.2
p2 <- 0.4

n1 <- 150
n2 <- 30

The three methods that I used were:

1. Using a function of pwr.2p2n.test in the pwr package.
2. Simulations using a function of prop.test.
3. Simulations using a function of fisher.test

### -------------------------------- ###
### --- Version 1: pwr.2p2n.test --- ###
### -------------------------------- ###
library(pwr)
pwr.2p2n.test(h = ES.h(p1 = p1, p2 = p2), 
              n1 = n1, n2 = n2,
              sig.level = 0.20,
              alternative = "less")

# power = 0.9145152

### ---------------------------- ###
### --- Version 2: Prop test --- ###
### ---------------------------- ###
nreps <- 10000
y1 <- rbinom(n = nreps, size = n1, p = p1)
y2 <- rbinom(n = nreps, size = n2, p = p2)

pval <- rep(NA, nreps)
for(i in 1:nreps) {
  pval[i] <- prop.test(c(y1[i], y2[i]), 
                       n= c(n1, n2), 
                       alternative = "less",
                       p = NULL, correct = TRUE)$p.value
}

power <- sum(pval < 0.20) / nreps 
power # [1] 0.8756

### -------------------------------------- ###
### --- Version 3: Fisher's exact test --- ###
### -------------------------------------- ###
pval.fin <- c()
for(i in 1:nreps){
  y1 <- rbinom(n1, size = 1, p = p1)
  y2 <- rbinom(n2, size = 1, p = p2)
  dat.comb <- rbind(data.frame(res = y1, group = "G1"),
                    data.frame(res = y2, group = "G2"))
  tab <- table(dat.comb$group, dat.comb$res)
  test.res <- fisher.test(tab, alternative = 'less')
  pval.fin[i] <- test.res$p.value
}
mean(pval.fin<0.2) # [1] 8e-04

Here are two things that I noticed:

1. The third approach that uses Fisher's exact test produces a power of 8e-0.4. However, if I change the 'alternative' option from 'less' to 'greater', the power becomes 0.8811, which is similar to what we obtain from the first and second approaches. The first and second approaches used the 'alternative' option of 'less' though... I do not understand what causes this huge difference.

2. The power estimates from Version 1 and Version 2 are also a bit off. I am not quite sure what causes this difference.

Can anyone help me understand why I am not getting similar power even with the identical setting?

Answer 1

As was raised in the comments:

pwr::pwr.2p2n.test (and its balanced sibling) does not apply a Yates correction, so you can get prop.test to match if you include correct = FALSE there.

Fisher's exact test is a different construction entirely, so it's unlikely that the power will ever match exactly. However, as stated in its documentation:

The alternative for a one-sided test is based on the odds ratio, so alternative = "greater" is a test of the odds ratio being bigger than or.

Since you are providing your data as non-response and G1 first, the calculated odds ratio is equal to $\frac{(1-p_1)/p_1}{(1-p_2)/p2}=\frac{8}{3}$ which you're currently testing against being less than the null value of $1$, and it quite clearly isn't. To get the intended test of $p_1<p_2$ you either have to swap rows, swap columns, or use alternative = "greater" in your input. Any one of these will effectively swap the numerator & denominator in your odds (ratio) calculation.