pdf("E:\\sta4211\\KNNL_TA19_07.pdf") castle <- read.table("http://www.stat.ufl.edu/~winner/sta4211/datasets/CH19TA07.txt",header=F, col.names=c("sales","trt_a","trt_b","obsnum")) attach(castle) trt_a <- factor(trt_a) trt_b <- factor(trt_b) n <- max(obsnum) (cell_mean <- as.vector(tapply(sales,list(trt_a,trt_b),mean))) (a_mean <- as.vector(tapply(sales,trt_a,mean))) a_lev <- length(a_mean) (b_mean <- as.vector(tapply(sales,trt_b,mean))) b_lev <- length(b_mean) ab_lev <- a_lev*b_lev trt_ab <- rep(1:ab_lev,each=n) trt_ab <- factor(trt_ab) (all_mean <- mean(sales)) castle.aov1 <- aov(sales ~ trt_ab) summary(castle.aov1) castle.aov2 <- aov(sales ~ trt_a + trt_b + trt_a:trt_b) summary(castle.aov2) interaction.plot(trt_a,trt_b,sales, ylim=c(0.7*min(sales),1.3*max(sales))) interaction.plot(trt_b,trt_a,sales, ylim=c(0.7*min(sales),1.3*max(sales))) e <- residuals(castle.aov2) yhat <- predict(castle.aov2) sse <- deviance(castle.aov2) df.e <- df.residual(castle.aov2) mse <- sse/df.e plot(yhat,e,xlab="Yht",ylab="Residual") qqnorm(e); qqline(e) ##### Comparison of Main Effects for Factors A,B (separate families) n <- max(obsnum) TukeyHSD(castle.aov2,"trt_a") TukeyHSD(castle.aov2,"trt_b") Bon_MSD_A <- qt(1-(0.025/(a_lev*(a_lev-1)/2)),df.e)*sqrt(2*mse/(b_lev*n)) Bon_MSD_B <- qt(1-(0.025/(b_lev*(b_lev-1)/2)),df.e)*sqrt(2*mse/(a_lev*n)) for(i1 in 1:(a_lev-1)) { for(i2 in (i1+1):a_lev) { print(cbind(i1,i2,(a_mean[i1]-a_mean[i2])-Bon_MSD_A,(a_mean[i1]-a_mean[i2])+Bon_MSD_A)) }} for(i1 in 1:(b_lev-1)) { for(i2 in (i1+1):b_lev) { print(cbind(i1,i2,(b_mean[i1]-b_mean[i2])-Bon_MSD_B,(b_mean[i1]-b_mean[i2])+Bon_MSD_B)) }} ### Comparison of Treatment Means - Interaction Model castle.aov3 <- aov(sales ~ trt_ab) summary(castle.aov3) TukeyHSD(castle.aov3,"trt_ab") Bon_MSD_AB <- qt(1-(0.025/(ab_lev*(ab_lev-1)/2)),df.e)*sqrt(2*mse/n) Schef_MSD_AB <- sqrt(ab_lev-1*qf(0.95,ab_lev-1,df.e))*sqrt(2*mse/n) for(i1 in 1:(ab_lev-1)) { for(i2 in (i1+1):ab_lev) { print(cbind(i1,i2,(cell_mean[i1]-cell_mean[i2])-Bon_MSD_AB, (cell_mean[i1]-cell_mean[i2])+Bon_MSD_AB)) }} for(i1 in 1:(ab_lev-1)) { for(i2 in (i1+1):ab_lev) { print(cbind(i1,i2,(cell_mean[i1]-cell_mean[i2])-Schef_MSD_AB, (cell_mean[i1]-cell_mean[i2])+Schef_MSD_AB)) }} dev.off()