> SSA <- sum(c_n*(c_mean-mean(tstrength))^2); MSA <- SSA/(t-1)
> print(c("SSA=",round(SSA,4),"MSA=",round(MSA,4)),quote=F)
[1] SSA=     147.8847 MSA=     73.9423 
> 
> SSW <- sum((c_n-1)*c_var); MSW <- SSW/(N-t)
> print(c("SSW=",round(SSW,4),"MSW=",round(MSW,4)),quote=F)
[1] SSW=    157.102 MSW=    5.8186 
> 
> sigma2_a <- (MSA-MSW)/r
> sigma2_e <- MSW
> 
> print(c("Point Estimate for sigma2_a:", round(sigma2_a,4)),quote=F)
[1] Point Estimate for sigma2_a: 6.8124                      
> print(c("Point Estimate for sigma2_e:", round(sigma2_e,4)),quote=F)
[1] Point Estimate for sigma2_e: 5.8186                      
> 
> 
> F_0 <- MSA/MSW; p_F_0 <- 1-pf(F_0,t-1,N-t)
> print(c("F_0=",round(F_0,4),"P-value=",round(p_F_0,4)),quote=F)
[1] F_0=     12.7079  P-value= 1e-04   
> 
> A <- qchisq(.95,N-t); B <- qchisq(.05,N-t)
> 
> sigma2_e_lo <- SSW/A;  sigma2_e_hi <- SSW/B
> 
> print(c("90% CI for sigma2_e:", round(sigma2_e_lo,4), round(sigma2_e_hi,4)),
+  quote=F)
[1] 90% CI for sigma2_e: 3.9165               9.7268              
> 
> C <- qchisq(.975,t-1);  D <- qchisq(.025,t-1)
> F_u <- qf(.975,t-1,N-t);   F_l <- qf(.025,t-1,N-t)
> 
> sigma2_a_lo <- SSA*(1-F_u/F_0)/(r*C)
> sigma2_a_hi <- SSA*(1-F_l/F_0)/(r*D)
> 
> print(c("90% CI for sigma2_a:", round(sigma2_a_lo,4), round(sigma2_a_hi,4)),
+  quote=F)
[1] 90% CI for sigma2_a: 1.3353               291.4742            
> 
> rho_ic <- sigma2_a/(sigma2_a+sigma2_e)
> 
> print(c("Point Estimate for rho_ic:", round(rho_ic,4)),quote=F)
[1] Point Estimate for rho_ic: 0.5393                    
> 
> F_u <- qf(.95,t-1,N-t);   F_l <- qf(.05,t-1,N-t)
> 
> rho_ic_lo <- (F_0-F_u)/(F_0+(r-1)*F_u)
> rho_ic_hi <- (F_0-F_l)/(F_0+(r-1)*F_l)
> 
> print(c("90% CI for rho_ic:", round(rho_ic_lo,4), round(rho_ic_hi,4)),
+  quote=F)
[1] 90% CI for rho_ic: 0.2181             0.961             
>