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

Ben Bolker bolker at zoo.ufl.edu
Tue Aug 28 03:49:07 CEST 2007


  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.

  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




More information about the R-sig-mixed-models mailing list