[R] help with glmmPQL

[Douglas Bates]

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.

I don't think that would be appropriate.  AFAIK the logLik generic 
applied to a glmmPQL object does not return the likelihood of the fitted 
model or an approximation to that likelihood.

That's why there is a difference between the values of the likelihood 
for the same model fit by glmmPQL and fit using the PQL method in GLMM 
from the lme4 package.

I will add anova methods for GLMM objects whenever I manage to dig 
myself out of the current rewrite of the representation of linear 
mixed-effects models.

> 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.
>      I looked briefly at the code for "anova.lme", and it looks to me 
> like "glmmPQL" could be added to the list of allowed options in the 
> following test roughly half way through the code:
>       if (!all(match(termsClass, c("gls", "gnls", "lm", "lmList",
>            "lme", "nlme", "nlsList", "nls"), 0))) {
>            stop(paste("Objects must inherit from classes \"gls\", 
> \"gnls\"",
>                "\"lm\",\"lmList\", \"lme\",\"nlme\",\"nlsList\", or 
> \"nls\""))
>        }
> 
>      hope this helps.      spencer graves
> 
> A.J. Rossini wrote:
> 
>> For better or worse, it's holidays in the states.  Very amusing for me
>> being in a non-Thanksgiving celebrating country.
>>
>> In addition, it's not a problem.  The complaint is valid.  Probably no
>> one has coded up the right solution yet for comparison.
>>
>> I can't recall if one would want those statistics for a binomial
>> random effects model, but I do recall some issues with model
>> comparison in that setting, though they are a bit dated (say, 2 years
>> or so).
>>
>> On Fri, 26 Nov 2004 09:31:40 +0700, Andrew Criswell
>> <r-stats at arcriswell.com> wrote:
>>  
>>
>>> Hello:
>>>
>>> Will someone PLEASE help me with this problem. This is the third time
>>> I've posted it.
>>>
>>> When I appply anova() to two equations estimated using glmmPQL, I get a
>>> complaint,
>>>
>>>   
>>>
>>>> anova(fm1, fm2)
>>>>     
>>>
>>> Error in anova.lme(fm1, fm2) : Objects must inherit from classes "gls",
>>> "gnls" "lm","lmList", "lme","nlme","nlsList", or "nls"
>>>
>>>   
>>> The two equations I estimated are these:
>>>
>>>   
>>>
>>>> fm1 <- glmmPQL(choice ~ day + stereotypy,
>>>>     
>>>
>>> +                random = ~ 1 | bear, data = learning, family = 
>>> binomial)
>>>
>>>   
>>>
>>>> fm2 <- glmmPQL(choice ~ day + envir + stereotypy,
>>>>     
>>>
>>> +                random = ~ 1 | bear, data = learning, family = 
>>> binomial)
>>>
>>> Individually, I get results from anova():
>>>
>>>   
>>>
>>>> anova(fm1)
>>>>     
>>>
>>>          numDF denDF   F-value p-value
>>> (Intercept)     1  2032   7.95709  0.0048
>>> day             1  2032 213.98391  <.0001
>>> stereotypy      1  2032   0.42810  0.5130
>>>
>>>   
>>>
>>>> anova(fm2)
>>>>     
>>>
>>>          numDF denDF   F-value p-value
>>> (Intercept)     1  2031   5.70343  0.0170
>>> day             1  2031 213.21673  <.0001
>>> envir           1  2031  12.50388  0.0004
>>> stereotypy      1  2031   0.27256  0.6017
>>>
>>>   
>>> I did look through the archives but didn't finding anything relevant to
>>> my problem.
>>>
>>> Hope someone can help.
>>>
>>> ANDREW
>>> ____________________________
>>>       _
>>> platform i586-mandrake-linux-gnu
>>> arch     i586
>>> os       linux-gnu
>>> system   i586, linux-gnu
>>> status
>>> major    2
>>> minor    0.0
>>> year     2004
>>> month    10
>>> day      04
>>> language R
>>>
>>> -- 
>>> Andrew R. Criswell, Ph.D.
>>> Graduate School, Bangkok University
>>>
>>> ______________________________________________
>>> 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
>>>
>>>   
>>
>>
>>
>>  
>>
>