[R-sig-ME] boneheaded (?) question about SS in anova (lmer vs lme)

[Douglas Bates]

On 8/27/07, Ben Bolker <bolker at zoo.ufl.edu> wrote:

>   I've been working through some basic nested analyses of variance
> just to try to understand the details, and where lme(r)-style analyses
> depart from classical F-ratio tests.

Thanks for the questions, Ben, and for the well-constructed example.

Before going into detail about how the analysis of variance table is
constructed for an lmer model, I thought that I would do a bit of
plotting of the data.  Some of the plots indicate that a linear mixed
model is probably not appropriate.  Consider

library(lattice)
print(dotplot(reorder(PATCH, ALGAE) ~ ALGAE,  urchins, ylab = "Patch"))

and the same plot for the TREAT variable (PDF files attached).

I would be much more inclined to use a generalized linear mixed model
with a Poisson conditional distribution than a normal conditional
distribution for these data.

>   I'm having trouble understanding the results of anova(lmer(...)) ,
> which don't match up in any way I can figure out with direct
> calculations or with the results of lme().
>
> Andrew and Underwood did a simple nested experiment on urchins
> and algal cover -- 4 treatments, 4 patches per treatment, 5 samples
> per patch.
>
> ## get data
> datafile <- "http://www.zoology.unimelb.edu.au/qkstats/chpt9/andrew.csv"
> urchins <- read.csv(file=datafile,
>                     colClasses=c(rep("factor",3),"numeric"))
>
> ## calculate SS/F statistic by hand
> attach(urchins)
> tmeans <- tapply(ALGAE,TREAT,mean)[c(1,4,3,2)]
> pmeans <- tapply(ALGAE,PATCH,mean)[c(1,9:16,2:8)]
> ss.treat <- 20*sum((tmeans-mean(ALGAE))^2)
> ss.patch <- 5*sum((rep(tmeans,each=4)-pmeans)^2)
> fratio = (ss.treat/3)/(ss.patch/12)
> c(ss.treat,ss.patch,fratio)
> ## results, rounded: 14429.14 21241.95     2.72
> detach(urchins)
>
> ## same model with lme
> library(nlme)
> lme1 = lme(ALGAE~TREAT,random=~1|PATCH,data=urchins,method="REML")
> anova(lme1)
> detach("package:nlme")
>
> ##             numDF denDF   F-value p-value
> ## (Intercept)     1    64 18.555081  0.0001
> ## TREAT           3    12  2.717102  0.0913
>
> F value agrees with hand calculation above
> (for what it's worth, line 2 of ?anova.lme
> states that result of applying anova.lme to
> a single lme object will include the sums of squares,
> which seems to be false)
>
>
> ## 3. lmer: all parameter estimates agree with
> ##    lme fit, but ANOVA table is very odd --
> ##    don't know where these SS numbers come from??
>
> library(lme4)
> anova(lmer(ALGAE~TREAT+(1|PATCH),data=urchins,method="REML"))
> detach("package:lme4")
>
> ## Analysis of Variance Table
> ##       Df  Sum Sq Mean Sq
> ## TREAT  3 2434.27  811.42
>
> why is treatment SS not the same as ss.treat?
> (ss.treat also matches the results of aov())
>
>   Hope I haven't done something boneheaded, nor included too much
> nor too little information.
>
>   thanks
>     Ben Bolker
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Urchins1.pdf
Type: application/pdf
Size: 8948 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20070831/aef3e375/attachment.pdf>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Urchins2.pdf
Type: application/pdf
Size: 7926 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20070831/aef3e375/attachment-0001.pdf>