[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