Skip to content

Questions on the results from glmmPQL(MASS)

8 messages · Ken Beath, zhijie zhang, David Duffy

Original
Ken Beath
#
On 07/12/2008, at 12:34 PM, zhijie zhang wrote:

            

      
      
      
    

        
PQL is not maximum likelihood (it is an approximation which uses lme  
internally and this is what generates the results) so the results  
should be NA. The R and S-Plus versions have obviously diverged.

SAS nlmixed uses maximum likelihood so a log likelihood is available.

            

      
      
      
    

        
I wonder if in PQL these have any real meaning, as they are obtained  
from the linear mixed effects step. Try using lmer in the lme4 package.

            

      
      
      
    

        
The standard answer is of the form "Just because SAS has it, doesn't  
mean it is a good idea".  There was a recent discussion on this list  
about significance of random effects.

Ken

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    
1 day later

    

      
      

        
      

    

    

      
      

        


  

        

      David Duffy
    

    

      
      #
    

  


  

      
On Sun, 7 Dec 2008, zhijie zhang wrote:

        

            

      
      
      
    

        
Someone answered your specific questions before.
Because this is a simple random intercept model, you can also use the
glmmML R package, which gives:

glmmML(y ~ trt + I(week > 2), cluster=ID, family=binomial, data=bacteria)

                    coef se(coef)      z Pr(>|z|)
(Intercept)      3.5482   0.6933  5.118 3.09e-07
trtdrug         -1.3667   0.6768 -2.020 4.34e-02
trtdrug+        -0.7826   0.6830 -1.146 2.52e-01
I(week > 2)TRUE -1.5986   0.4754 -3.363 7.71e-04

Standard deviation in mixing distribution:  1.242
Std. Error:                                 0.4024

Residual deviance: 192.3 on 215 degrees of freedom      AIC: 202.3

You can then construct a test for the random effect.

  
    

      
      
    

  


  


      

    

    

      
      

        
      

    

    

      
      

        


  

        

      David Duffy
    

    

      
      #
    

  


  

      
On Mon, 8 Dec 2008, zhijie zhang wrote:

        

            

      
      
      
    

        
anova does not have a method for glmmML, but the deviances seem to be 
calculated the same (see model0$cluster.null.deviance etc):.

model0 Residual deviance: 192.3  on 215 degrees of freedom      AIC: 202.3
model1 Residual deviance: 199.18 on 216  degrees of freedom

LRTS = 6.88.  We will assume that the test statistic is distributed 1/2 
X2(0) and 1/2 X2(1), so P ~ 0.004.  Comparing to a Wald test using the SE 
on the random effect SD, I get:

Z = 1.242/0.4024 = 3.08, P=0.001

  
    

      
      
    

  


  


      

    

    

      
      

        
      

    

    

      
      

        


  

        

      Ken Beath
    

    

      
      #
    

  


  

      
On 08/12/2008, at 6:42 PM, zhijie zhang wrote:

        

            

      
      
      
    

        
X2 is chi-square. Because the test is on the boundary the test  
statistic is distributed as the weighted sum of chi-square rather than  
the usual chi-square. Verbeke and Molenberghs cover this in their books.

Simulation (parametric bootstrap) seems a better way of doing this.  
Information criteria (AIC or BIC) can also be used. Most usual  
justification for a random effect is that there is expected to be one,  
so provided it can be estimated it is included.

            

      
      
      
    

        
P value for z test is 2*dnorm(3.08) which is close enough to 0.001  
given that the test is only an approximation.

Ken