[R] help with glmmPQL
Deepayan Sarkar
deepayan at stat.wisc.edu
Fri Nov 26 20:09:30 CET 2004
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
More information about the R-help
mailing list