[R] lmer, p-values and all that

[Ben Bolker]

On 13-03-27 10:10 PM, David Winsemius wrote:
> 
> On Mar 27, 2013, at 7:00 PM, Ben Bolker wrote:
> 
>> Michael Grant <michael.grant <at> colorado.edu> writes:
>>> Dear Help:
>>
>>> I am trying to follow Professor Bates' recommendation, quoted by
>>> Professor Crawley in The R Book, p629, to determine whether I should
>>> model data using the 'plain old' lm function or the mixed model
>>> function lmer by using the syntax anova(lmModel,lmerModel).
>>> Apparently I've not understood the recommendation or the proper
>>> likelihood ratio test in question (or both) for I get this error
>>> message: Error: $ operator not defined for this S4 class.
>>
>>  I don't have the R Book handy (some more context would be extremely
>> useful!  I would think it would count as "fair use" to quote the
>> passage you're referring to ...)
> 
> This is the quoted Rhelp entry:
> 
> http://tolstoy.newcastle.edu.au/R/help/05/01/10006.html
> 
> (I'm unable to determine whether it applies to the question at hand.)

  OK, I misspoke -- sorry.  I think the lmer()/lme() likelihoods are
actually comparable; it's GLMMs (glmer(), with no analogue in lme()-land
except
for MASS::glmmPQL, which doesn't give you log-likelihoods at all)
where the problem arises.

  You can (1) use lme(), (2)  look at http://glmm.wikidot.com/faq for
suggestions about testing random-effects terms (including "perhaps
don't do it at all"), or (3) construct the likelihood ratio test
yourself as follows:

library("nlme")
data(Orthodont)
fm1 <- lme(distance~age,random=~1|Subject,data=Orthodont)
fm0 <- lm(distance~age,data=Orthodont)
anova(fm1,fm0)[["p-value"]]
detach("package:nlme",unload=TRUE)
library("lme4.0") ## stable version of lme4
fm2 <- lmer(distance~age+(1|Subject),data=Orthodont,REML=FALSE)
anova(fm2,fm0) ## fails
ddiff <- -2*c(logLik(fm0)-logLik(fm2))
pchisq(ddiff,1,lower.tail=FALSE)
## not identical to above but close enough

> 
>>
>>> Would someone be kind enough to point out my blunder?
>>
>>  You should probably repost this to the r-sig-mixed-models at r-project.org
>> mailing list.
>>
>>  My short answer would be: (1) I don't think you can actually
>> use anova() to compare likelihoods between lm() and lme()/lmer()
>> fits in the way that you want: *maybe* for lme() [don't recall],
>> but almost certainly not for lmer().  See http://glmm.wikidot.com/faq
>> for methods for testing significance/inclusion of random factors
>> (short answer: you should *generally* try to make the decision
>> whether to include random factors or not on _a priori_ grounds,
>> not on the basis of statistical tests ...)
>>
>>  Ben Bolker
>>
>