Skip to content

help with glmmPQL

8 messages · Andrew Criswell, A.J. Rossini, Spencer Graves +3 more

Original
Andrew Criswell
# 
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,
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:
+                random = ~ 1 | bear, data = learning, family = binomial)
+                random = ~ 1 | bear, data = learning, family = binomial)

Individually, I get results from anova():
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
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
A.J. Rossini
# 
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:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Spencer Graves
    

    

      
      #
    

  


  

      
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. 

      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:

        

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Deepayan Sarkar
    

    

      
      #
    

  


  

      
On Friday 26 November 2004 12:35, Spencer Graves wrote:

            

      
      
      
    

        
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:

            

      
      
      
    

        
+              family = binomial, data = bacteria)

            

      
      
      
    

        
+              family = binomial, data = bacteria)

            

      
      
      
    

        
`log Lik.' -104.3780 (df=5)

            

      
      
      
    

        
`log Lik.' -97.68162 (df=6)

Deepayan

  
  


  


      

    

    

      
      

        


  

        

      Spencer Graves
    

    

      
      #
    

  


  

      
Hi, Deepayan: 

      Thanks much.  That's very helpful. 

ANDREW:  How difficult would it be for you to generate a Monte Carlo to 
simulate data according to your two models, e.g., "y ~ trt" and "y ~ trt 
+ I(week > 2)"?  If you did that, you could then fit both models to each 
set of simulated data, and compute and store logLR <- 
fit2$logLik-fit1$logLik for each one.  This will give you a reference 
distribution, from which you can estimate both the p-value and the 
statistical power of this analysis against your chosen alternatives. 

      If you do that, I suggest you use "GLMM" in library(lme4), not 
glmmPQL.  The "logLik" produced by glmmPLQ for model 2 was LESS THAN 
that for model 1.  If the function were maximizing the likelihood, 
fit2$logLik should be greater than fit1$logLik, not less. 

      Of course, of you simulate both models and compute the 
distribution of your favorite test statistic, you can get a p-value that 
is as good as your simulation.  I've done this kind of thing before, and 
it should be relatively easy in R using rnorm and rbinom. 

      hope this helps. 
      spencer graves

        
Deepayan Sarkar wrote:

        

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Douglas Bates
    

    

      
      #
    

  


  

      
Spencer Graves wrote:

            

      
      
      
    

        
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.

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      John Fox
    

    

      
      #
    

  


  

      
Dear Doug,

I assume that in the absence of a suitable anova() method, there's nothing
wrong with testing the difference between two GLMM objects nested in their
random effects by manually comparing the log-likelihoods as printed or
returned by logLik(). Is that correct?

Regards,
 John

--------------------------------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox 
--------------------------------

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Douglas Bates
    

    

      
      #
    

  


  

      
John Fox wrote:

            

      
      
      
    

        
Yes.

I should have thought of that myself.  At this time of year those in 
Canada where they celebrate Thanksgiving in October are obviously 
thinking more clearly than those in the United States whose minds are 
fogged by too much turkey and televised football games.