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")
Your self-contained example interested me, because both lme{nlme} and
lmer{lme4} returned only error messages for my saturated models, but
'aov' solved it fine. I'd be interested in hearing from someone
else if
it's actually possible to solve your problem using either 'lme' or
'lmer'; typical failed attempts appear below. Both stopped with a
singular error message, and neither has a 'singular.ok' option.
If I've read your post correctly, you are asking, "why aren't
the
contrasts that R is using to estimate the SS in the ANOVA table also not
orthogonal?" Are you confused, because when you look at summary.lm(m1)
when 'm1' is fit with and without one of the interactions, you see
different numbers? Part of the magic of least squares is that it works
by thinking of N observations as a point in N-dimensional space and then
projecting the vector of observations on different sub-spaces. The
particular parameterization, which is what you see with summary.lm(m1),
does not matter as long as you have specified the correct subspaces.
You can find more on this in almost any book that discusses regression
and least squares from this perspective.
hope this help.
Spencer Graves
p.s. Your self-contained example helped me a lot. I don't know if I
answered your question, if I did, I likely would not have without that
example. I almost never use aov, because I rarely have such a simple,
well-balanced experiment. Failed attempts to solve the same problem
with lme{nlme} and lmer{lme4} appear below. If your example had NOT
been so balanced, it would not have been feasible to get a sensible
answer without something like lme or lmer.
> library(nlme)
> m1.lme <-lme(asin.Stimulus.ER ~
+ responseFinger + oom +
+ (stimulusName + mapping.code)^2,
+ random=~stimulusName:mapping.code|Subject,
+ data=w.tmp)
Error in MEEM(object, conLin, control$niterEM) :
Singularity in backsolve at level 0, block 1
In addition: Warning message:
Fewer observations than random effects in all level 1 groups in:
lme.formula(asin.Stimulus.ER ~ responseFinger + oom + (stimulusName +
library(lme4)
> m1 <- lmer(asin.Stimulus.ER ~
+ responseFinger + oom +
+ (stimulusName + mapping.code)^2
+ +(mapping.code:stimulusName|Subject),
+ data=w.tmp)
Error in lmer(asin.Stimulus.ER ~ responseFinger + oom + (stimulusName +
:Leading minor of order 17 in downdated X'X is not positive definite
sessionInfo()
Version 2.3.0 (2006-04-24)
i386-pc-mingw32
attached base packages:
[1] "methods" "stats" "graphics"
"grDevices" "utils" "datasets"
[7] "base"
other attached packages:
lme4 lattice Matrix
"0.995-2" "0.13-8" "0.995-10"
Steven Lacey wrote:> 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 at 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
## Thank you for a very nice example. I would like to make two
## expansions on the discussion that Spencer Graves gave.
## Rich
## 1. Multiple regression is a sequential algorithm. In the
## sequential anova table the sums of squares for non-orthogonal
## predictor variables depends on the order in which they are brought
## into the model. As part of the regression algorithm each of the
## predictor variables, in this case the dummy variables corresponding
## to the contrasts, is orthogonalized against all previous dummy
## variables.
## We can see that the algorithm has done its task by looking at the
## projection of Y=asin.Stimulus.ER onto the orthogonalized subspaces
## associated with each of the multi-degree-of-freedom terms in the
## model. We do so by projecting Y onto each of those subspaces with
## the proj() function.
m1 <- aov(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)
m1.proj <- proj(m1)
tmp <- crossprod(m1.proj)
## zapsmall(tmp)
dimnames(tmp) <- rep(list(abbreviate(dimnames(tmp)[[1]])), 2)
## The projections onto the subspaces are orthogonal to each other.
zapsmall(tmp, digits=3)
## > zapsmall(tmp, digits=3)
## (In) stmN rspF oom mpp. m.:S stN:. sN:.:
## (In) 0.2561 0.000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## stmN 0.0000 0.005 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## rspF 0.0000 0.000 0.0337 0.0000 0.0000 0.0000 0.0000 0.0000
## oom 0.0000 0.000 0.0000 0.0061 0.0000 0.0000 0.0000 0.0000
## mpp. 0.0000 0.000 0.0000 0.0000 0.0282 0.0000 0.0000 0.0000
## m.:S 0.0000 0.000 0.0000 0.0000 0.0000 0.0204 0.0000 0.0000
## stN:. 0.0000 0.000 0.0000 0.0000 0.0000 0.0000 0.0223 0.0000
## sN:.: 0.0000 0.000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0139
## >
## The sums of squares of the projections are exactly the sums of
## squares in the ANOVA table.
as.matrix(apply(m1.proj, 2, function(x) sum(x^2)))
## > as.matrix(apply(m1.proj, 2, function(x) sum(x^2)))
## [,1]
## (Intercept) 0.256144760
## stimulusName 0.005001724
## responseFinger 0.033730254
## oom 0.006146134
## mapping.code 0.028189996
## mapping.code:Subject.n 0.020374325
## stimulusName:mapping.code 0.022317815
## stimulusName:mapping.code:Subject.n 0.013902509
## >
## 2. We can ask the summary function to do all the right tests.
## There is no need to do them by hand. The Error() function works
## directly with the orthogonal subspaces and recognizes which terms
## are the numerators and denominators for the F tests.
m1e <- aov(asin.Stimulus.ER ~
(stimulusName + mapping.code + responseFinger + oom)^2 +
Error(mapping.code:Subject.n),
data=w.tmp)
summary(m1e)
## > m1e <- aov(asin.Stimulus.ER ~
## + (stimulusName + mapping.code + responseFinger + oom)^2 +
## + Error(mapping.code:Subject.n),
## + data=w.tmp)
## Warning message:
## Error() model is singular in: aov(asin.Stimulus.ER ~ (stimulusName +
mapping.code +
responseFinger +
## >
## > summary(m1e)
##
## Error: mapping.code:Subject.n
## Df Sum Sq Mean Sq F value Pr(>F)
## mapping.code 3 0.0281900 0.0093967 1.8448 0.2794
## Residuals 4 0.0203743 0.0050936
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## stimulusName 3 0.005002 0.001667 1.4391 0.280137
## responseFinger 3 0.033730 0.011243 9.7048 0.001569 **
## oom 3 0.006146 0.002049 1.7684 0.206609
## stimulusName:mapping.code 3 0.022318 0.007439 6.4212 0.007677 **
## Residuals 12 0.013903 0.001159
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05
'.' 0.1 ' ' 1
## >
## The warning about the singular model is an indicator that more
## dummy variables have been generated than are needed. This is
## consistent with your observation that any of the twoway
## interactions could be used to capture those interactions.