[R] lme: reproducing example

[Pascal A. Niklaus]

Karl Knoblick wrote:

>Dear R-community!
>
>I still have the problem reproducing the following
>example using lme.
>
>id<-factor(rep(rep(1:5,rep(3,5)),3))
>factA <- factor(rep(c("a1","a2","a3"),rep(15,3)))
>factB <- factor(rep(c("B1","B2","B3"),15))
>Y<-numeric(length=45)
>Y[ 1: 9]<-c(56,52,48,57,54,46,55,51,51)
>Y[10:18]<-c(58,51,50,54,53,46,54,50,49)
>Y[19:27]<-c(53,49,48,56,48,52,52,52,50)
>Y[28:36]<-c(55,51,46,57,49,50,55,51,47)
>Y[37:45]<-c(56,48,51,58,50,48,58,46,52)
>df<-data.frame(id, factA, factB, Y) 
>df.aov <- aov(Y ~ factA*factB + Error(factA:id),
>data=df)
>summary(df.aov)
>
>Is there a way to get the same results with lme as 
>with aov with Error()? HOW???
>
>One idea was the following:
>df$factAid=factor(paste(as.character(df$factA),":",as.character(df$id),sep=""))
>df.lme <-
>lme(Y~factA*factB,df,random=~1|factAid,method="REML") 
>The degrees of freedom look right, but the F values
>don't match aov.
>
>Hope somebody can help! Thanks!!
>
>Karl
>  
>
Hmmm, strange, it works if I use factB:id as plot... it also works when 
I use factA:id as plot and replace your Y's by random numbers... is this 
a problem with convergence?

Pascal


 > df$Y=rnorm(45)
 > summary(aov(Y ~ factB*factA + Error(id:factA),data=df))

Error: id:factA
          Df  Sum Sq Mean Sq F value Pr(>F)
factA      2  2.9398  1.4699  0.9014 0.4318
Residuals 12 19.5675  1.6306

Error: Within
            Df  Sum Sq Mean Sq F value   Pr(>F)
factB        2  7.1431  3.5716  7.4964 0.002956 **
factB:factA  4  4.2411  1.0603  2.2254 0.096377 .
Residuals   24 11.4345  0.4764
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

 > anova(lme(Y ~ factB*factA ,data=df, random = ~ 1 | plot))
            numDF denDF  F-value p-value
(Intercept)     1    24 0.014294  0.9058
factB           2    24 7.496097  0.0030
factA           2    12 0.901489  0.4318
factB:factA     4    24 2.225317  0.0964

Pascal


 > summary(aov(Y ~ factA*factB + Error(factB:id)))

Error: factB:id
          Df Sum Sq Mean Sq F value    Pr(>F)
factB      2 370.71  185.36  51.488 1.293e-06 ***
Residuals 12  43.20    3.60
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1

Error: Within
            Df Sum Sq Mean Sq F value  Pr(>F)
factA        2  9.911   4.956  1.6248 0.21788
factA:factB  4 45.556  11.389  3.7341 0.01686 *
Residuals   24 73.200   3.050

 > df$plot <- factor(paste(df$factB,df$id))
 > anova(lme(Y ~ factB*factA , data=df, random = ~1 | plot))
            numDF denDF  F-value p-value
(Intercept)     1    24 33296.02  <.0001
factB           2    12    51.47  <.0001
factA           2    24     1.63  0.2178
factB:factA     4    24     3.73  0.0168
 >