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.k...@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.k...@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: 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.
--
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.
