[R] anova on binomial LMER objects

From: Robert Bagchi <bagchi.r_at_gmail.com>
Date: Fri 23 Sep 2005 - 05:51:45 EST


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@sheffield.ac.uk
   bagchi.r@gmail.com

http://www.shef.ac.uk/aps/apsrtp/bagchi-r

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Received on Fri Sep 23 06:00:42 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:40:25 EST