> (ybar <- mean(a_acid))
[1] 579.61
> (ybar.trt <- as.vector(tapply(a_acid,variety,mean)))
 [1] 709.9 523.9 659.9 463.5 844.8 590.2 411.8 510.9 537.3 543.9
> (var.trt <- as.vector(tapply(a_acid,variety,var)))
 [1]  662.76667 1436.54444  101.21111  225.61111  229.28889 1049.28889
 [7] 2394.17778   98.76667  610.67778   91.87778
> (rep.trt <- as.vector(tapply(a_acid,variety,length)))
 [1] 10 10 10 10 10 10 10 10 10 10
> (num.trt <- length(ybar.trt))
[1] 10
> 
> sse <- 0;    sstr <- 0
> 
> for (i1 in 1:num.trt) {
+ sse <- sse + (rep.trt[i1]-1)*var.trt[i1]
+ sstr <- sstr + rep.trt[i1]*((ybar.trt[i1]-ybar)^2)
+ }
> 
> df.trt <- num.trt-1
> df.err <- sum(rep.trt)-num.trt
> 
> (mse <- sse/df.err)
[1] 690.0211
> (mstr <- sstr/df.trt)
[1] 162657.3
> 
> (f_star <- mstr/mse)
[1] 235.7280
> (p_f_star <- 1-pf(f_star,df.trt,df.err))
[1] 0
> 
> ### 95% CI for mu_dot
> 
> (ybar + c(qt(.025,df.trt)*sqrt(mstr/sum(rep.trt)),qt(.975,df.trt)*sqrt(mstr/sum(rep.trt))))
[1] 488.3754 670.8446
> 
> #### 95% CI for rho_I
> 
> L <- (1/rep.trt[1])*((mstr/mse)*(1/qf(.975,df.trt,df.err))-1)
> U <- (1/rep.trt[1])*((mstr/mse)*(1/qf(.025,df.trt,df.err))-1)
> 
> (c(L/(1+L),U/(1+U)))
[1] 0.9117866 0.9877013
> 
> #### 95% CI for sigma^2
> 
> (c(sse/qchisq(.975,df.err),sse/qchisq(.025,df.err)))
[1] 525.6819 946.0030
> 
> #### 95% CI for sigma_mu^2
> 
> s2.mu <- (mstr-mse)/rep.trt[1]
> 
> 
> df.num <- (rep.trt[1]*s2.mu)^2
> df.den <- (mstr^2/df.trt) + (mse^2/df.err)
> (df <- round(df.num/df.den))
[1] 9
> 
> (c(df*s2.mu/qchisq(.975,df),df*s2.mu/qchisq(.025,df)))
[1]  7662.953 53981.313