Difference between ANOVA power simulation and power calculation

I've been running some power simulations for a one-way ANOVA in R, and my problem is that the results from the simulation doesn't match the result from g*power or pwr.anova.test from the "pwr"-package. As an example, let's compare simulated power to analytical power using these values:

group_size <- c(40,40,40)
means <- c(0.2,0,-0.2)
sds <- c(1,1,1)


Analytical power analysis, f = 0.1632993 is calculated from the means and standard deviations above.

size <- 10
plot_df <- data.frame()
power <- 0
while(power < 0.80) {
power <- pwr.anova.test(k=3, n=size, f= 0.1632993, sig.level=0.05)$power plot_df <- rbind(plot_df, data.frame("n" = size, "power" = power)) size <- size + 2 print(power) }  Simulated power analysis set.seed(1001) run_sim <- function() { # generate all data create_sim_data <- function(i) { #stdev bias correction c4 <- (sqrt(2/(group_size[1] - 1))) * (gamma(group_size[1]/2)/gamma((group_size[1] - 1)/2)) sds2 <- sds / c4 # pre-allocate matrix test_matrix <- matrix(nrow=sims, ncol=sum(group_size)) # nested loops to create simulated data for all runs for(j in 1:sims) { for(i in 1:length(group_size)) { # col_start & cold_end is used to have the different groups on the same row col_start <- sum(group_size[1:i])-(group_size[i]-1) col_end <- cumsum(group_size)[i] # generate data with rnorm test_matrix[j,col_start:col_end] <- rnorm(group_size[i], mean = means[i], sd = sds2[i]) } } return(test_matrix) } # extract results from simulations get_power <- function() { sig <- rep(NA, sims) eta_2 <- rep(NA, sims) omega_2 <- rep(NA, sims) for(i in 1:sims) { # perform ANOVA on data result <- summary(aov(test_matrix[i,] ~ group)) # calculate effect size eta_2[i] <- result[[1]]$'Sum Sq'[1] / sum(result[[1]]$'Sum Sq') omega_2[i] <- (result[[1]]$'Sum Sq'[1] - (result[[1]]$Df[1] * result[[1]]$'Mean Sq'[2])) / (sum(result[[1]]$'Sum Sq') + result[[1]]$'Mean Sq'[2])
# get p-value from ANOVA
sig_result <- result[[1]]$'Pr(>F)'[1] # check sig.level sig[i] <- sig_result < 0.05 } out <- list("power" = mean(sig), "eta_2" = mean(eta_2), "omega" = mean(omega_2)) } power <- 0 plot_df <- data.frame() eta <- NULL omega <- NULL # repeat the simulation until the desired power is found while(power < 0.8) { # regenerate grouping as group_size increases group <- c(rep(1, group_size[1]), rep(2,group_size[2]), rep(3,group_size[3])) # create data matrix test_matrix <- create_sim_data() # get anova power result <- get_power() # extract power value power <- result$power
# save eta-squared each iteration
eta <- rbind(eta, result$eta_2) # save eta-squared each iteration omega <- rbind(omega, result$omega)
cat("power =", power, "group size =", group_size,"\n\n")
# save group size and power for each iteration
plot_df <- rbind(plot_df, data.frame("group_n" = group_size[1], "power" = power))
# increase group size with 2
group_size <- group_size + 2
}

out <- list("power" = plot_df, "f" = sqrt(eta / (1 - eta)), "omega" = omega)
return(out)
}
sims <- 1000
sim <- run_sim()


This will generate a difference of about 10 % between the two methods (the short line is from simulation)

My thoughts are that the difference is due to Cohen's f being an biased estimator of the population effect size. But how should I interpret my results, are my simulations overestimating the power? If so, how can I get it to match the output from the analytical power estimation.

To summarize my question: why doesn't the two methods give the same output when fed with the same means and standard deviations?

I'd be glad for any pointers were I wen't wrong. Thanks in advance!

• You should comment your code some, get across what you really wanted to accomplish because it seems inordinately complicated. – John May 21 '12 at 16:32
• My first thought is to check whether or not you are comparing apples to apples. Power depends on the alternative you look at. Ar you matching alternative hypotheses in the simulation with the analytical? – Michael R. Chernick May 21 '12 at 17:09
• I'm using the same alternative hypotheses, and I've checked so that simulations will actually converge towards the specified means and standard deviations. – celt-Ail May 21 '12 at 17:20

It seems like you simply hit a specific R peculiarity: When you analyze linear models, and you have a predictor with numerical values, you have to tell R whether it really represents data from a numerical variable (default, leading to a regression model), or whether it actually is a factor (leading to an ANOVA).

In your case you just have to change result <- summary(aov(test_matrix[i,] ~ group)) to result <- summary(aov(test_matrix[i,] ~ factor(group))) to get close to correct results.

In addition, I don't understand your correction to the standard deviations. The sds are the true standard deviations that are required to simulate data with rnorm(). Leave out your correction, and you get even closer to the correct result.

When you have time to explore R some more, you might want to look at some features / strategies that make simulations like yours somewhat easier. E.g.

• rnorm() is vectorized, supply a vector of $\mu$s and $\sigma$s, each with length = number of simulated values, and you can eliminate the double loop in create_sim_data()
• anova(lm()) returns a data frame that is a lot easier to index than the result of summary(aov())
• rep() accepts a vector for its times argument that simplifies its use for your purpose.

Here's your simulation stripped to the bare bones, just for the group size of 40, giving us the p-values.

Nj     <- c(40, 40, 40)                 # group sizes for 3 groups
mu     <- c(0.2, 0, -0.2)               # expected values in groups
sigma  <- c(1, 1, 1)                    # true standard deviations in groups
mus    <- rep(mu, times=Nj)             # for use in rnorm(): vector of mus
sigmas <- rep(sigma, times=Nj)          # for use in rnorm(): vector of true sds
IV     <- factor(rep(1:3, times=Nj))    # factor for ANOVA
nsims  <- 1000                          # number of simulations

# reference: correct power
power.anova.test(groups=3, n=Nj[1], between.var=var(mu), within.var=sigma[1]^2)\$power

doSim <- function() {                   # function to run one ANOVA on simulated data
DV <- rnorm(sum(Nj), mus, sigmas)   # data from all three groups
anova(lm(DV ~ IV))["IV", "Pr(>F)"]  # p-value from ANOVA
}

pVals  <- replicate(nsims, doSim())     # run the simulation nsims times
(power <- sum(pVals < 0.05) / nsims)    # fraction of significant ANOVAs

• Many thanks! I had a feeling the answer would be something simple like this. And regarding the correction to the standard deviations, it grew out of desperation since it did produce more correct results in the original script. Also thanks for the addition about rnorm() being vectorized, that will be really useful! – celt-Ail May 22 '12 at 5:45