[R] help with glmmPQL

[Deepayan Sarkar]

On Friday 26 November 2004 12:35, Spencer Graves wrote:
> HI, DOUG & JOSE:
>
>       Is there some reason that "anova.lme" should NOT accept an
> object of class "glmmPQL" in the example below?  If you don't see one
> either, then I suggest you consider modifying the code as described
> below.
>
> HI, ANDREW:
>
>       I couldn't find your data "learning" in my Windows installation
> of R 2.0.0pat, which meant that I had to take the time to find
> another example before I could get the error message you described. 
> I got it from modifying the example in the documentation for
> "glmmPQL" as follows:
>
> fit1 <- glmmPQL(y ~ trt, random = ~ 1 | ID,
>                      family = binomial, data = bacteria)
> fit2 <- glmmPQL(y ~ trt + I(week > 2), random = ~ 1 | ID,
>                      family = binomial, data = bacteria)
> anova(fit1, fit2)
>     Error in anova.lme(fit1, fit2) : Objects must inherit from
> classes "gls", "gnls" "lm","lmList", "lme","nlme","nlsList", or "nls"
>
>       I then checked the class of fit1 and fit2:
>  > class(fit1)
>
> [1] "glmmPQL" "lme"
>
>  > class(fit2)
>
> [1] "glmmPQL" "lme"
>
>       As an experiment, I changed the class of fit1 and fit2:
>  > class(fit1) <- "lme"
>  > class(fit2) <- "lme"
>  > anova(fit1, fit2)
>
>      Model df      AIC      BIC    logLik   Test  L.Ratio p-value
> fit1     1  5 1054.623 1071.592 -522.3117
> fit2     2  6 1113.622 1133.984 -550.8111 1 vs 2 56.99884  <.0001
>
>       Unless someone like Doug or Jose tells us, "Don't do that", I
> would use these answers.

These likelihoods are (AFAIK) NOT the likelihoods of the models fitted, 
they are likelihoods of an lme model that approximates it. Thus, the 
test may not be appropriate. Having 'anova(fit1, fit2)' silently 
producing an answer would certainly be misleading.

E.g., lme4 produces the following:

> library(lme4)
> data(bacteria, package = "MASS")
> fit1 <- GLMM(y ~ trt, random = ~ 1 | ID, 
+              family = binomial, data = bacteria) 
> fit2 <- GLMM(y ~ trt + I(week > 2), random = ~ 1 | ID, 
+              family = binomial, data = bacteria) 
> logLik(fit1)
`log Lik.' -104.3780 (df=5)
> logLik(fit2)
`log Lik.' -97.68162 (df=6)

Deepayan