Re: [R] lmer and glmmPQL

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Sun 11 Dec 2005 - 11:33:00 EST

Hello, Stephen:

          I can't match Doug's understanding of these issues, but I will offer a couple of comments. Regarding your question 1, the difference is in the denominator of degrees of freedom: 6620 for lmer and 10 for lme.

  pf(2.10459, 4, c(6610, 10), lower.tail=FALSE) [1] 0.07752865 0.15499674

          This may not tell you more than you already know. However, from earlier remarks on the listserve, I know that getting accurate p values for these kinds of models is still a substantive research issue. If you want an accurate p value, there may be no substitute for Monte Carlo. I just got 34 hits for 'RSiteSearch("simulate.lme")' and 10 for 'RSiteSearch("simulate lmer")'. If you are not already familiar with Pinheiro and Bates (2000) Mixed-Effects Models in S and S-Plus (Springer), I suspect you will find at least partial answers to many of your mixed-effects questions there. Certainly, I found it quite illuminating.

	  Best Wishes,
	  spencer graves

Cox, Stephen wrote:

> 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@gmail.com] 
>>Sent: Wednesday, December 07, 2005 8:28 AM
>>To: Cox, Stephen
>>Cc: r-help@stat.math.ethz.ch
>>Subject: Re: [R] lmer and glmmPQL
>>
>>On 12/5/05, Cox, Stephen <stephen.cox@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@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
>>>
>>
>>

>
> ______________________________________________
> R-help@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
-- 
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

spencer.graves@pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915

______________________________________________
R-help@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
Received on Sun Dec 11 11:40:53 2005

This archive was generated by hypermail 2.1.8 : Sun 11 Dec 2005 - 14:50:56 EST