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
boneheaded (?) question about SS in anova (lmer vs lme)
2 messages · Ben Bolker, Douglas Bates
3 days later
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>