# How to calculate the power of a test that compares two proportions

I'm trying to compute the power of a proportion test by hand. The null and alternative hypotheses are below and I'm using a significance level, $\alpha$ = 0.05

$H_0: p_1 = p_2$

$H_a: p_1 \neq p_2$

Let's assume the true value of $p_1$ = 0.69 and that the true value of $p_2$ = 0.68. Based on this, we can calculate the estimator of $p_i$ as:

$\hat{P}_i$ = $X_i/n$ where $X_i \sim~ B(n, p_i)$. Here $B$ is the Binomial distribution parameterized by the sample size, $n$ and the probability, $p_i$

I went ahead and simulated both distributions as follows:

p_1 <- 0.69
p_2 <- 0.68
num_samples <- 10000
sample_size <- 10000
set.seed(1231421)

p_1_hat <- rbinom(num_samples, size = sample_size, prob = p_1)/sample_size
p_2_hat <- rbinom(num_samples, size = sample_size, prob = p_2)/sample_size


Now to calculate power, I first find out the 0.025 and 0.975 quantile of $\hat{P}_1$ as follows:

> quantile(p_1_hat, c(0.025, 0.975))
2.5%  97.5%
0.6809 0.6990


At this point, the power calculation should be straightforward:

> pnorm(quantile(p_1_hat, 0.025), p_2, sd(p_2_hat)) +
+ pnorm(quantile(p_1_hat, 0.975), p_2, sd(p_2_hat), lower.tail = FALSE)
0.576146


Unfortunately, this does not match with the output from the pwr library.

> library(pwr)
> pwr.2p.test(ES.h(p_1, p_2), n = sample_size)

Difference of proportion power calculation for binomial distribution (arcsine transformation)

h = 0.0215284
n = 10000
sig.level = 0.05
power = 0.3310592
alternative = two.sided

NOTE: same sample sizes


Where am I going wrong? A plot of the way I'm trying to calculate power is shown below. I'm trying to find the area to the left of the vertical line near 0.68 and to the right of the vertical line near 0.70 (practically zero). Isn't that what power is?

Edit: Calculations using the arcsine transform based on gammer's comments.

p_1_hat_asin <- asin(sqrt(p_1_hat))
p_2_hat_asin <- asin(sqrt(p_2_hat))
> pnorm(quantile(p_1_hat_asin, 0.025), mean(p_2_hat_asin), sd(p_2_hat_asin))
0.5752689

• Harvey Motulsky, you changed the meaning of the question when you altered the title. The question clearly describes comparing two proportions. It's not about comparing a sample proportion to a specific value. Did you even read the question, look at code to simulate this scenario, or notice the call to pwr.2p.test? Be more careful. Jan 29, 2017 at 1:59
• I don't understand the logic of your calculation. Under the null hypothesis (of equal proportions), the difference between the transformed proportions has mean zero and $\sigma=\sqrt{1/20000}$, in which case the critical value (for $\alpha=0.05$) is .0138; given the values you supplied, the power should be approximately p_1 = 0.69; p_2 = 0.68; greater = 1-pnorm(.0138, mean=asin(sqrt(p_1))-asin(sqrt(p_2)), sd=sqrt(1/20000)); less = pnorm(-.0138, mean=ap_1-ap_2, sd=sqrt(1/20000)); greater + less; Jan 30, 2017 at 2:35
• That's a bit off what the power calc gave you (but still very close), so maybe something slightly different is happening under the hood. I suggest doing some googling, or trying to get your hands on Cohen's textbook, where the details of the test may be more clearly described. My guess that it was just based on the difference of arcsine (square root) transformed proportions was based just on what i know about cohen's H. Good luck. Jan 30, 2017 at 2:38
• p.s. ap_1 and ap_2 are the same as asin(sqrt(p_1)) and asin(sqrt(p_2)). Not sure why I re-calculated them in the first call to pnorm. Oh well. Anyway, hope that helps. Jan 30, 2017 at 6:28

The power calculator you're using operates on the arcsin transformed scale (which is what the cohen's H is based on). So, you'd need to use the approximate normality of the arcsine transformed proportions and do the analytic calculation on that scale to recover the result of pwr.2p.test. I'll leave that analysis to you. It shouldn't be too complicated. Just use the fact that difference between two independent approximately normal variables is approximately normal and calculate areas under that distribution. Give it a shot and post the answer if you figure it out. Good luck.

As a sidebar, below is some code to estimate the power of the $\chi^2$ test of equality of proportions by simulation by calculating the proportion of times (out of a large number, 1000 in this case) the null hypothesis is rejected.

