1
$\begingroup$

Suppose I generate two independent binomial samples, each of size $N$, with success probabilities $p_A$ and $p_B$. I then use R's t.test to test whether the difference in means is statistically significantly different from zero. What is the correct way (in R) to do a power calculation for such a test? I thought I should use power.t.test, but my simulation suggests that I am misunderstanding something:

library(pwr)

conversion_rate_A <- 0.10  # Success probability in group A
conversion_rate_B <- 0.15  # Success probability in group B
var_A <- conversion_rate_A * (1 - conversion_rate_A)  # Bernoulli variance = p*(1-p)
var_B <- conversion_rate_B * (1 - conversion_rate_B)
effect_size <- abs(conversion_rate_A - conversion_rate_B) / sqrt(var_A + var_B)

power <- 0.60  # Probability of rejecting null of no difference between A and B

pwr.t.test(d=effect_size, power=power, type="two.sample")
# Returns n = 853 for 0.05 significance
power.t.test(delta=conversion_rate_A - conversion_rate_B,
             sd=sqrt(var_A + var_B),
             sig.level=0.05,
             strict=TRUE,
             power=power)
# Returns n=853, agrees with pwr.t.test

## Sanity check: calculate power from simulation
test_is_significant <- function(conversion_rate_A, conversion_rate_B, n_obs=853) {
  outcomes_A <- rbinom(n=n_obs, size=1, prob=conversion_rate_A)
  outcomes_B <- rbinom(n=n_obs, size=1, prob=conversion_rate_B)
  test <- t.test(x=outcomes_A, y=outcomes_B,
                 alternative="two.sided",
                 paired=FALSE, var.equal=FALSE)
  return(test$p.value <= 0.05)  # Significance indicator (reject null of no difference)
}

test_is_significant(conversion_rate_A, conversion_rate_B)  # Single simulation

n_simulations <- 10^4
mean(replicate(n=n_simulations, 
               test_is_significant(conversion_rate_A, conversion_rate_B, n_obs=853)))
# Power estimate from simulation

My simulation is essentially a Monte Carlo power calculation. It should return something close to $0.60$, but I am getting numbers around $0.89$. Where is my mistake?


I have partly answered my own question: I should have used power.prop.test(p1=conversion_rate_A, p2=conversion_rate_B, power=power) # 428. If I set n_obs=428 in my replicate call, the simulation agrees with the power calculation.

Updated question: Why are these power calculations (power.prop.test versus power.t.test) so different? Does power.t.test assume a continuous distribution?

$\endgroup$
  • $\begingroup$ Related to stats.stackexchange.com/questions/83700/… ? $\endgroup$ – Adrian Nov 3 '17 at 0:08
  • 1
    $\begingroup$ The t-test does assume a continuous distribution (indeed a normal distribution -- which will have some impact in small samples), but more importantly the default two sample t-test in R is a Welch-Satterthwaite test (see the help, which makes that clear) while power.t.test assumes the equal-variance t-test -- and the impact of that difference is not necessarily small. First be clear about which test you want to do and then some advice may be possible. $\endgroup$ – Glen_b Nov 3 '17 at 0:41
  • $\begingroup$ @Glen_b I think what you're hinting at is that the correct test and power calculations, given the rbinom data, are prop.test and power.prop.test? Is using t.test(x, y, paired=FALSE, var.equal=FALSE) a mistake in this setting? My rough sense is that it would be fine given that the sample mean of binomial data (or the difference in two sample means, from two independent binomials) will be "close" to normally distributed, but I'm not being very rigorous. $\endgroup$ – Adrian Nov 4 '17 at 18:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.