On 23/08/12 08:34, Alan Haynes wrote:
If you have a copy available, Zuur et al 2009 Mixed effects models and extensions in ecology with R is a good book and describes a procedure well. Almost the whole book is based on lme and also has examples of variance/correlation structures which might be useful to you (although you already seem to what what youre doing with them...).
Great - I'll look into it. I have luckily access via the University.
The book suggests doing something like: mod1 <- gls(noBefore~pHarv*year, data= dat) # model without random term
OK - makes sense.
mod2 <- lme(noBefore~pHarv*year, data= dat, random=~1|plant, method="REML") # random intercept
OK.
mod3 <- lme(noBefore~pHarv*year, data= dat, random=~year|plant, method="REML") # random intercept
Question concerning the "random intercept" you mention - I assume this should be "random effect of year" ?
anova(mod1, mod2, mod3) then if you accept mod2: mod2.1 <- update(mod2, method="ML") # ML for fixed effects mod2.2 <- update(mod2.1, .~. - pHarv:year) # create the nested model of mod2.1 anova(mod2.1, mod2.2 beyond that you should create two more nested models (each with a fixed effect removed) and compare them back to mod2.2 (assuming you dont need the interaction).
OK - makes sense.
Where exactly testing the correlation structures would come in im not sure though. Also, you need to be aware of "testing on the boundary." I forget exactly where it comes in though (testing for the random effect I think). Thats covered by Zuur et al too.
I'll check it out - thanks. Cheers, Rainer
HTH Alan -------------------------------------------------- Email: aghay...@gmail.com <mailto:aghay...@gmail.com> Mobile: +41794385586 Skype: aghaynes On 22 August 2012 18:04, Rainer M Krug <r.m.k...@gmail.com <mailto:r.m.k...@gmail.com>> wrote: Further discussed on r-sig-mixed-models Rainer On 22/08/12 17:04, Bert Gunter wrote: 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 <mailto: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 <mailto: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 <tel:%2B33%20-%20%280%299%2053%2010%2027%2044> Cell: +33 - (0)6 85 62 59 98 <tel:%2B33%20-%20%280%296%2085%2062%2059%2098> Fax : +33 - (0)9 58 10 27 44 <tel:%2B33%20-%20%280%299%2058%2010%2027%2044> Fax (D): +49 - (0)3 21 21 25 22 44 <tel:%2B49%20-%20%280%293%2021%2021%2025%2022%2044> email: rai...@krugs.de <mailto:rai...@krugs.de> Skype: RMkrug ________________________________________________ R-help@r-project.org <mailto:R-help@r-project.org> mailing list https://stat.ethz.ch/mailman/__listinfo/r-help <https://stat.ethz.ch/mailman/listinfo/r-help> PLEASE do read the posting guide http://www.R-project.org/__posting-guide.html <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 <tel:%2B33%20-%20%280%299%2053%2010%2027%2044> Cell: +33 - (0)6 85 62 59 98 <tel:%2B33%20-%20%280%296%2085%2062%2059%2098> Fax : +33 - (0)9 58 10 27 44 <tel:%2B33%20-%20%280%299%2058%2010%2027%2044> Fax (D): +49 - (0)3 21 21 25 22 44 <tel:%2B49%20-%20%280%293%2021%2021%2025%2022%2044> email: rai...@krugs.de <mailto:rai...@krugs.de> Skype: RMkrug _________________________________________________ R-sig-mixed-models@r-project.__org <mailto:r-sig-mixed-mod...@r-project.org> mailing list https://stat.ethz.ch/mailman/__listinfo/r-sig-mixed-models <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>
-- 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.
