Question concerning anova()

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...).

The book suggests doing something like:
mod1 <- gls(noBefore~pHarv*year, data= dat) # model without random term

mod2 <- lme(noBefore~pHarv*year, data= dat, random=~1|plant, method="REML") # random intercept

mod3 <- lme(noBefore~pHarv*year, data= dat, random=~year|plant, method="REML") # random intercept

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).

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.

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