Introduction to Power Analysis
The post is about statistical power analysis in R. First, define the meaning of power in statistics. The power is the probability ($1-\beta$) of detecting an effect given that the effect is here. Power is the probability of correctly rejecting the null hypothesis when it is false.
Suppose, a simple study of a drug-A and a placebo. Let the drug be truly effective. The power is the probability of finding a difference between two groups (drug-A and placebo group). Imagine that a power of $1-\beta=0.8$ (having a power of 0.8 means that 80% of the time, there will be statistically significant differences between the drug-A and the placebo group, whereas there are 20% of the time, the statistically significant effect will not be obtained between two groups). Also, note that this study was conducted many times. Therefore, the probability of a Type-II error is $\beta=0.2$.
One-Sample Power
The following plot is for a one-sample one-tailed greater than t-test. In the graph below, let the null hypothesis $H_0:\mu = \mu_0$ be true, and the test statistic $t$ follows the null distribution indicated by the hashed area. Under the specific alternative hypothesis, $H_1:\mu = \mu_1$, the test statistic $t$ follows the distribution shown by solid area.
The $\alpha$ is the probability of making a type-I error (that is rejecting $H_0$ when it is true), and the “crit. Val” is the location of the $t_{crit}$ value associated with $H_0$ on the scale of the data. The rejection region is the area under $H_0$ at least as far as $crit. val.” is from $\mu_0$.
The test’s power ($1-\beta$) is the green area, the area under $H_1$ in the rejection region. A type-II error is made when $H_1$ is true, but we fail to reject $H_0$ in the red region.
Type-II Error and Power Analysis in R
#One Sample Power
x <- seq(-4, 4, length = 1000)
hx <- dnorm(x, mean = 0, sd = 1)
plot(x, hx, type = "n", xlim = c(-4, 8), ylim = c(0, 0.5),
main = expression (paste("Type-II Error (", beta, ") and Power (", 1 - beta, ")")),
axes = FALSE)
# one-tailed shift
shift = qnorm (1 - 0.05, mean=0, sd = 1 )*1.7
xfit2 = x + shift
yfit2 = dnorm(xfit2, mean=shift, sd = 1 )
axis (1, at = c(-qnorm(0.05), 0, shift), labels = expression("crit. val.", mu[0], mu[1]))
axis(1, at = c(-4, 4 + shift), labels = expression(-infinity, infinity),
lwd = 1, lwd.tick = FALSE)
# The alternative hypothesis area
# the red - underpowered area
lb <- min(xfit2) # lower bound
ub <- round(qnorm(0.95), 2) # upper bound
col1 = "#CC2222"
i <- xfit2 >= lb & xfit2 <= ub
polygon(c(lb, xfit2[i], ub), c(0, yfit2[i],0), col = col1)
# The green area where the power is
col2 = "#22CC22"
i <- xfit2 >= ub
polygon(c(ub, xfit2[i], max(xfit2)), c(0, yfit2[i], 0), col = col2)
# Outline the alternative hypothesis
lines(xfit2, yfit2, lwd = 2)
# Print null hypothesis area
col_null = "#AAAAAA"
polygon (c(min(x), x, max(x)), c(0, hx, 0), col = col_null,
lwd = 2, density = c(10, 40), angle = -45, border = 0)
lines(x, hx, lwd = 2, lty = "dashed", col=col_null)
axis(1, at = (c(ub, max(xfit2))), labels = c("", expression(infinity)), col = col2,
lwd = 1, lwd.tick = FALSE)
#Legend
legend("topright", inset = 0.015, title = "Color",
c("Null Hypothesis", "Type-II error", "Power"), fill = c(col_null, col1, col2),
angle = -45, density = c(20, 1000, 1000), horiz = FALSE)
abline(v=ub, lwd=2, col="#000088", lty = "dashed")
arrows(ub, 0.45, ub+1, 0.45, lwd=3, col="#008800")
arrows(ub, 0.45, ub-1, 0.45, lwd=3, col="#880000")
