Question concerning anova()

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.

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.

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