My question is based on This Paper (see its Appendix for R code). The paper is trying to basically obtain the empirical distribution of a statistic called $\eta^2$ (eta-squared) via large replications. If we assume a one-way ANOVA with a normal population for each level (i.e., group), the population parameter of this statistic ($\eta^2$) is defined as:

$\eta^2 = \frac{\sigma_{\text{between normal levels' means}}^2}{\sigma_{\text{between normal levels' means}}^2 + ~\sigma_{\text{within each normal level}}^2}$ (see pop.eta.sq below)

where $\sigma_{\text{between normal levels' means}}^2$ is the variance between levels means, and $\sigma_{\text{within each normal level}}^2$ is the variance common to all levels (I set this to $1$).


Based on "page 137" of This Paper, as sample size increases the expected value (the filled black circles in picture below) of the replications should get closer to the pop.eta.sq. But why I see the opposite in my simulation below in R (see picture below code): is my simulation legitimate?

Anova = function(n, means = c(0, .4, .8, 1.2), n.sim = 1e3){
    k = length(means)
 sigb = var(means)*(k-1)/k # "*(k-1)/k" changes sample variance to pop. variance 
pop.eta.sq = round(sigb / (sigb + 1), 3) # I have set Sigma.within.each.normal.level = 1

eta.sim = Vectorize(function(n){
    eta = numeric(n.sim)
for(i in 1:n.sim){ 
     y = as.vector(unlist(mapply(FUN = rnorm, n = rep(n, k), mean = means)))
groups = gl(k, 1, k*n)
   sim = anova(aov(y ~ groups))
eta[i] = sim[, 2][1] / (sim[, 2][1] + sim[, 2][2]) # sample estimate of eta.sq
  hist(eta, main = NA)
  mtext(paste0("Sample Size = ", n, "\n", "Pop.eta^2 = ", pop.eta.sq), font = 2)
  points(mean(eta), 0, pch = 19, xpd = NA, cex = 2)
  text(mean(eta), 0, round(mean(eta), 3), font = 2, pos = 3, xpd = NA)
  }, "n")
  par(mfrow = c(2, length(n)))
Anova(n = c(5, 50)) # Two different group sample sizes, 5 & 50

enter image description here


You have an bug in your code; your group variable is not correct. Replace that line with this:

groups = factor(rep(1:k,each = n))

and it should work properly.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.