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