[R] lme: reproducing example

[Karl Knoblick]

Thanks!
I think the minor differences taking the values with
rnorm result of the homogen distribution without an
effect. But the results of aov and lme should be
similiar for data with effects, too (at least for
simple and balanced designs).

Karl

 --- "Pascal A. Niklaus" <Pascal.Niklaus at unibas.ch>
schrieb: > 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
>  >
>