[R-sig-ME] anova (lm, lmer ) question

[Ben Bolker]

On 14-10-06 12:36 PM, Farrar, David wrote:
> If I understood, the following should have worked?
> 
>> m.ML <- lmer(
> + log10(TWA) ~ ns(Wind.Speed,3) + isLamar + isLime1p + isLime5p  + log10(TWAuw)
> + + (1|BiosolidSource) + (1|sample) + (1|sample.trial),
> + REML=F,
> + data=da.regr
> + )
>> m.lm <- lm(
> + log10(TWA) ~ ns(Wind.Speed,3) + isLamar + isLime1p + isLime5p + log10(TWAuw),
> + data=da.regr
> + )
>> anova(m.ML, m.lm)
> Error in UseMethod("isREML") : 
>   no applicable method for 'isREML' applied to an object of class "lm"

  I didn't mean to imply that anova(m.ML,m.lm) would actually work, but
rather that the equivalent calculation would be appropriate, e.g.

    library(lme4)
    fm1 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy, REML=FALSE)
    fm0 <- lm(Reaction ~ Days, sleepstudy)
    NL1 <- -logLik(fm1)
    NL0 <- -logLik(fm0)
    devdiff <- 2*(NL0-NL1)
    dfdiff <- attr(NL1,"df")-attr(NL0,"df")
    pchisq(devdiff,dfdiff,lower.tail=FALSE)

The p-value is very small in this case, but that's consistent
with a large/well-determined variance estimate ...

    pp <- profile(fm1)
    library(lattice)
    xyplot(logProf(pp))

Keep in mind that the likelihood ratio test also has some theoretical
problems in this case, mainly with boundary issues (see
http://glmm.wikidot.com/faq for more info)

> 
> -----Original Message-----
> From: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of Ben Pelzer
> Sent: Saturday, October 04, 2014 10:13 AM
> To: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] anova (lm, lmer ) question
> 
> Dear romunov, Ben and Ken,
> 
> Thanks for your replies. From these I conclude that:
> - for linear (lmer vs. lm) models there's no problem in using the deviance difference
> - for generalized  linear models (glmer vs. glm) it's ok to use the deviance difference as long as nAGQ=1.
> Would you agree with me? Best regards,
> 
> Ben.
> 
> On 4-10-2014 2:48, Ben Bolker wrote:
>> Thanks for checking.  The comparison with Stata isn't necessarily 
>> relevant though -- or question is whether `lm` and `lmer` (or `glm` 
>> and `glmer`) include/exclude the same additive constants, so that 
>> their log-likelihoods are directly comparable.
>>
>> On Fri, Oct 3, 2014 at 8:38 PM, Ken Beath <ken.beath at mq.edu.au> wrote:
>>
>>> nAGQ=1 and greater than 1 give different results, and the nAGQ=1 
>>> matches fairly closely the log likelihood from Stata for 3 quadrature 
>>> points, so presumably is correct. Stata's Laplace didn't converge with my data.
>>>
>>>
>>> Ken
>>>
>>>
>>>
>>> On 4 October 2014 09:06, Ben Bolker <bbolker at gmail.com> wrote:
>>>
>>>> romunov <romunov at ...> writes:
>>>>
>>>>> FWIW, this is from the glmm faq site <http://glmm.wikidot.com/faq>.
>>>>>
>>>>> How can I test whether a random effect is significant?
>>>>>
>>>>    ...
>>>>
>>>>>     - *do not* compare lmer models with the corresponding lm fits, or
>>>>>     glmer/glm; the log-likelihoods are not commensurate (i.e., they
>>>> include
>>>>>     different additive terms)
>>>>    For what it's worth, I believe this is out of date, _except_ for 
>>>> glmer fits with nAGQ>1.  It should be possible to implement
>>>> anova(<merMod>,<lm>/<glm>) -- it's only a nuisance (sadly, if we 
>>>> were still using S4 classes at this level it would be easier ...)
>>>>
>>>>    Ben Bolker
>>>>
>>>> _______________________________________________
>>>> R-sig-mixed-models at r-project.org mailing list 
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>>
>>>
>>>
>>> --
>>>
>>> *Ken Beath*
>>> Lecturer
>>> Statistics Department
>>> MACQUARIE UNIVERSITY NSW 2109, Australia
>>>
>>> Phone: +61 (0)2 9850 8516
>>>
>>> Building E4A, room 526
>>> http://stat.mq.edu.au/our_staff/staff_-_alphabetical/staff/beath,_ken
>>> /
>>>
>>> CRICOS Provider No 00002J
>>> This message is intended for the addressee named and 
>>> m...{{dropped:11}}
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list 
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>