[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