[R] Question concerning anova()
Rainer M Krug
r.m.krug at gmail.com
Wed Aug 22 18:04:55 CEST 2012
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> 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> 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
>
>
>
More information about the R-help
mailing list