[R] Using anova() with glmmPQL()

[Toby Marthews]

Dear R HELP,

ABOUT glmmPQL and the anova command. Here is an example of a repeated-measures ANOVA focussing on the way starling masses vary according to (i) roost situation and (ii) time (two time points only).

library(nlme);library(MASS)
stmass=c(78,88,87,88,83,82,81,80,80,89,78,78,85,81,78,81,81,82,76,74,79,73,79,75,77,78,80,78,83,84,77,68,75,70,74,84,80,75,76,75,85,88,86,95,100,87,98,86,89,94,84,88,91,96,86,87,93,87,94,96,91,90,87,84,86,88,92,96,83,85,90,87,85,81,84,86,82,80,90,77)
roostsitu=factor(c("tree","tree","tree","tree","tree","tree","tree","tree","tree","tree","nest-box","nest-box","nest-box","nest-box","nest-box","nest-box","nest-box","nest-box","nest-box","nest-box","inside","inside","inside","inside","inside","inside","inside","inside","inside","inside","other","other","other","other","other","other","other","other","other","other","tree","tree","tree","tree","tree","tree","tree","tree","tree","tree","nest-box","nest-box","nest-box","nest-box","nest-box","nest-box","nest-box","nest-box","nest-box","nest-box","inside","inside","inside","inside","inside","inside","inside","inside","inside","inside","other","other","other","other","other","other","other","other","other","other"),levels=c("tree","nest-box","inside","other"))
mnth=factor(c(rep("Nov",times=40),rep("Jan",times=40)),levels=c("Nov","Jan"))
subjectnum=c(1:10,1:10,1:10,1:10,1:10,1:10,1:10,1:10)
subject=factor(paste(roostsitu,subjectnum,sep=""))
dataf=data.frame(mnth,roostsitu,subjectnum,subject,stmass)
lmeres=lme(fixed=stmass~mnth*roostsitu,random=~1|subject/mnth,na.action=na.exclude)
anova(object=lmeres,test="Chisq")

              numDF denDF   F-value p-value
(Intercept)        1    36 31143.552  <.0001
mnth               1    36    95.458  <.0001
roostsitu          3    36    10.614  <.0001
mnth:roostsitu     3    36     0.657  0.5838

I can conclude from this that variation with both roost situation and month are significant, but with no interaction term. So far so good. However, say I were interested only in whether or not those starlings were heavier or lighter than 78g: seemingly, I could change my response variable and analyse like this -

stmassheavy=ifelse(stmass>78,1,0)
lmeres1=lme(fixed=stmassheavy~mnth*roostsitu,random=~1|subject/mnth,na.action=na.exclude,family=binomial)
anova(object=lmeres1,test="Chisq")

but I get errors doing that. After a certain amount of web searching, I find that I'm supposed to use glmmPQL for this so I tried:

lmeres2=glmmPQL(fixed=stmassheavy~mnth*roostsitu,random=~1|subject/mnth,na.action=na.exclude,family=binomial)
anova(object=lmeres2,test="Chisq")

The glmmPQL command runs, but I get "Error in anova.glmmPQL(object = lmeres, test = "Chisq") : 'anova' is not available for PQL fits". Looking into this, I find that I am not supposed to use the anova command in conjunction with glmmPQL (several posts from Brian Ripley http://r.789695.n4.nabble.com/R-glmmPQL-in-2-3-1-td808574.html and http://www.biostat.wustl.edu/archives/html/s-news/2002-06/msg00055.html and http://www.mail-archive.com/r-help@stat.math.ethz.ch/msg46894.html ) even though it appears that the earlier versions of glmmPQL did allow the anova command to work (before ~2004).
   However, I couldn't find any other way to run a repeated-measures ANOVA with famiy=binomial. After a while longer on Google, I found a 'workaround' from Spencer Graves (on http://markmail.org/message/jddj6aq66wdidrog#query:how%20to%20use%20anova%20with%20glmmPQL+page:1+mid:jddj6aq66wdidrog+state:results ):

class(lmeres2)="lme"
anova(object=lmeres2,test="Chisq")

               numDF denDF   F-value p-value
(Intercept)        1    36 182.84356  <.0001
mnth               1    36 164.57288  <.0001
roostsitu          3    36  17.79263  <.0001
mnth:roostsitu     3    36   3.26912  0.0322

Which does give me a result and tells me that the interaction term is significant here. HOWEVER, on that link Douglas Bates told Spencer Graves that this wasn't an approprate method.

I haven't found any other workarounds for this except some general advice that I should move onto using the lmer command (which I can't do because I need to get p-values for my fits and according to https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html I won't get those from lmer).

My questions are: (1) Is lmer the only way to do a binomial repeated-measures ANOVA in R? (which means that there's no way to do such an ANOVA in R without losing the p-values) and (2) if I am supposed to be using glmmPQL for this simple situation, what am I doing wrong?

Thanks very much for any help anyone can give me.

best,
Toby Marthews