## Power analysis with simulation
## sample size per group
n <- 7 ## sample size: 4, 5, 6, 7
## assumptions about the model under the alternative
means <- c(7, 13, 10, 10, 10)
sigma <- sqrt(7.5)
## create data-frame
group <- factor(rep(LETTERS[1:5], each = n))
group
## test run: single data-set ####
y <- rnorm(n * 5, mean = rep(means, each = n), sd = sigma)
data <- data.frame(y = y, group = group)
## visualization
stripchart(y ~ group, data = data, vertical = TRUE, pch = 1)
## fit one-way ANOVA model
fit <- aov(y ~ group, data = data)
summary(fit)
summary(fit)[[1]][1, "Pr(>F)"] ## extract p-value of global test
## this was *one* data-set
## now have a look at *many* of such data-sets
## we simulate 1000 data-sets
## for every data-set we check whether the global null hypothesis
## is being rejected or not
## store the test decision in a vector of length 1000
results <- numeric(1000)
for(i in 1:1000){
## simulate new response
y <- rnorm(n * 5, mean = rep(means, each = n), sd = sigma)
data <- data.frame(y = y, group = group)
fit <- aov(y ~ group, data = data)
results[i] <- summary(fit)[[1]][1, "Pr(>F)"] < 0.05
}
mean(results)
## = proportion of cases where we actually reject H_0
## Results
# n = 4: 0.54
# n = 5: 0.69
# n = 6: 0.80
# n = 7: 0.87