> 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 -- Restricted Model > > (s2b.r <- (msb/(a*n)) - (mse/(a*n))) [1] 26.62048 > > (s2b_df.r <- (s2b.r)^2/(((msb/(a*n))^2/(b-1)) + ((mse/(a*n))^2/(a*b*(n-1))))) [1] 6.780392 > > (c(s2b_df.r*s2b.r/qchisq(.975,s2b_df.r),s2b_df.r*s2b.r/qchisq(.025,s2b_df.r))) [1] 11.5162 113.7333 > > #### 95% CI for sigma_beta^2 -- Unrestricted Model > > (s2b.u <- (msb/(a*n)) - (msab/(a*n))) [1] 23.58456 > > (s2b_df.u <- (s2b.u)^2/(((msb/(a*n))^2/(b-1)) + ((msab/(a*n))^2/((a-1)*(b-1))))) [1] 5.082424 > > (c(s2b_df.u*s2b.u/qchisq(.975,s2b_df.u),s2b_df.u*s2b.u/qchisq(.025,s2b_df.u))) [1] 9.243141 138.951478 > > > #### Tukey's HSD and Bonferroni MSD for comparing all pairs of Oil Patterns > > q_crit <- qtukey(.95,a,(a-1)*(b-1)) > > se_trt_diff <- sqrt((2*msab)/(b*n)) > > (tukey_hsd <- (q_crit/sqrt(2))*se_trt_diff) [1] 7.473484 > > (bon_msd <- qt(1-(0.05/(2*a*(a-1)/2)),(a-1)*(b-1))*se_trt_diff) [1] 7.736271 > > 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 -1.954436 [2,] 1 2 12.992532 i1 i2 [1,] 1 2 -2.217224 [2,] 1 2 13.255319 i1 i2 [1,] 1 3 1.221754 [2,] 1 3 16.168722 i1 i2 [1,] 1 3 0.9589666 [2,] 1 3 16.4315096 i1 i2 [1,] 1 4 -4.773484 [2,] 1 4 10.173484 i1 i2 [1,] 1 4 -5.036271 [2,] 1 4 10.436271 i1 i2 [1,] 2 3 -4.297293 [2,] 2 3 10.649674 i1 i2 [1,] 2 3 -4.560081 [2,] 2 3 10.912462 i1 i2 [1,] 2 4 -10.292532 [2,] 2 4 4.654436 i1 i2 [1,] 2 4 -10.555319 [2,] 2 4 4.917224 i1 i2 [1,] 3 4 -13.468722 [2,] 3 4 1.478246 i1 i2 [1,] 3 4 -13.731510 [2,] 3 4 1.741033 > > ##### 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 > > > # install.packages("lmerTest") > library(lmerTest) Loading required package: lme4 Loading required package: Matrix Attaching package: ‘lmerTest’ The following object is masked from ‘package:lme4’: lmer The following object is masked from ‘package:stats’: step > > wpba.mod2 <- lmer(score ~ oil + (1|bowler) + (1|oil:bowler)) > summary(wpba.mod2) Linear mixed model fit by REML. t-tests use Satterthwaite's method [ lmerModLmerTest] Formula: score ~ oil + (1 | bowler) + (1 | oil:bowler) REML criterion at convergence: 7834.5 Scaled residuals: Min 1Q Median 3Q Max -2.6911 -0.6694 -0.0191 0.6234 3.2689 Random effects: Groups Name Variance Std.Dev. oil:bowler (Intercept) 12.14 3.485 bowler (Intercept) 23.58 4.856 Residual 649.59 25.487 Number of obs: 840, groups: oil:bowler, 60; bowler, 15 Fixed effects: Estimate Std. Error df t value Pr(>|t|) (Intercept) 210.090 2.340 44.893 89.786 < 2e-16 *** oilChameleon -5.519 2.794 42.000 -1.975 0.05482 . oilScorpion -8.695 2.794 42.000 -3.112 0.00333 ** oilShark -2.700 2.794 42.000 -0.966 0.33938 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation of Fixed Effects: (Intr) olChml olScrp oilChamelen -0.597 oilScorpion -0.597 0.500 oilShark -0.597 0.500 0.500 > anova(wpba.mod2) Type III Analysis of Variance Table with Satterthwaite's method Sum Sq Mean Sq NumDF DenDF F value Pr(>F) oil 6962.8 2320.9 3 42 3.5729 0.02172 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > difflsmeans(wpba.mod2) Least Squares Means table: Estimate Std. Error df t value lower upper oilViper - oilChameleon 5.51905 2.79387 42 1.9754 -0.11921 11.15731 oilViper - oilScorpion 8.69524 2.79387 42 3.1123 3.05698 14.33350 oilViper - oilShark 2.70000 2.79387 42 0.9664 -2.93826 8.33826 oilChameleon - oilScorpion 3.17619 2.79387 42 1.1368 -2.46207 8.81445 oilChameleon - oilShark -2.81905 2.79387 42 -1.0090 -8.45731 2.81921 oilScorpion - oilShark -5.99524 2.79387 42 -2.1459 -11.63350 -0.35698 Pr(>|t|) oilViper - oilChameleon 0.054817 . oilViper - oilScorpion 0.003335 ** oilViper - oilShark 0.339375 oilChameleon - oilScorpion 0.262049 oilChameleon - oilShark 0.318747 oilScorpion - oilShark 0.037710 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Confidence level: 95% Degrees of freedom method: Satterthwaite