> 
> (ybar <- mean(rating))
[1] 71.45
> (ybar.trt <- as.vector(tapply(rating,officer,mean)))
[1] 75.00 70.50 54.75 79.75 77.25
> (var.trt <- as.vector(tapply(rating,officer,var)))
[1] 67.33333 91.66667 72.25000 74.25000 60.91667
> (rep.trt <- as.vector(tapply(rating,officer,length)))
[1] 4 4 4 4 4
> (num.trt <- length(ybar.trt))
[1] 5
> 
> 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] 73.28333
> (mstr <- sstr/df.trt)
[1] 394.925
> 
> (f_star <- mstr/mse)
[1] 5.389015
> (p_f_star <- 1-pf(f_star,df.trt,df.err))
[1] 0.006802857
> 
> ### 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] 59.11238 83.78762
> 
> #### 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.09431972 0.91943643
> 
> #### 95% CI for sigma^2
> 
> (c(sse/qchisq(.975,df.err),sse/qchisq(.025,df.err)))
[1]  39.98961 175.53909
> 
> #### 90% 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] 3
> 
> (c(df*s2.mu/qchisq(.95,df),df*s2.mu/qchisq(.05,df)))
[1]  30.8688 685.6154
> 
> stripchart(rating ~ officer, xlab="Rating", ylab="Officer")
>