s_a <- 10 s_e <- 5 t <- 3 r <- 6 mu <- rep(100,t*r) numsim <- 1 N <- t*r F.alpha <- qf(.95,t-1,N-t) reject.H0 <- 0 rho.I <- rep(0,numsim) for (i in 1:numsim) { a <- rep(s_a*rnorm(t),each=r) e <- s_e*rnorm(t*r) y <- mu+a+e y1 <- y[1:r] y2 <- y[(r+1):(2*r)] y3 <- y[(2*r+1):(3*r)] ybar.1 <- mean(y1); var.1 <- var(y1) ybar.2 <- mean(y2); var.2 <- var(y2) ybar.3 <- mean(y3); var.3 <- var(y3) ybar.0 <- (ybar.1+ybar.2+ybar.3)/3 SS.TRT <- r * ( (ybar.1-ybar.0)^2 + (ybar.2-ybar.0)^2 + (ybar.3-ybar.0)^2 ) SS.ERR <- (r-1) * (var.1 + var.2 + var.3) F.obs <- (SS.TRT/(t-1))/(SS.ERR/(N-t)) if (F.obs >= F.alpha) reject.H0 <- reject.H0 + 1 s2_e.est <- SS.ERR/(N-t) s2_a.est <- ((SS.TRT/(t-1))-s2_e.est)/r rho.I[i] <- s2_a.est/(s2_a.est + s2_e.est) } print(c("Empirical Rejection Rate:", reject.H0/numsim),quote=F) mu a e y