[R] Using anova(f1, f2) to compare lmer models yields seemingly erroneous Chisq = 0, p = 1

Viechtbauer Wolfgang (STAT) Wolfgang.Viechtbauer at STAT.unimaas.nl
Mon Sep 7 10:23:41 CEST 2009


Are the two models (f1 and f2) actually nested?

Aside from that, it is strange that the output is exactly the same after you used REML=FALSE. The log likelihoods should have changed.

Best,

--
Wolfgang Viechtbauer
 Department of Methodology and Statistics
 School for Public Health and Primary Care
 University of Maastricht, The Netherlands
 http://www.wvbauer.com/




----Original Message----
From: r-help-bounces at r-project.org
[mailto:r-help-bounces at r-project.org] On Behalf Of Matt Killingsworth
Sent: Friday, September 04, 2009 22:29 To: Bert Gunter
Cc: r-help at r-project.org
Subject: Re: [R] Using anova(f1, f2) to compare lmer models yields
seemingly erroneous Chisq = 0, p = 1

> Hi Bert,
>
> Thank you for your note!  I tried changing the REML default, and it
> still produces the same result (see below).  Is that what you meant
> for me to try?
>
> Incidentally, I am using lmer() not lme()
>
>  ### ORIGINAL ###
>> f1 <- (lmer(outcome ~ predictor.1 + (1 | person), data=i))
>> f2 <- (lmer(outcome ~ predictor.2 + (1 | person), data=i)) anova(f1,
>> f2)
>
> Data: i
> Models:
> f1: outcome ~ predictor.1 + (1 | person)
> f2: outcome ~ predictor.2 + (1 | person)
>     Df    AIC     BIC    logLik Chisq Chi Df Pr(>Chisq)
> f1  6  45443  45489 -22715
> f2 25  47317  47511 -23633     0     19          1
>
> ### DO NOT USE REML ###
>> f1 <- (lmer(outcome ~ predictor.1 + (1 | person), data=i, REML =
>> FALSE)) f2 <- (lmer(outcome ~ predictor.2 + (1 | person), data=i,
>> REML = FALSE)) anova(f1, f2)
>
> Data: i
> Models:
> f1: outcome ~ predictor.1 + (1 | person)
> f2: outcome ~ predictor.2 + (1 | person)
>     Df    AIC     BIC    logLik Chisq Chi Df Pr(>Chisq)
> f1  6  45443  45489 -22715
> f2 25  47317  47511 -23633     0     19          1
>
>
> On Fri, Sep 4, 2009 at 4:18 PM, Bert Gunter <gunter.berton at gene.com>
> wrote:
>
>> My guess would be:
>>
>> "Likelihood comparisons are not meaningful for objects fit using
>> restricted maximum likelihood and with different fixed effects. "
>>
>> (from ?anova.lme in the nlme package).
>>
>> Are you using the REML = TRUE default?
>>
>> Bert Gunter
>> Genentech Nonclinical Statistics
>>
>> -----Original Message-----
>> From: r-help-bounces at r-project.org
>> [mailto:r-help-bounces at r-project.org]
>> On
>> Behalf Of rapton
>> Sent: Friday, September 04, 2009 9:10 AM
>> To: r-help at r-project.org
>> Subject: [R] Using anova(f1, f2) to compare lmer models yields
>> seemingly erroneous Chisq = 0, p = 1
>>
>>
>> Hello,
>>
>> I am using R to analyze a large multilevel data set, using
>> lmer() to model my data, and using anova() to compare the fit of
>> various models.  When I run two models, the output of each model is
>> generated correctly as far as I can tell (e.g. summary(f1) and
>> summary(f2) for the multilevel model output look perfectly
>> reasonable), and in this case (see
>> below) predictor.1 explains vastly more variance in outcome than
>> predictor.2 (R2 = 15% vs. 5% in OLS regression, with very large N).
>> What I am utterly puzzled by is that when I run an anova comparing
>> the two multilevel model fits, the Chisq comes back as 0, with p =
>> 1.  I am pretty sure that fit #1 (f1) is a much better predictor of
>> the outcome than f2, which is reflected in the AIC, BIC , and logLik
>> values.  Why might anova be giving me this curious output?  How can
>> I fix it?  I am sure I am making a dumb error somewhere, but I
>> cannot figure out what it is.  Any help or suggestions would be
>> greatly appreciated!
>>
>> -Matt
>>
>>
>>> f1 <- (lmer(outcome ~ predictor.1 + (1 | person), data=i)) f2 <-
>>> (lmer(outcome ~ predictor.2 + (1 | person), data=i)) anova(f1, f2)
>>
>> Data: i
>> Models:
>> f1: outcome ~ predictor.1 + (1 | person)
>> f2: outcome ~ predictor.2 + (1 | person)
>>   Df    AIC      BIC    logLik   Chisq Chi Df Pr(>Chisq)
>> f1  6  45443  45489 -22715
>> f2 25  47317  47511 -23633     0     19          1
>> --




More information about the R-help mailing list