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

Kingsford Jones kingsfordjones at gmail.com
Mon Sep 7 22:26:43 CEST 2009


On Mon, Sep 7, 2009 at 10:34 AM, Alain Zuur<highstat at highstat.com> wrote:
>
>
>
> rapton wrote:
>>
>> 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.

And, unless I'm missing something, also by the (misspecified) test.  A
large p-value indicates you have no evidence that the additional 19
parameters in f2 improve fit, which matches what the other methods
suggested.  However, as has been pointed out, the lack of nesting
makes this a faulty LRT.  This is made apparent by the fact that you
get a test statistic outside the support of the chi-squared
distribution (positive reals)

> (lambda <- (-2)*(-22715 - (-23633)))
[1] -1836

and since the test is uses right-tail probability, anova is not
changing anything by moving the statistic to 0.

> pchisq(lambda, 19, lower = FALSE)
[1] 1
> pchisq(0, 19, lower = FALSE)
[1] 1

To do the test properly the restricted (null) model must be a special
case of the general (alternative) model (e.g., with the additional 19
parameters set to zero) which will result in the null model having a
smaller likelihood, leading to a positive tests statistic.  When that
statistic is small you get a large p-value indicating a lack of
evidence that the additional parameters improve fit...

hth,

Kingsford




 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
>>
>
>
> ** NOT  ** nested       ....sorry....the brain is going faster than the
> fingers.
>
>
>
>
>
> -----
> --------------------------------------------------------------------
> Dr. Alain F. Zuur
> First author of:
>
> 1. Analysing Ecological Data (2007).
> Zuur, AF, Ieno, EN and Smith, GM. Springer. 680 p.
>
> 2. Mixed effects models and extensions in ecology with R. (2009).
> Zuur, AF, Ieno, EN, Walker, N, Saveliev, AA, and Smith, GM. Springer.
>
> 3. A Beginner's Guide to R (2009).
> Zuur, AF, Ieno, EN, Meesters, EHWG. Springer
>
>
> Statistical consultancy, courses, data analysis and software
> Highland Statistics Ltd.
> 6 Laverock road
> UK - AB41 6FN Newburgh
> Email: highstat at highstat.com
> URL: www.highstat.com
>
>
>
> --
> View this message in context: http://www.nabble.com/Using-anova%28f1%2C-f2%29-to-compare-lmer-models-yields-seemingly-erroneous-Chisq-%3D-0%2C-p-%3D-1-tp25297254p25333148.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>




More information about the R-help mailing list