[R] lme vs. aov with Error term again

array chip arrayprofile at yahoo.com
Wed Oct 1 18:41:26 CEST 2003


Hi all,

Sent the following question yesterday, but haven't got
any suggestions yet. So just trying again, can anyone
comment on the problem that I have? Thank you!

-------------

Hi,

I have a question about using "lme" and "aov" for the
following dataset. If I understand correctly, using
"aov" with an Error term in the formula is equivalent
to using "lme" with default settings, i.e. both assume
compound symmetry correlation structure. And I have
found that equivalency in the past. However, with the
follwing dataset, I got different answers with using
"lme" and using "aov", can anyone explain what
happened here? I have 2 differnt response variables
"x" and "y" in the following dataset, they are
actually slightly different (only 3 values of them are
different). With "y", I achieved the equivalency
between "lme" and "aov"; but with "x", I got different
p values for the ANOVA table.

-------

x<-c(-0.0649,-0.0923,-0.0623,0.1809,0.0719,0.1017,0.0144,-0.1727,-0.1332,0.0986,0.304,-0.4093,0.2054,0.251,-0.1062,0.3833,0.0649,0.2908,0.1073,0.0919,0.1167,0.2369,0.306,0.1379)

y<-c(-0.0649,-0.0923,0.32,0.08,0.0719,0.1017,0.05,-0.1727,-0.1332,0.15,0.304,-0.4093,0.2054,0.251,-0.1062,0.3833,0.0649,0.2908,0.1073,0.0919,0.1167,0.2369,0.306,0.1379)

treat<-as.factor(c(1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2))
time<-as.factor(c(1,1,1,1,2,2,2,2,3,3,3,3,1,1,1,1,2,2,2,2,3,3,3,3))
sex<-as.factor(c('F','F','M','M','F','F','M','M','F','F','M','M','F','F','M','M','F','F','M','M','F','F','M','M'))
subject<-as.factor(c(rep(1:4,3),rep(5:8,3)))
xx<-cbind(x=data.frame(x),y=y,treat=treat,time=time,sex=sex,subject=subject)

######## using x as dependable variable

xx.lme<-lme(x~treat*sex*time,random=~1|subject,xx)
xx.aov<-aov(x~treat*sex*time+Error(subject),xx)

summary(xx.aov)

Error: subject
          Df   Sum Sq  Mean Sq F value  Pr(>F)  
treat      1 0.210769 0.210769  6.8933 0.05846 .
sex        1 0.005775 0.005775  0.1889 0.68627  
treat:sex  1 0.000587 0.000587  0.0192 0.89649  
Residuals  4 0.122304 0.030576                  
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.'
0.1 ` ' 1 

Error: Within
               Df  Sum Sq Mean Sq F value Pr(>F)
time            2 0.00102 0.00051  0.0109 0.9891
treat:time      2 0.00998 0.00499  0.1066 0.9002
sex:time        2 0.02525 0.01263  0.2696 0.7704
treat:sex:time  2 0.03239 0.01619  0.3458 0.7178
Residuals       8 0.37469 0.04684 

anova(xx.lme)
               numDF denDF  F-value p-value
(Intercept)        1     8 3.719117  0.0899
treat              1     4 5.089022  0.0871
sex                1     4 0.139445  0.7278
time               2     8 0.012365  0.9877
treat:sex          1     4 0.014175  0.9110
treat:time         2     8 0.120538  0.8880
sex:time           2     8 0.304878  0.7454
treat:sex:time     2     8 0.391012  0.6886

#### using y as dependable variable

xx.lme2<-lme(y~treat*sex*time,random=~1|subject,xx)
xx.aov2<-aov(y~treat*sex*time+Error(subject),xx)

summary(xx.aov2)

Error: subject
          Df   Sum Sq  Mean Sq F value Pr(>F)
treat      1 0.147376 0.147376  2.0665 0.2239
sex        1 0.000474 0.000474  0.0067 0.9389
treat:sex  1 0.006154 0.006154  0.0863 0.7836
Residuals  4 0.285268 0.071317               

Error: Within
               Df   Sum Sq  Mean Sq F value Pr(>F)
time            2 0.009140 0.004570  0.1579 0.8565
treat:time      2 0.012598 0.006299  0.2177 0.8090
sex:time        2 0.043132 0.021566  0.7453 0.5049
treat:sex:time  2 0.069733 0.034866  1.2050 0.3488
Residuals       8 0.231480 0.028935               

anova(xx.lme2)
               numDF denDF   F-value p-value
(Intercept)        1     8 3.0667809  0.1180
treat              1     4 2.0664919  0.2239
sex                1     4 0.0066516  0.9389
time               2     8 0.1579473  0.8565
treat:sex          1     4 0.0862850  0.7836
treat:time         2     8 0.2177028  0.8090
sex:time           2     8 0.7453185  0.5049
treat:sex:time     2     8 1.2049883  0.3488




More information about the R-help mailing list