[R] lmer and glmmPQL

Wed Dec 7 16:16:37 CET 2005

Thanks for the reply Doug!

A follow up question and comment ...

1) If I understand correctly, looking at a simple situation in which
SITES are nested in ZONES, the following should be similar.  However,
despite the same F values, the p-value from lmer is 1/2 the other
methods.  Why is this true?

> anova(lmer(RICH ~ ZONE + (1|SITE:ZONE), data))
Analysis of Variance Table
     Df Sum Sq Mean Sq  Denom F value  Pr(>F)  
ZONE  4   97.8    24.5 6610.0  2.1046 0.07753 .

> # make the nesting explicit
> data$SinZ = with(data, ZONE:SITE)[drop=TRUE]
> anova(lme(RICH ~ ZONE, data, random = ~1 | SinZ))
            numDF denDF   F-value p-value
(Intercept)     1  6600 100.38331  <.0001
ZONE            4    10   2.10459   0.155

> summary(aov(RICH ~ ZONE + Error(SITE:ZONE), data))

Error: SITE:ZONE
          Df Sum Sq Mean Sq F value Pr(>F)
ZONE       4  29669    7417  2.1046  0.155
Residuals 10  35243    3524   

2) I think the anova problems with lmer may also apply to poisson.
Compare the following (which includes a covariate).  Based on the
parameter estimates, the covariate should be significant.  However, its
anova p-value is .998:

> lmer(RICH ~ ZONE + lANPP + (1|SITE:ZONE), family = poisson, data)
Generalized linear mixed model fit using PQL 
Formula: RICH ~ ZONE + lANPP + (1 | SITE:ZONE) 
   Data: data 
 Family: poisson(log link)
      AIC      BIC    logLik deviance
 9700.252 9754.628 -4842.126 9684.252
Random effects:
     Groups        Name    Variance    Std.Dev. 
  SITE:ZONE (Intercept)    0.069493     0.26361 
# of obs: 6615, groups: SITE:ZONE, 15

Estimated scale (compare to 1)  1.183970 

Fixed effects:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.5169605  0.1533564  9.8917  < 2e-16 ***
ZONE2        0.4034169  0.2156956  1.8703  0.06144 .  
ZONE3       -0.1772011  0.2158736 -0.8209  0.41173    
ZONE4       -0.2368290  0.2158431 -1.0972  0.27254    
ZONE5       -0.1011186  0.2158114 -0.4686  0.63939    
lANPP        0.2201926  0.0081857 26.8995  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

> anova(lmer(RICH ~ ZONE + lANPP + (1|SITE:ZONE), family = poisson,
data))
Analysis of Variance Table
      Df    Sum Sq   Mean Sq Denom   F value Pr(>F)
ZONE   4 2.809e-05 7.022e-06  6609 4.298e-06 1.0000
lANPP  1 5.229e-06 5.229e-06  6609 3.200e-06 0.9986

Thanks again for any insight you may be able to provide!!

> -----Original Message-----
> From: Douglas Bates [mailto:dmbates at gmail.com] 
> Sent: Wednesday, December 07, 2005 8:28 AM
> To: Cox, Stephen
> Cc: r-help at stat.math.ethz.ch
> Subject: Re: [R] lmer and glmmPQL
> 
> On 12/5/05, Cox, Stephen <stephen.cox at ttu.edu> wrote:
> > I have been looking into both of these approaches to conducting a 
> > GLMM, and want to make sure I understand model 
> specification in each.  
> > In particular - after looking at Bates' Rnews article and searching 
> > through the help archives, I am unclear on the 
> specification of nested 
> > factors in lmer.  Do the following statements specify the same mode 
> > within each approach?
> >
> > m1 = glmmPQL(RICH ~ ZONE, family = poisson, data, random = ~ YEAR | 
> > SITE / QUADRAT)
> > m2 = lmer(RICH ~ ZONE +(YEAR|SITE)+ (YEAR|QUADRAT), family 
> = poisson,
> > data)
> 
> If you want to ensure that QUADRAT is nested within SITE then 
> use the interaction operator explicitly
> 
> m2 <- lmer(RICH ~ ZONE +(YEAR|SITE)+ (YEAR|SITE:QUADRAT), 
> family = poisson, data)
> 
> For the grouping factors nested versus non-nested depends on 
> the coding.  If QUADRAT has a distinct level for each 
> SITE:QUADRAT combination then the nesting will automatically 
> be detected.  However, if the nesting is implicit (that is, 
> if levels of QUADRAT are repeated at different SITES) then it 
> is necessary to use the interaction operator.  There is no 
> harm in using the interaction operator when the nesting is explicit.
> >
> > As a follow up - what would be the most appropriate model formula 
> > (using glmmPQL syntax) to specify both a nested facor and repeated 
> > observations?  Specifically, I am dealing with experimental 
> data with 
> > three factors.  ZONE is a fixed effect.  Three sites (SITE) 
> are nested 
> > within each ZONE.  Multiple quadrats within each SITE are measured 
> > across multiple years.  I want to represent the nesting of 
> SITE within 
> > ZONE and allow for repeated observations within each 
> QUADRAT over time 
> > (the YEAR | QUADRAT random effect).  -- I am assuming that 
> glmmPQL is 
> > the best option at this point because of recent discussion on Rhelp 
> > about issues associated with the Matrix package used in lmer (i.e., 
> > the anova results do not seem to match parameter tests).
> >
> 
> I believe the anova problems only occur with a binomial response. 
> They are caused by my failure to use the prior.weights appropriately. 
> For a Poisson model this should not be a problem.
> 
> > Any information would be very much appreciated!
> >
> > Regards
> >
> > Stephen
> >
> > ______________________________________________
> > 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
> >
> 
>