[R] [R-sig-ME] Question concerning anova()
Rainer M Krug
r.m.krug at gmail.com
Thu Aug 23 14:19:19 CEST 2012
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: aghaynes at gmail.com <mailto:aghaynes at gmail.com>
> Mobile: +41794385586
> Skype: aghaynes
>
>
> On 22 August 2012 18:04, Rainer M Krug <r.m.krug at gmail.com <mailto:r.m.krug at 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.krug at gmail.com
> <mailto:r.m.krug at 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.krug at gmail.com
> <mailto:r.m.krug at 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: Rainer at krugs.de <mailto:Rainer at krugs.de>
>
> Skype: RMkrug
>
> ________________________________________________
> R-help at r-project.org <mailto:R-help at 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: Rainer at krugs.de <mailto:Rainer at krugs.de>
>
> Skype: RMkrug
>
>
>
>
>
> _________________________________________________
> R-sig-mixed-models at r-project.__org <mailto:R-sig-mixed-models at 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: Rainer at krugs.de
Skype: RMkrug
More information about the R-help
mailing list