[R] model comparison with mixed effects glm

[Spencer Graves]

	  It's not clear from your email what you tried, but "anova" to compare 
two glmmPQL fits would not work for me.  I switched to lmer and got 
reasonable answers.  The first includes what worked for me then what I 
tried unsuccessfully with glmmPQL:

library(MASS) # for the "bacteria" data used in this example
library(lme4)
Fit.ID0. <- lmer(y~trt+(1|ID),
                family=binomial, data=bacteria,
                  method="Laplace")
Fit.ID1. <- lmer(y~trt+I(week>2)+(1|ID),
                family=binomial, data=bacteria,
                  method="Laplace")
 > anova(Fit.ID1., Fit.ID0.)
Data: bacteria
Models:
Fit.ID0.: y ~ trt + (1 | ID)
Fit.ID1.: y ~ trt + I(week > 2) + (1 | ID)
          Df      AIC      BIC   logLik  Chisq Chi Df Pr(>Chisq)
Fit.ID0.  4  214.324  227.898 -103.162
Fit.ID1.  5  202.261  219.230  -96.131 14.062      1  0.0001768 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

### However, I could not get the same command to work with a "glm" object:
fit0 <- glm(y ~ trt + I(week > 2),
          family = binomial, data = bacteria)
 > anova(Fit.ID1., fit0)
Error in FUN(X[[1]], ...) : no slot of name "call" for this object of 
class "glm"
In addition: Warning message:
extra arguments discarded in: logLik.glm(X[[2]], ...)

### To get around that, I computed 2*log(likelihood ratio) manually:

lglk0 <- logLik(fit0)
lglk.ID1. <- logLik(Fit.ID1.)
pchisq(as.numeric(chisq.ID.), 1, lower=FALSE)
 > [1] 0.008545848

###### AFTER QUITTING AND RESTARTING R
#*****anova(glmmPQL)
library(MASS)
library(nlme) # will be loaded automatically if omitted
fit.ID0 <- glmmPQL(y ~ trt, random = ~ 1 | ID,
             family = binomial, data = bacteria)
fit.ID1 <- glmmPQL(y ~ trt + I(week > 2), random = ~ 1 | ID,
             family = binomial, data = bacteria)
anova(fit.ID1)
             numDF denDF  F-value p-value
(Intercept)     1   169 35.01574  <.0001
trt             2    47  1.58393  0.2159
I(week > 2)     1   169 20.11802  <.0001
 > anova(fit.ID1, fit.ID0)
Error in anova.lme(fit.ID1, fit.ID0) : Objects must inherit from classes 
"gls", "gnls" "lm","lmList", "lme","nlme","nlsList", or "nls"
 > sessionInfo()
R version 2.2.1, 2005-12-20, i386-pc-mingw32

attached base packages:
[1] "methods"   "stats"     "graphics"  "grDevices" "utils"     "datasets"
[7] "base"

other attached packages:
       nlme       MASS
"3.1-68.1"   "7.2-24"
#####################

	  Hope this helps.
	  spencer graves
p.s.  One of the most active contributors to R recently commented that 
he rarely replies to emails that do NOT follow the posting guide 
"www.R-project.org/posting-guide.html".  This is not just a silly 
bureaucratic requirement but a process for helping people with questions 
work through their problem and either find an answer themselves or help 
them express their question in terms that potential repsondents can more 
easily be understood.  I know neither your modelA nor modelB, so I'm 
forced to guess.  You almost certainly would have gotten a more 
informative answer quicker if your question had been more consistent 
with the posting guide.

Paula den Hartog wrote:

> I use model comparison with glms without mixed effects with
> anova(modelA,modelB), 
> with mixed effects glm (glmmPQL), this doesn't work. Is there a way to
> compare model fits with glmmPQL's?
>  
> Paula M. den Hartog
> Behavioural Biology
> Institute of Biology Leiden
> Leiden University
>  
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at 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