[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