[R] Using anova() with glmmPQL()

Toby Marthews toby.marthews at ouce.ox.ac.uk
Mon Jan 17 16:38:48 CET 2011

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).


              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 -


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:


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 ):


               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.

Toby Marthews

More information about the R-help mailing list