Oops -- missed that. OTOH, my reply demonstrates the value of the mixed models list recommendation.
-- Bert On Wed, Aug 22, 2012 at 7:55 AM, Rainer M Krug <r.m.k...@gmail.com> wrote: > On 22/08/12 16:36, Bert Gunter wrote: >> >> Models with different fixed effects estimated by REML cannot be >> compared by anova. > > > I have seen that much in "Modern Applied Statistics in S", and therefore > have chosen the model = "ML" > >> >> In future, please post questions on mixed effects models on the >> r-sig-mixed-effects mailing lists. You're likely to receive more >> informative replies there, too. > > > Thanks - wasn't aware of this sig - I'll send the reply there as well. > > Thanks, > > Rainer > >> >> -- Bert >> >> On Wed, Aug 22, 2012 at 7:23 AM, Rainer M Krug <r.m.k...@gmail.com> wrote: >>> >>> Hi >>> >>> I am comparing four different linear mixed effect models, derived from >>> updating the original one. To compare these, I want to use anova(). I >>> therefore do the following (not reproducible - just to illustration >>> purpose!): >>> >>> dat <- loadSPECIES(SPECIES) >>> subs <- expression(dead==FALSE & recTreat==FALSE) >>> feff <- noBefore~pHarv*year # fixed effect in the model >>> reff <- ~year|plant # random effect in the model, where year >>> is >>> the >>> corr <- corAR1(form=~year|plant) # describing the within-group >>> correlation >>> structure >>> # >>> dat.lme <- lme( >>> fixed = feff, # fixed effect in >>> the >>> model >>> data = dat, >>> subset = eval(subs), >>> method = "ML", >>> random = reff, # random effect in >>> the >>> model >>> correlation = corr, >>> na.action = na.omit >>> ) >>> dat.lme.r1 <- update(dat.lme, random=~1|plant) >>> dat.lme.f1 <- update(dat.lme, fixed=noBefore~year) >>> dat.lme.r1.f1 <- update(dat.lme.r1, fixed=noBefore~year) >>> >>> >>> The anova is as follow: >>> >>>> anova(dat.lme, dat.lme.r1, dat.lme.f1, dat.lme.r1.f1) >>> >>> Model df AIC BIC logLik Test L.Ratio >>> p-value >>> dat.lme 1 9 1703.218 1733.719 -842.6089 >>> dat.lme.r1 2 7 1699.218 1722.941 -842.6089 1 vs 2 1.019230e-07 >>> 1 >>> dat.lme.f1 3 7 1705.556 1729.279 -845.7779 >>> dat.lme.r1.f1 4 5 1701.556 1718.501 -845.7779 3 vs 4 8.498318e-08 >>> 1 >>> >>> I have two questions: >>> 1) I am wondering why the "2 vs 3" does not give the Test values? >>> Is this because the two models are considered as "identical", which would >>> be >>> strange, due to the different logLik values. >>> >>> 2) If I want to compare all models among each other - is there a "best" >>> way? >>> I would be reluctant to do several ANOVA's, due to necessary corrections >>> for >>> multple tests (although this should not be a problem here?) >>> >>> I can obviously select the best model based on the AIC. >>> >>> Thanks in advance, >>> >>> Rainer >>> >>> -- >>> Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation >>> Biology, >>> UCT), Dipl. Phys. (Germany) >>> >>> Centre of Excellence for Invasion Biology >>> Stellenbosch University >>> South Africa >>> >>> Tel : +33 - (0)9 53 10 27 44 >>> Cell: +33 - (0)6 85 62 59 98 >>> Fax : +33 - (0)9 58 10 27 44 >>> >>> Fax (D): +49 - (0)3 21 21 25 22 44 >>> >>> email: rai...@krugs.de >>> >>> Skype: RMkrug >>> >>> ______________________________________________ >>> R-help@r-project.org mailing list >>> https://stat.ethz.ch/mailman/listinfo/r-help >>> PLEASE do read the posting guide >>> http://www.R-project.org/posting-guide.html >>> and provide commented, minimal, self-contained, reproducible code. >> >> >> >> > > > -- > Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology, > UCT), Dipl. Phys. (Germany) > > Centre of Excellence for Invasion Biology > Stellenbosch University > South Africa > > Tel : +33 - (0)9 53 10 27 44 > Cell: +33 - (0)6 85 62 59 98 > Fax : +33 - (0)9 58 10 27 44 > > Fax (D): +49 - (0)3 21 21 25 22 44 > > email: rai...@krugs.de > > Skype: RMkrug -- Bert Gunter Genentech Nonclinical Biostatistics Internal Contact Info: Phone: 467-7374 Website: http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
