[R] anova on binomial LMER objects
Robert Bagchi
bagchi.r at gmail.com
Thu Sep 22 21:51:45 CEST 2005
Dear R users,
I have been having problems getting believable estimates from anova on a
model fit from lmer. I get the impression that F is being greatly
underestimated, as can be seen by running the example I have given below.
First an explanation of what I'm trying to do. I am trying to fit a glmm
with binomial errors to some data. The experiment involves 10
shadehouses, divided between 2 light treatments (high, low). Within each
shadehouse there are 12 seedlings of each of 2 species (hn & sl). 3
damage treatments (0, 0.1, 0.25 leaf area removal) were applied to the
seedlings (at random) so that there are 4 seedlings of each
species*damage treatment in each shadehouse. There maybe a shadehouse
effect, so I need to include it as a random effect. Light is applied to
a shadehouse, so it is outer to shadehouse. The other 2 factors are
inner to shadehouse.
We want to assess if light, damage and species affect survival of
seedlings. To test this I fitted a binomial mixed effects model with
lmer (actually with quasibinomial errors). THe summary function suggests
a large effect of both light and species (which agrees with graphical
analysis). However, anova produces F values close to 0 and p values
close to 1 (see example below).
Is this a bug, or am I doing something fundamentally wrong? If anova
doesn't work with lmer is there a way to perform hypothesis tests on
fixed effects in an lmer model? I was going to just delete terms and
then do liklihood ratio tests, but according to Pinheiro & Bates (p. 87)
that's very untrustworthy. Any suggestions?
I'm using R 2.1.1 on windows XP and lme4 0.98-1
Any help will be much appreciated.
many thanks
Robert
###############################
The data are somewhat like this
#setting up the dataframe
bm.surv<-data.frame(
house=rep(1:10, each=6),
light=rep(c("h", "l"), each=6, 5),
species=rep(c("sl", "hn"), each=3, 10),
damage=rep(c(0,.1,.25), 20)
)
bm.surv$survival<-ifelse(bm.surv$light=="h", rbinom(60, 4, .9),
rbinom(60, 4, .6)) # difference in probablility should ensure a
light effect
bm.surv$death<-4-bm.surv$survival
# fitting the model
m1<-lmer(cbind(survival, death)~light+species+damage+(1|house),
data=bm.surv, family="quasibinomial")
summary(m1) # suggests that light is very significant
Generalized linear mixed model fit using PQL
Formula: cbind(survival, death) ~ light + species + damage + (1 | table)
Data: bm.surv
Family: quasibinomial(logit link)
AIC BIC logLik deviance
227.0558 239.6218 -107.5279 215.0558
Random effects:
Groups Name Variance Std.Dev.
table (Intercept) 1.8158e-09 4.2613e-05
Residual 3.6317e+00 1.9057e+00
# of obs: 60, groups: table, 10
Fixed effects:
Estimate Std. Error DF t value Pr(>|t|)
(Intercept) 2.35140 0.36832 56 6.3841 3.581e-08 ***
lightl -1.71517 0.33281 56 -5.1535 3.447e-06 ***
speciessl -0.57418 0.30085 56 -1.9085 0.06145 .
damage 1.49963 1.46596 56 1.0230 0.31072
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) lightl spcssl
lightl -0.665
speciessl -0.494 0.070
damage -0.407 -0.038 -0.017
anova(m1) # very low F value for light, corresponding to p
values approaching 1
Analysis of Variance Table
Df Sum Sq Mean Sq Denom F value Pr(>F)
light 1 0.014 0.014 56.000 0.0018 0.9661
species 1 0.002 0.002 56.000 0.0002 0.9887
damage 1 0.011 0.011 56.000 0.0014 0.9704
--
Robert Bagchi
Animal & Plant Science
Alfred Denny Building
University of Sheffield
Western Bank
Sheffield S10 2TN
UK
t: +44 (0)114 2220062
e: r.bagchi at sheffield.ac.uk
bagchi.r at gmail.com
http://www.shef.ac.uk/aps/apsrtp/bagchi-r
More information about the R-help
mailing list