[R-sig-ME] [R] Question concerning anova()
Rainer M Krug
r.m.krug at gmail.com
Wed Aug 22 17:56:11 CEST 2012
On 22/08/12 17:36, Alan Haynes wrote:
> Hi Rainer,
>
> 1) I *think* you dont get tests for your 2 vs 3 because your models arent nested which is a
> condition for using anova() in this way I think. I would suggest you write out your models rather
> than use update as you have been. The you'll know exactly whats nested and whats not.
OK - tried it out, and it makes sense.
>
> 2) Im not sure about a best way (you'll probably get different answers depending on who you ask). In
> any case, I believe its recommended to sort out your random effects before you start dealing with
> fixed effects. Then reduce your fixed effects if needs be.
That is what I was thinking about, effectively going along with the steps described in the Modern
Applied Statistics book. But further suggestions are welcome.
Rainer
>
> HTH
>
> Alan
>
>
>
> --------------------------------------------------
> Email: aghaynes at gmail.com <mailto:aghaynes at gmail.com>
> Mobile: +41794385586
> Skype: aghaynes
>
>
> On 22 August 2012 16:55, 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-sig-mixed-models
mailing list