#pdf("C:\\Rmisc\\graphs\\kuehl_05_01a.pdf") casting <- c(rep(1:3,each=10)) tstrength <- c(88.0,88.0,94.8,90.0,93.0,89.0,86.0,92.9,89.0,93.0, 85.9,88.6,90.0,87.1,85.6,86.0,91.0,89.6,93.0,87.5, 94.2,91.5,92.0,96.5,95.6,93.8,92.5,93.2,96.2,92.5) alloy <- data.frame(casting,tstrength) attach(alloy) casting <- factor(casting) c_mean <- tapply(tstrength,casting,mean) c_var <- tapply(tstrength,casting,var) c_n <- tapply(tstrength,casting,length) t <- length(c_n) r <- c_n[[1]] N <- sum(c_n) 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) SSW <- sum((c_n-1)*c_var); MSW <- SSW/(N-t) print(c("SSW=",round(SSW,4),"MSW=",round(MSW,4)),quote=F) sigma2_a <- (MSA-MSW)/r sigma2_e <- MSW print(c("Point Estimate for sigma2_a:", round(sigma2_a,4)),quote=F) print(c("Point Estimate for sigma2_e:", round(sigma2_e,4)),quote=F) 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) 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) 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) rho_ic <- sigma2_a/(sigma2_a+sigma2_e) print(c("Point Estimate for rho_ic:", round(rho_ic,4)),quote=F) 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) #dev.off()