This is also a link Suppose we toss a coin (Bernoulli Experiment) a fixed number \(n\) times, say 40. We wish to test: \[ H_0: p = .5 \] vs \[ H_1: p > .5 \]
Let \(X\) be the outcome of this binomial experiment. Under \(H_0\), \(X \sim Bin(40, .5)\)
Note that this section is entirely optional if we just want to use binom.test function in R and extract p-values.
We reject \(H_0\) when \(X >= x_{crit}\). What is \(x_{crit}\)? For some alpha, we want \(x_{crit}\) such that: \[ \alpha = \mathbb{P}(\text{Reject } H_0 | H_0 \text{ True}) \] \[ \alpha = \mathbb{P}(X \geq x_{crit} | p = .5) \] So, \[ \alpha = \sum_{k = x_{crit}}^{40} \binom{40}{k} .5^k (1-.5)^{40-k} \] We can generate a list of possible cutoffs for the rejection region and find the corresponding values of \(\alpha\):
n <- 40
xcrits <- 1:40
round(sapply(xcrits, function(x) sum(dbinom(x:40, n, .5))), 3)## [1] 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.999 0.997 0.992 0.981 0.960 0.923 0.866 0.785 0.682 0.563 0.437
## [22] 0.318 0.215 0.134 0.077 0.040 0.019 0.008 0.003 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000Although we are not able to achieve an alpha of exactly \(.05\), rejecting when \(x_{obs} \geq 26\) yields a significance level of \(.04\) Under this rejection rule, we can compute the power of the test analytically.
Suppose, in truth, that \(p = .7\). Then the power of the test is actually: \[ 1 - \beta = \mathbb{P}(\text{Reject }H_0 | H_1 \text{ True}) \] \[ = \mathbb{P}(X \geq 26 | p = .7) \] \[ = \sum_{k=26}^{40} \binom{40}{k} .7^k (1-.7)^{40-k} \] Which we find in R:
sum(dbinom(26:40, 40, .7))## [1] 0.8074482Now, we run simulations where the true p is indeed \(.7\). We reject using the above rule. We replicated the experiment of \(40\) coin tosses \(10000\) times. The empirical power is:
reps <- 10000
sims_H1 <- rbinom(reps, n, .7)
mean(sims_H1 >= 26)## [1] 0.8164This matches with our analytical result above.
Note that if we did not want to manually find the rejection region as we did above, we could have used the built in R function binom.test():
mean(sapply(sims_H1 , function(x) binom.test(x, 40, alternative = "greater")$p.value) <.05)## [1] 0.8164Similarly, if we were running a linear model, we could use a similar approach: simulate linear models under different data, and extract those p-values. The analytical work is not neccessary to simulate the power. In this example we calculated the power knowing that \(n = 40\) in advance. However, we could have worked in the opposite direction, assumed we needed a power of say, \(.9\), and find the required sample size.