> wpba.mod1 <- aov(score ~ oil*bowler) > anova(wpba.mod1) Analysis of Variance Table Response: score Df Sum Sq Mean Sq F value Pr(>F) oil 3 8785 2928.37 4.5080 0.003817 ** bowler 14 29965 2140.34 3.2949 3.844e-05 *** oil:bowler 42 34423 819.60 1.2617 0.127108 Residuals 780 506679 649.59 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > ##### Direct computation of ANOVA > > (ybar <- mean(score)) [1] 205.8619 > (ybar.bowler <- as.vector(tapply(score,bowler,mean))) [1] 209.6607 198.8750 212.9464 216.1786 208.1964 198.6250 209.5357 208.4286 [9] 198.9643 209.3036 202.4286 194.7500 202.7321 211.5893 205.7143 > (ybar.oil <- as.vector(tapply(score,oil,mean))) [1] 210.0905 204.5714 201.3952 207.3905 > (ybar.trt <- as.vector(tapply(score,list(bowler,oil),mean))) [1] 223.2857 211.7857 218.1429 219.4286 210.5714 211.0000 223.3571 209.5714 [9] 199.5714 205.8571 202.5000 206.2143 198.5000 212.0000 199.5714 208.7857 [17] 195.7143 208.5000 216.6429 198.4286 203.2857 199.2857 214.2143 198.5714 [25] 213.7143 205.2857 182.6429 207.8571 205.8571 209.7857 195.5714 196.4286 [33] 209.6429 212.4286 204.7143 193.0714 194.4286 208.6429 193.2857 198.3571 [41] 194.3571 196.1429 210.7143 208.2857 204.8571 211.0000 191.5714 215.5000 [49] 216.2143 219.0714 187.1429 221.0714 201.2857 204.4286 219.2857 207.5714 [57] 194.0000 193.8571 220.2143 208.6429 > (var.trt <- as.vector(tapply(score,list(bowler,oil),var))) [1] 1800.8352 827.8736 600.5934 350.2637 431.4945 955.3846 1173.6319 [8] 633.3407 783.0330 1090.1319 722.5769 1022.9505 240.4231 901.6923 [15] 783.0330 528.1813 732.0659 688.2692 472.7088 353.4945 285.6044 [22] 917.2967 1000.7967 469.4945 429.7582 209.9121 665.6319 501.2088 [29] 649.5165 572.7967 743.4945 590.8791 771.6319 757.3407 823.4505 [36] 562.8407 409.8022 431.0165 344.8352 975.0165 457.9396 279.9780 [43] 669.7582 503.1429 410.7473 832.3077 391.0330 198.5769 1269.1044 [50] 615.9176 950.1319 521.7637 672.8352 717.1868 774.9890 312.1099 [57] 707.2308 1131.5165 169.2582 185.4780 > (rep.trt <- as.vector(tapply(score,list(bowler,oil),length))) [1] 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 [26] 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 14 [51] 14 14 14 14 14 14 14 14 14 14 > (b <- length(ybar.bowler)) [1] 15 > (a <- length(ybar.oil)) [1] 4 > (n <- rep.trt[1]) [1] 14 > > sse <- 0; sstr <- 0; ssa <- 0; ssb <- 0 > > for (i1 in 1:(a*b)) { + sse <- sse + (n-1)*var.trt[i1] + sstr <- sstr + n*((ybar.trt[i1]-ybar)^2) + } > > for(i1 in 1:a) { + ssa <- ssa + b*n*(ybar.oil[i1]-ybar)^2 + } > > for(i1 in 1:b) { + ssb <- ssb + a*n*(ybar.bowler[i1]-ybar)^2 + } > > ssab <- sstr-ssa-ssb > > df.a <- a-1; df.b <- b-1; df.ab <- (a-1)*(b-1) > df.err <- a*b*(n-1) > > (mse <- sse/df.err) [1] 649.5885 > (msab <- ssab/df.ab) [1] 819.5998 > (msa <- ssa/df.a) [1] 2928.365 > (msb <- ssb/df.b) [1] 2140.335 > > (f_star_ab <- msab/mse) [1] 1.261722 > (p_f_star_ab <- 1-pf(f_star_ab,df.ab,df.err)) [1] 0.1271085 > > (f_star_a <- msa/msab) [1] 3.572921 > (p_f_star_a <- 1-pf(f_star_a,df.a,df.ab)) [1] 0.02171787 > > (f_star_b <- msb/mse) [1] 3.29491 > (p_f_star_b <- 1-pf(f_star_b,df.b,df.err)) [1] 3.844295e-05 > > #### 95% CI for sigma_beta^2 > > (s2b <- (msb/(a*n)) - (mse/(a*n))) [1] 26.62048 > > (s2b_df <- round((s2b)^2/(((msb/(a*n))^2/(b-1)) + ((mse/(a*n))^2/(a*b*(n-1)))))) [1] 7 > > (c(s2b_df*s2b/qchisq(.975,s2b_df),s2b_df*s2b/qchisq(.025,s2b_df))) [1] 11.63718 110.27088 > > #### Tukey's HSD and Bonferroni MSD for comparing all pairs of Oil Patterns > > q_crit <- 3.44 > > se_trt_diff <- sqrt(msab/(b*n)) > > (tukey_hsd <- (q_crit/sqrt(2))*se_trt_diff) [1] 4.805457 > > (bon_msd <- qt(1-(0.05/(2*a*(a-1)/2)),(a-1)*(b-1))*se_trt_diff) [1] 5.47037 > > for(i1 in 1:(a-1)) { + for (i2 in (i1+1):a) { + print(cbind(i1,i2,(ybar.oil[i1]-ybar.oil[i2])+c(-tukey_hsd,tukey_hsd))) + print(cbind(i1,i2,(ybar.oil[i1]-ybar.oil[i2])+c(-bon_msd,bon_msd))) + }} i1 i2 [1,] 1 2 0.7135904 [2,] 1 2 10.3245048 i1 i2 [1,] 1 2 0.04867761 [2,] 1 2 10.98941763 i1 i2 [1,] 1 3 3.889781 [2,] 1 3 13.500695 i1 i2 [1,] 1 3 3.224868 [2,] 1 3 14.165608 i1 i2 [1,] 1 4 -2.105457 [2,] 1 4 7.505457 i1 i2 [1,] 1 4 -2.77037 [2,] 1 4 8.17037 i1 i2 [1,] 2 3 -1.629267 [2,] 2 3 7.981648 i1 i2 [1,] 2 3 -2.29418 [2,] 2 3 8.64656 i1 i2 [1,] 2 4 -7.624505 [2,] 2 4 1.986410 i1 i2 [1,] 2 4 -8.289418 [2,] 2 4 2.651322 i1 i2 [1,] 3 4 -10.800695 [2,] 3 4 -1.189781 i1 i2 [1,] 3 4 -11.4656081 [2,] 3 4 -0.5248681 > > ##### Simultaneous 95% CI's for Oil Pattern Means (Bonferroni) > > c1 <- (a-1)/(n*a*b); c2 <- 1/(n*a*b) > > df_num <- (c1*msab + c2*msb)^2 > > df_den <- ((c1*msab)^2/((a-1)*(b-1))) + ((c2*msb)^2/(b-1)) > > (df.mean <- round(df_num/df_den)) [1] 45 > > for (i1 in 1:a) { + print(cbind(i1,ybar.oil[i1] + + c(qt(.025/(a*(a-1)/2),df.mean)*sqrt(c1*msab+c2*msb), + qt(1-.025/(a*(a-1)/2),df.mean)*sqrt(c1*msab+c2*msb)))) + } i1 [1,] 1 203.6325 [2,] 1 216.5485 i1 [1,] 2 198.1135 [2,] 2 211.0294 i1 [1,] 3 194.9373 [2,] 3 207.8532 i1 [1,] 4 200.9325 [2,] 4 213.8485 >