[R] greco-latin square
Spencer Graves
spencer.graves at pdf.com
Thu May 18 16:41:28 CEST 2006
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
More information about the R-help
mailing list