# The R code used in our study for large effect magnitudes and maximum variability
# is presented below. Note that this code can be easily adapted for other conditions by
# changing the value of muvec.

muvec <- c(0.00,0.00,1.20,1.20)
meanmu <- mean(muvec)
sigb <- sum((muvec-meanmu)^2)/4
eta2p <- sigb/(sigb+1)
k <- length(muvec)
nsim <- 100
njs <- c(10,20,30,40,50,60,70,80,90,100)
BIASmat <- matrix(NA,nrow=length(njs),ncol=3)
rownames(BIASmat) <- njs
colnames(BIASmat) <- c("eta2","epsilon2","omega2")
RMSEmat <- matrix(NA,nrow=length(njs),ncol=3)
rownames(RMSEmat) <- njs
colnames(RMSEmat) <- c("eta2","epsilon2","omega2")
SDmat <- matrix(NA,nrow=length(njs),ncol=3)
rownames(SDmat) <- njs
colnames(SDmat) <- c("eta2","epsilon2","omega2")
niter <- 1
for (nj in njs){
	x <- matrix(NA,nrow=nj,ncol=4)
	eta2 <- rep(NA,nsim)
	epsilon2 <- rep(NA,nsim)
	omega2 <- rep(NA,nsim)
	for (ii in 1: nsim){
		y <- c(rnorm(n=nj,mean=muvec[1],sd=1),
		rnorm(n=nj,mean=muvec[2],sd=1),
		rnorm(n=nj,mean=muvec[3],sd=1),
		rnorm(n=nj,mean=muvec[4],sd=1))
		x <- as.factor(c(rep("mu1",nj),rep("mu2",nj),
		rep("mu3",nj),rep("mu4",nj)))
		res <- anova(aov(y~x))
		res <- as.matrix(res)
		SSb <- res[1,2]
		SSt <- res[1,2] + res[2,2]
		MSw <- res[2,3]
		eta2[ii] <- SSb/SSt
		epsilon2[ii] <- (SSb - 3*MSw)/SSt
		omega2[ii] <- (SSb - 3*MSw)/(SSt+MSw)
	}
	BIASmat[niter,1] <- mean(eta2) - eta2p
	BIASmat[niter,2] <- mean(epsilon2) - eta2p
	BIASmat[niter,3] <- mean(omega2) - eta2p
	RMSEmat[niter,1] <- sqrt(sum((eta2-eta2p)^2)/nsim)
	RMSEmat[niter,2] <- sqrt(sum((epsilon2-eta2p)^2)/nsim)
	RMSEmat[niter,3] <- sqrt(sum((omega2-eta2p)^2)/nsim)
	SDmat[niter,1] <- sqrt(sum((eta2-mean(eta2))^2)/nsim)
	SDmat[niter,2] <- sqrt(sum((epsilon2-mean(epsilon2))^2)
/nsim)
	SDmat[niter,3] <- sqrt(sum((omega2-mean(omega2))^2)
/nsim)
	niter <- niter+1
}