[R] Question concerning anova()

[Bert Gunter]

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



-- 

Bert Gunter
Genentech Nonclinical Biostatistics

Internal Contact Info:
Phone: 467-7374
Website:
http://pharmadevelopment.roche.com/index/pdb/pdb-functional-groups/pdb-biostatistics/pdb-ncb-home.htm