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 at 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
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.5838By 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 ANOVAwith 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 issignificant 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)