## 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