5
$\begingroup$

I am trying to test whether the sex ratio of some sampled individuals significantly differs from the expected sex ratio of 1. I have n= 64, of which female=34 and male=30.

I ran a binomial test:

succ <- c(34,30) 
binom.test(succ,0.5)

data:  succ
number of successes = 34, number of trials = 64, p-value = 0.708
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.4023098 0.6572035
sample estimates:
probability of success 
               0.53125 

I would like to calculate the statistical power of this test, and I know that power = 1-β, where β is the type II error.

I am getting confused when reading this explanation. I don't understand how to adapt this formula (for different choices of n) to my case:

enn = 1:2000
critical = qbinom(.025, enn, .5)
beta = pbinom(enn-critical,enn,.55) - pbinom(critical-1,enn,.55)

What I did was

1-(pbinom(34,64,0.5)- pbinom(30, 64, .5))
[1] 0.7410237

but I am not sure if it is correct to use 0.5 as probability. Moreover, I tried a different method, and I get a completely different result

pwr.p.test(ES.h(.53125,.5),n=64, power=NULL, alternative = "two.sided")

 proportion power calculation for binomial distribution (arcsine transformation) 

              h = 0.06254076
              n = 64
      sig.level = 0.05
          power = 0.07913605
    alternative = two.sided

Is one of these two tests correct and why?

Thanks for your help!

$\endgroup$
7
$\begingroup$

In order to find 'power', you need to have a specific alternative in mind. Suppose your null hypothesis is $H_0: p = 0.5$ vs. $H_a: p > 0.5,$ where $p = P(\mathrm{Female}).$ Also suppose you have $n = 64$ and you want the power of a test at level $\alpha = 0.05$ against the specific alternative $p = 0.6.$

For an exact binomial test, you need to find the critical value $c$ such that $P(X \ge c\,|\,n=64, p=.5)$ is maximized, but still below $0.05.$ In R, where dbinom, pbinom, and qbinom denote binomial PDF, CDF, and quantile function (inverse CDF), respectively, we see that the critical value is $c = 40.$ Notice that, because of the discreteness of binomial distributions, the so-called `5% level' actually rejects with probability $P(\mathrm{Rej}\, H_0 | H_0\, \mathrm{True}) \approx 3\%.$

qbinom(.95, 64, .5)
[1] 39
sum(dbinom(39:64, 64, .5))
[1] 0.05171094
sum(dbinom(40:64, 64, .5))
[1] 0.02997059
1 - pbinom(39, 64, .5)
[1] 0.02997059

Then the power of this test against alternative value $p = 0.6$ is given by $P(X \ge 40\,|\,n=64, p=0.6) = 0.3927.$

1 - pbinom(39, 64, .6)
[1] 0.392654

We can make a 'power curve' for this test by looking at a sequence of alternative values p.a between $0.5$ and $.75.$ The first block of R code below makes the solid black line in the plot below.

p.a = seq(.50, .75, by=.01)
p.rej = 1 - pbinom(39, 64, p.a)
plot(p.a, p.rej, type="l", main="Power Curve")
 abline(h=c(.03,1), col="green2")

enter image description here

If we look at a level $\alpha = 0.05$ test of $H_0: p = 0.5$ vs $H_a: p > 0.5$ with $n = 256$ subjects, then the critical value is $c = 141,$ the rejection probability when $H_0$ is true is $0.046,$ and the power against various alternative values of $p$ is greater, as shown by the dotted blue line in the figure.

c.256 = qbinom(.95, 256, .5); c.256
[1] 141
1 - pbinom(c.256, 256, .5)
[1] 0.04565604
p.rej.256 = 1 - pbinom(c.256, 256, p.a)
lines(p.a, p.rej.256, col="blue", lty="dotted")

Notes: Because $n = 64$ is sufficiently large to use normal approximations, you might want to try using normal approximations. A disadvantage is that this ignores the issue of discreteness, so it may appear that your test rejects exactly 5% of the time when $H_0$ is true. Also, you'd need to use a continuity correction for best results.

One relevant computation for the significance level in R is:

1 - pnorm(39.5, 32, 4)
[1] 0.03039636

(Approximate) power is $0.3895:$

mu.a = 64*.6;  sg.a = sqrt(64*.6*.4)
mu.a; sg.a
[1] 38.4
[1] 3.919184

1 - pnorm(39.5, mu.a, sg.a)     # Using NORM(mu.a, sg.a)
[1] 0.3894815
1 - pnorm((39.5 - mu.a)/sg.a)   # Standardizing and using NORM(0,1).
[1] 0.3894815
$\endgroup$
1
$\begingroup$

It's also important to ask yourself why you are calculating the power. Since you already have the data, you are calculating a 'post-hoc' power statistic (as opposed to a pre-hoc power statistic). It's worth nothing that a number of authors have criticised the use of post-hoc power statistics from a Frequentist point of view - see, for example, https://gpsych.bmj.com/content/32/4/e100069

$\endgroup$

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.