[R] Using anova() with glmmPQL()

[Ben Bolker]

Toby Marthews <toby.marthews <at> ouce.ox.ac.uk> writes:

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

[snip]
>  lmeres=lme(fixed=stmass~mnth*roostsitu,
> random=~1|subject/mnth,na.action=na.exclude)
> anova(object=lmeres,test="Chisq")

  By the way, 'test="Chisq"' has no effect -- it's silently
ignored -- in an anova run on an lme() object the F test
is always used (as your output suggests).

>               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

  By the way, I'm a bit puzzled/suspicious about
your setup here. You seem to measure each subject
once in each month, so your grouping variable subject/month
is confounded with the residual error.  I would suggest
random=~1|subject (which gives a very nearly identical
ANOVA table but makes more sense).  In this case, because
you have only two months, you could do almost the same
analysis by taking differences or averages of subjects
across months: a paired t-test on subjects for the effect
of month, a one-way ANOVA on subject averages for the
effect of roost, and a one-way ANOVA on subject differences
for the month x roost interaction.

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

  [snip lme with family="binomial"]

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

  [snip]

>    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?
> 

  Hmm.  This feels harder than it ought to be, and I've already
spent longer on it than I should, so I'm going to give you some
thoughts and encourage you to send this on to R-sig-mixed-models.

  I think the *best* way to do this would probably be to do
the analogue of the paired/grouped tests I suggested above:
for the month effect, test whether more individuals go from
light to heavy or heavy to light, and so forth.

  Treating this as a binomial problem:

 you can use glmmPQL and use the Wald test (see wald.test in the
aod package, as below) to test groups of parameters.  In principle
you could do this for glmer, too, but beware the Hauck-Donner
effect (e.g., the 'Jan x other' parameter in the glmer fit
has estimate -13.6, standard error 6840 [!])

 you can use (g)lmer and anova() to do Likelihood Ratio Tests
on the effects: this is a little bit dicey for small sample
sizes.

  I also question your approach a little bit here: it looks
like you have chosen the cutoff weight as the minimum weight
for the 'tree/Nov' category (control?).  That may make biological
sense, but it means that there is essentially no scope for
changing categories across months in the 'tree' treatment,
which is bound to induce a a treatment x month interaction
if there are changes in the other treatments ...


==============
library(ggplot2)
boxplot(stmass~roostsitu*mnth,data=dataf)
boxplot(stmass~mnth*roostsitu,data=dataf)

library(ggplot2)
ggplot(dataf,aes(mnth,stmass))+
  geom_point()+
  geom_line(aes(group=subject))+
  facet_grid(.~roostsitu)+
  geom_hline(yintercept=78,colour="red")

## or 
ggplot(dataf,aes(mnth,stmass,colour=roostsitu))+
  geom_point()+
  geom_line(aes(group=subject))
## I prefer the first

lmeres=lme(fixed=stmass~mnth*roostsitu,
  random=~1|subject/mnth,na.action=na.exclude,data=dataf)

lmeres2=lme(fixed=stmass~mnth*roostsitu,
  random=~1|subject,na.action=na.exclude,data=dataf)
anova(lmeres2)

## or equivalently
summary(aov(stmass~mnth*roostsitu+Error(subject),data=dataf))

plot(lmeres)
library(lattice)
plot(lmeres,mnth~resid(.))
plot(lmeres,interaction(mnth,roostsitu)~resid(.))

dataf$stmassheavy <- as.numeric(stmass>78)

library(lme4)
glmerfit <- glmer(stmassheavy~mnth*roostsitu+(1|subject),
                 family=binomial,
                 na.action=na.exclude,data=dataf)
g2 <- update(glmerfit,.~.-mnth:roostsitu)
g3 <- update(g2,.~.-mnth)
g4 <- update(g2,.~.-roostsitu)

## Could use likelihood ratio tests -- but these are
anova(glmerfit,g2)
anova(g2,g3)
anova(g2,g4)

library(MASS)


PQLfit <- glmmPQL(fixed=stmassheavy~mnth*roostsitu,
                  random=~1|subject,na.action=na.exclude,data=dataf,
                  family=binomial)
library(aod)

wald.test(vcov(PQLfit),fixef(PQLfit),
          Terms=3:5)
wald.test(vcov(PQLfit),fixef(PQLfit),
          Terms=6:8)


df2 <- as.data.frame.table(with(dataf,
     table(stmassheavy,mnth,roostsitu)))
ggplot(subset(df2,stmassheavy==1),
       aes(mnth,Freq/10))+
  geom_line(aes(group=1))+
  geom_point()+
  facet_grid(.~roostsitu)