[R] greco-latin square
Steven Lacey
slacey at umich.edu
Thu May 11 17:15:06 CEST 2006
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")
More information about the R-help
mailing list