[R-sig-ME] [R] Question concerning anova()

[Rainer M Krug]

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> 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:      Rainer at krugs.de
>>
>> Skype:      RMkrug
>>
>> ______________________________________________
>> R-help at 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:      Rainer at krugs.de

Skype:      RMkrug