Hi, I am analyzing a repeated-measures Greco-Latin Square with the aov command. I am using aov to calculate the MSs and then picking by hand the appropriate neumerator and denominator terms for the F tests.
The data are the following: responseFinger mapping.code Subject.n index middle ring little ---------------------------------------------------------------------------- 1 1 green(1) yellow(4) blue(2) red(3) 1 2 green(1) yellow(4) blue(2) red(3) 2 3 red(2) blue(3) yellow(1) green(4) 2 4 red(2) blue(3) yellow(1) green(4) 3 5 yellow(3) green(2) red(4) blue(1) 3 6 yellow(3) green(2) red(4) blue(1) 4 7 blue(4) red(1) green(3) yellow(2) 4 8 blue(4) red(1) green(3) yellow(2) There are 4 factors: factor levels type ----------------------------------------------------------------- responseFinger index, middle, ring, little within-subject stimulusName green, yellow, blue, red within-subject oom 1, 2, 3, 4 within-subject mapping.code 1, 2, 3, 4 between-subjects Subject.n 1, 2, 3, 4, 5, 6, 7, 8 nested within mapping.code DV = asin.Stimulus.ER There are 32 observations and 31 total dfs. I fit the following model: m1 <- asin.Stimulus.ER ~ stimulusName + responseFinger + oom + mapping.code + mapping.code:Subject.n + stimulusName:mapping.code + stimulusName:mapping.code:Subject.n, data=w.tmp) summary(m1) Df Sum Sq Mean Sq stimulusName 3 0.005002 0.001667 responseFinger 3 0.033730 0.011243 oom 3 0.006146 0.002049 mapping.code 3 0.028190 0.009397 mapping.code:Subject.n 4 0.020374 0.005094 stimulusName:mapping.code 3 0.022318 0.007439 stimulusName:mapping.code:Subject.n 12 0.013903 0.001159 There are 3 dfs for each main effect of stimulusName, responseFinger, oom, and mapping.code. That leaves 3 df to estimate any higher-order interactions involving these 4 factors. To capture this variance I use stimulusName:mapping.code, but I believe I could use any of the two-way interactions. The variance aassociated with this effect should be orthogonal to the variance for the main effects. However, when I look at the contrasts with summary.lm it seems as though the estimates of contrasts involving responseFinger, for example, depend on whether stimulusName:mapping.code is in the model. That suggests to me that variance contributing to stimulusName:mapping.code in the ANOVA table is not orthogonal to the main effect of responseFinger, as it should be. Yet, I have calculated the MS by hand and am confident that the SS in the ANOVA table for stimulusName:mapping.code is orthogonal to other terms in the model. If that is true, then why aren't the contrasts that R is using to estimate the SS in the ANOVA table also not orthogonal? m2 <- aov(asin.Stimulus.ER ~ stimulusName + responseFinger + oom + mapping.code + mapping.code:Subject.n + stimulusName:mapping.code, data=w.tmp, contrasts=list(stimulusName="contr.poly", responseFinger="contr.poly", oom="contr.poly",mapping.code="contr.poly", Subject.n="contr.poly")) summary.lm(m2) Call: aov(formula = asin.Stimulus.ER ~ stimulusName + responseFinger + oom + mapping.code + mapping.code:Subject.n, data = w.tmp, contrasts = list(stimulusName = "contr.poly", responseFinger = "contr.poly", oom = "contr.poly", mapping.code = "contr.poly", Subject.n = "contr.poly")) Residuals: Min 1Q Median 3Q Max -0.063950 -0.027306 -0.002311 0.033117 0.055900 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.0894680 0.0086867 10.299 3.38e-08 *** stimulusName.L 0.0095592 0.0173734 0.550 0.59027 stimulusName.Q 0.0205827 0.0173734 1.185 0.25456 stimulusName.C 0.0104971 0.0173734 0.604 0.55474 responseFinger.L 0.0246606 0.0173734 1.419 0.17622 responseFinger.Q -0.0571795 0.0173734 -3.291 0.00495 ** responseFinger.C -0.0184023 0.0173734 -1.059 0.30626 oom.L 0.0243220 0.0173734 1.400 0.18187 oom.Q -0.0126019 0.0173734 -0.725 0.47940 oom.C -0.0042307 0.0173734 -0.244 0.81090 mapping.code.L 0.0465695 0.0173734 2.681 0.01711 * mapping.code.Q -0.0300432 0.0173734 -1.729 0.10428 mapping.code.C 0.0212706 0.0173734 1.224 0.23971 mapping.code13:Subject.n.L -0.0154836 0.0245697 -0.630 0.53805 mapping.code14:Subject.n.L -0.0642852 0.0245697 -2.616 0.01945 * mapping.code15:Subject.n.L -0.0001106 0.0245697 -0.005 0.99647 mapping.code16:Subject.n.L -0.0268560 0.0245697 -1.093 0.29162 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.04914 on 15 degrees of freedom Multiple R-Squared: 0.7207, Adjusted R-squared: 0.4227 F-statistic: 2.419 on 16 and 15 DF, p-value: 0.0474 m3 <- aov(asin.Stimulus.ER ~ stimulusName + responseFinger + oom + mapping.code + mapping.code:Subject.n + stimulusName:mapping.code, data=w.tmp, contrasts=list(stimulusName="contr.poly", responseFinger="contr.poly", oom="contr.poly",mapping.code="contr.poly", Subject.n="contr.poly")) summary.lm(m3) Call: aov(formula = asin.Stimulus.ER ~ stimulusName + responseFinger + oom + mapping.code + mapping.code:Subject.n + stimulusName:mapping.code, data = w.tmp, contrasts = list(stimulusName = "contr.poly", responseFinger = "contr.poly", oom = "contr.poly", mapping.code = "contr.poly", Subject.n = "contr.poly")) Residuals: Min 1Q Median 3Q Max -4.332e-02 -1.370e-02 -9.216e-19 1.370e-02 4.332e-02 Coefficients: (6 not defined because of singularities) Estimate Std. Error t value Pr(>|t|) (Intercept) 0.0894680 0.0060170 14.869 4.3e-09 *** stimulusName.L 0.0095592 0.0120340 0.794 0.442421 stimulusName.Q 0.0205827 0.0120340 1.710 0.112901 stimulusName.C 0.0104971 0.0120340 0.872 0.400168 responseFinger.L 0.0483851 0.0163680 2.956 0.012008 * responseFinger.Q -0.1217915 0.0269089 -4.526 0.000694 *** responseFinger.C -0.0089211 0.0142389 -0.627 0.542703 oom.L 0.0319155 0.0131826 2.421 0.032257 * oom.Q -0.0465613 0.0269089 -1.730 0.109180 oom.C -0.0080275 0.0123312 -0.651 0.527323 mapping.code.L 0.0465695 0.0120340 3.870 0.002229 ** mapping.code.Q -0.0300432 0.0120340 -2.497 0.028094 * mapping.code.C 0.0212706 0.0120340 1.768 0.102532 mapping.code13:Subject.n.L -0.0154836 0.0170187 -0.910 0.380838 mapping.code14:Subject.n.L -0.0642852 0.0170187 -3.777 0.002636 ** mapping.code15:Subject.n.L -0.0001106 0.0170187 -0.006 0.994921 mapping.code16:Subject.n.L -0.0268560 0.0170187 -1.578 0.140543 stimulusName.L:mapping.code.L 0.0848984 0.0601701 1.411 0.183649 stimulusName.Q:mapping.code.L -0.1444769 0.0538178 -2.685 0.019869 * stimulusName.L:mapping.code.Q -0.0853739 0.0269089 -3.173 0.008029 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.03404 on 12 degrees of freedom Multiple R-Squared: 0.8928, Adjusted R-squared: 0.723 F-statistic: 5.259 on 19 and 12 DF, p-value: 0.002629 The dataframe is below. Thanks for any insight you can provide, Steve w.tmp <- structure(list(activeMapping = structure(as.integer(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)), .Label = "letter", class = "factor"), mapping.code = structure(as.integer(c(1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4)), .Label = c("13", "14", "15", "16"), class = "factor"), oom = structure(as.integer(c(1, 1, 4, 4, 2, 2, 3, 3, 2, 2, 3, 3, 1, 1, 4, 4, 3, 3, 2, 2, 4, 4, 1, 1, 4, 4, 1, 1, 3, 3, 2, 2)), .Label = c("1", "2", "3", "4"), class = "factor"), stimulusName = structure(as.integer(c(1, 1, 2, 2, 3, 3, 4, 4, 4, 4, 3, 3, 2, 2, 1, 1, 2, 2, 1, 1, 4, 4, 3, 3, 3, 3, 4, 4, 1, 1, 2, 2)), .Label = c("Q", "J", "X", "B"), class = "factor"), responseFinger = structure(as.integer(c(1, 1, 2, 2, 3, 3, 4, 4, 1, 1, 2, 2, 3, 3, 4, 4, 1, 1, 2, 2, 3, 3, 4, 4, 1, 1, 2, 2, 3, 3, 4, 4)), .Label = c("index", "middle", "ring", "little"), class = "factor"), Subject.a = structure(as.integer(c(1, 5, 1, 5, 1, 5, 1, 5, 2, 6, 2, 6, 2, 6, 2, 6, 3, 7, 3, 7, 3, 7, 3, 7, 4, 8, 4, 8, 4, 8, 4, 8)), .Label = c("17", "18", "19", "20", "21", "22", "23", "24"), class = "factor"), Subject = structure(as.integer(c(1, 5, 1, 5, 1, 5, 1, 5, 2, 6, 2, 6, 2, 6, 2, 6, 3, 7, 3, 7, 3, 7, 3, 7, 4, 8, 4, 8, 4, 8, 4, 8)), .Label = c("25", "26", "27", "28", "29", "30", "36", "41"), class = "factor"), Subject.n = structure(as.integer(c(1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2)), .Label = c("1", "2" ), class = "factor"), Block = structure(as.integer(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)), .Label = "1", class = "factor"), Stimulus.ACC = c(0.9875, 0.995833333333333, 0.9875, 0.975, 0.9625, 0.995833333333333, 0.975, 0.9875, 0.954166666666667, 0.983333333333333, 0.9125, 0.95, 0.908333333333333, 0.983333333333333, 0.958333333333333, 0.979166666666667, 0.9875, 0.979166666666667, 0.9625, 0.954166666666667, 0.904166666666667, 0.904166666666667, 0.975, 0.991666666666667, 0.979166666666667, 0.975, 0.970833333333333, 0.954166666666667, 0.9, 0.958333333333333, 0.929166666666667, 0.958333333333333), Stimulus.ER = c(0.0125, 0.00416666666666667, 0.0125, 0.025, 0.0375, 0.00416666666666667, 0.025, 0.0125, 0.0458333333333333, 0.0166666666666667, 0.0875, 0.05, 0.0916666666666667, 0.0166666666666667, 0.0416666666666667, 0.0208333333333333, 0.0125, 0.0208333333333333, 0.0375, 0.0458333333333333, 0.0958333333333333, 0.0958333333333333, 0.025, 0.00833333333333333, 0.0208333333333333, 0.025, 0.0291666666666667, 0.0458333333333333, 0.1, 0.0416666666666667, 0.0708333333333333, 0.0416666666666667), asin.Stimulus.ER = c(0.0315400751462974, 0.0105133583820991, 0.0315400751462974, 0.0630801502925948, 0.0714356562826538, 0.0105133583820991, 0.0630801502925948, 0.0259003510988588, 0.115646942203090, 0.0420534335283965, 0.209501077929205, 0.114880852490312, 0.190564470141143, 0.0420534335283965, 0.0994938597735527, 0.0525667919104956, 0.0315400751462974, 0.0525667919104956, 0.0889805013914536, 0.110007218155652, 0.219248346598526, 0.218622673632042, 0.0630801502925948, 0.0210267167641983, 0.0525667919104956, 0.0511750292312336, 0.0735935086746939, 0.110007218155652, 0.229761704980625, 0.105133583820991, 0.161807920353369, 0.0994938597735527), cell = structure(as.integer(c(1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16)), .Label = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16"), class = c("ordered", "factor"))), .Names = c("activeMapping", "mapping.code", "oom", "stimulusName", "responseFinger", "Subject.a", "Subject", "Subject.n", "Block", "Stimulus.ACC", "Stimulus.ER", "asin.Stimulus.ER", "cell"), row.names = c("13/1/Q/index/17/25/1", "13/1/Q/index/21/29/2", "13/4/J/middle/17/25/1", "13/4/J/middle/21/29/2", "13/2/X/ring/17/25/1", "13/2/X/ring/21/29/2", "13/3/B/little/17/25/1", "13/3/B/little/21/29/2", "14/2/B/index/18/26/1", "14/2/B/index/22/30/2", "14/3/X/middle/18/26/1", "14/3/X/middle/22/30/2", "14/1/J/ring/18/26/1", "14/1/J/ring/22/30/2", "14/4/Q/little/18/26/1", "14/4/Q/little/22/30/2", "15/3/J/index/19/27/1", "15/3/J/index/23/36/2", "15/2/Q/middle/19/27/1", "15/2/Q/middle/23/36/2", "15/4/B/ring/19/27/1", "15/4/B/ring/23/36/2", "15/1/X/little/19/27/1", "15/1/X/little/23/36/2", "16/4/X/index/20/28/1", "16/4/X/index/24/41/2", "16/1/B/middle/20/28/1", "16/1/B/middle/24/41/2", "16/3/Q/ring/20/28/1", "16/3/Q/ring/24/41/2", "16/2/J/little/20/28/1", "16/2/J/little/24/41/2"), class = "data.frame") ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html