set.seed(1122)
p1 = 0.69
p2 = 0.68
n = 10000
n.reject = 0
n.rep = 1000
for(k in 1:n.rep)
{
x1 = sum(runif(n)<p1)
x2 = sum(runif(n)<p2)
ptest = prop.test( c(x1,x2), c(n,n), correct=FALSE)
if(  ptest$p.value < 0.05 ) n.reject = n.reject + 1 } n.reject/n.rep [1] 0.322  • Couple of follow up questions. a. Why is the transformation required? Is it because$\hat{P}$is not normal? When I do a QQ plot against the standard normal, it seems pretty close to normal. b. I'm still getting the same answer when I do the arcsin transform. Since the transform is monotonic, wouldn't that be expected? I've added the code that I used to calculate the power via the transform in the OP. Jan 30, 2017 at 1:39 • I think the transformation is used so that you can use a single effect size measure (cohen's H) to capture differences between proportions that are, in some sense, invariant to the location (similar to cohen's D), because the variance, on the transformed scale, is (approximately) the same regardless of the value for$p\$. I'm not sure why you're getting the wrong answer but make sure you didn't make any calculation errors. See the wikipedia page about the arcsin transform of proportions: en.wikipedia.org/wiki/… Jan 30, 2017 at 1:48

Addressing this part of the question first:

Where am I going wrong? A plot of the way I'm trying to calculate power is shown below. I'm trying to find the area to the left of the vertical line near 0.68 and to the right of the vertical line near 0.70 (practically zero). Isn't that what power is?

Power should be calculated relative to the alternate hypothesis, rather than the null hypothesis, as we're calculating the probability of correctly rejecting the null hypothesis assuming that the alternate hypothesis is true. So the plot should look more like this:

And we're looking to calculate the area shaded in Pink (to the left of the black dashed line that is $$-Z_{\alpha/2}$$, and to the right of the other black dashed line that is $$Z_{\alpha/2}$$, that as you've mentioned, is close to 0)

As for the calculations, there are 2 ways to calculate it from the graph:

1. We can find the equivalent Z values on the Null hypothesis curve: $$Z = \frac{p_2 - p_1}{\sqrt{p_{pooled} * (1 - p_{pooled}) * (\frac{1}{n} + \frac{1}{n})}}$$ = $$\frac{0.68 - 0.69}{\sqrt{0.685 * (1 - 0.685) * (\frac{1}{10000} + \frac{1}{10000})}}$$ = $$-1.52224$$

$$Z_{\beta L} = -Z_{\alpha/2} + Z = -1.95996 + 1.52224 = -0.43772$$

$$Z_{\beta R} = Z_{\alpha/2} + Z = 1.95996 + 1.52224 = 3.4822$$

Power = Pink Area to the left + Pink Area to the right = scipy.stats.norm(0, 1).cdf(-0.43772) + (1 - scipy.stats.norm(0, 1).cdf(3.4822)) = 0.33104329

2. Or directly find the area under the Alternate hypothesis curve: scipy.stats.norm(-1.52224, 1).cdf(-1.95996) + (1 - scipy.stats.norm(-1.52224, 1).cdf(1.95996)) = 0.33104329

This is really close to R's pwr.2p.test answer of 0.3310592, and the reason for the difference is because of the methods used to calculate Power. As mentioned by @gammer, pwr uses the arcsin transformation, and calculates $$Z_\beta$$ with the following formula:

$$Z_\beta = Z_\alpha - \frac{h\sqrt{n}}{\sigma}$$,
where $$h = 2arcsin(\sqrt{p_1} - 2arcsin(\sqrt{p_2})$$,
and $$n = n/2$$ for equal $$n$$'s,
and $$n = \frac{1}{\frac{1}{n_1}+\frac{1}{n_2}}$$ for unequal $$n$$'s

Using this method:

$$h = 0.0215$$

$$Z_{\beta L} = 1.96 - \frac{0.0215 * \sqrt{5000}}{1} = -0.437$$

$$Z_{\beta R} = -1.96 - \frac{0.0215 * \sqrt{5000}}{1} = 3.482$$

Power = 1 - scs.norm(0, 1).cdf(3.482) + (scs.norm(0, 1).cdf(-0.437)) = 0.331059

Note: numbers have been truncated for brevity. But if you use the full number string, you would get the same numbers as pwr:

• h = 0.02152840303762038
• Z_alpha = ± 1.9599639845400545
• $$Z_{\beta L}$$ = 3.4822519621418984
• $$Z_{\beta R}$$ = -0.43767600693820996

Resources:
1. Documentation on pwr using the arcsin transformation
2. This other SE question outlining the 2 methods of calculating power
3. pwr source code: getAnywhere(pwr.2p.test)

I think that the n in the pwr.2p.test refers to the size of one of the two samples instead of the total size of the two samples.

If you double the sample size (i.e if you write: pwr.2p.test(ES.h(p_1, p_2), n = 2*sample_size) then you get a power of 0.5764 which is reasonably close to your other estimate.