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
>