> 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))) [1] 45 65 40 43 69 44 > > (a_mean <- as.vector(tapply(sales,trt_a,mean))) [1] 44 67 42 > a_lev <- length(a_mean) > > (b_mean <- as.vector(tapply(sales,trt_b,mean))) [1] 50 52 > 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)) [1] 51 > > castle.aov1 <- aov(sales ~ trt_ab) > summary(castle.aov1) Df Sum Sq Mean Sq F value Pr(>F) trt_ab 5 1580 316.000 30.581 0.0003384 *** Residuals 6 62 10.333 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > castle.aov2 <- aov(sales ~ trt_a + trt_b + trt_a:trt_b) > summary(castle.aov2) Df Sum Sq Mean Sq F value Pr(>F) trt_a 2 1544 772.00 74.7097 5.754e-05 *** trt_b 1 12 12.00 1.1613 0.3226 trt_a:trt_b 2 24 12.00 1.1613 0.3747 Residuals 6 62 10.33 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > interaction.plot(trt_a,sales,trt_b, ylim=c(0.8*min(sales),1.2*max(sales))) There were 12 warnings (use warnings() to see them) > interaction.plot(trt_b,sales,trt_a, ylim=c(0.8*min(sales),1.2*max(sales))) There were 11 warnings (use warnings() to see them) > > 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") Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = sales ~ trt_a + trt_b + trt_a:trt_b) $trt_a diff lwr upr p adj 2-1 23 16.02572 29.974281 0.0001335 3-1 -2 -8.97428 4.974281 0.6714131 3-2 -25 -31.97428 -18.025719 0.0000829 > TukeyHSD(castle.aov2,"trt_b") Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = sales ~ trt_a + trt_b + trt_a:trt_b) $trt_b diff lwr upr p adj 2-1 2 -2.541276 6.541276 0.3226055 > > 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)) + }} i1 i2 [1,] 1 2 -30.47249 -15.52751 i1 i2 [1,] 1 3 -5.472485 9.472485 i1 i2 [1,] 2 3 17.52751 32.47249 > > 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)) + }} i1 i2 [1,] 1 2 -6.541276 2.541276 > > ### Comparison of Treatment Means - Interaction Model > > castle.aov3 <- aov(sales ~ trt_ab) > summary(castle.aov3) Df Sum Sq Mean Sq F value Pr(>F) trt_ab 5 1580 316.000 30.581 0.0003384 *** Residuals 6 62 10.333 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > TukeyHSD(castle.aov3,"trt_ab") Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = sales ~ trt_ab) $trt_ab diff lwr upr p adj 2-1 -2 -14.793417 10.793417 0.9848729 3-1 20 7.206583 32.793417 0.0060698 4-1 24 11.206583 36.793417 0.0023243 5-1 -5 -17.793417 7.793417 0.6488665 6-1 -1 -13.793417 11.793417 0.9993784 3-2 22 9.206583 34.793417 0.0036950 4-2 26 13.206583 38.793417 0.0015059 5-2 -3 -15.793417 9.793417 0.9236156 6-2 1 -11.793417 13.793417 0.9993784 4-3 4 -8.793417 16.793417 0.8034156 5-3 -25 -37.793417 -12.206583 0.0018644 6-3 -21 -33.793417 -8.206583 0.0047154 5-4 -29 -41.793417 -16.206583 0.0008246 6-4 -25 -37.793417 -12.206583 0.0018644 6-5 4 -8.793417 16.793417 0.8034156 > > 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)) + }} i1 i2 [1,] 1 2 -35.10171 -4.898292 i1 i2 [1,] 1 3 -10.10171 20.10171 i1 i2 [1,] 1 4 -13.10171 17.10171 i1 i2 [1,] 1 5 -39.10171 -8.898292 i1 i2 [1,] 1 6 -14.10171 16.10171 i1 i2 [1,] 2 3 9.898292 40.10171 i1 i2 [1,] 2 4 6.898292 37.10171 i1 i2 [1,] 2 5 -19.10171 11.10171 i1 i2 [1,] 2 6 5.898292 36.10171 i1 i2 [1,] 3 4 -18.10171 12.10171 i1 i2 [1,] 3 5 -44.10171 -13.89829 i1 i2 [1,] 3 6 -19.10171 11.10171 i1 i2 [1,] 4 5 -41.10171 -10.89829 i1 i2 [1,] 4 6 -16.10171 14.10171 i1 i2 [1,] 5 6 9.898292 40.10171 > > > 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)) + }} i1 i2 [1,] 1 2 -24.08213 -15.91787 i1 i2 [1,] 1 3 0.9178682 9.082132 i1 i2 [1,] 1 4 -2.082132 6.082132 i1 i2 [1,] 1 5 -28.08213 -19.91787 i1 i2 [1,] 1 6 -3.082132 5.082132 i1 i2 [1,] 2 3 20.91787 29.08213 i1 i2 [1,] 2 4 17.91787 26.08213 i1 i2 [1,] 2 5 -8.082132 0.0821318 i1 i2 [1,] 2 6 16.91787 25.08213 i1 i2 [1,] 3 4 -7.082132 1.082132 i1 i2 [1,] 3 5 -33.08213 -24.91787 i1 i2 [1,] 3 6 -8.082132 0.0821318 i1 i2 [1,] 4 5 -30.08213 -21.91787 i1 i2 [1,] 4 6 -5.082132 3.082132 i1 i2 [1,] 5 6 20.91787 29.08213 > >