[R-sig-ME] anova function gives (irrelevant) results on two models with different numbers of observations

[Olivier Renaud]

Hi,

I have encountered what I believe is a bug or an extremely unwanted 
"feature" in lme4. In short, the anova function does not check that the 
two models are based on the same number of observations and can thus 
give irrelevant results, as exemplified below. Note that lme (at least 
on S+) correctly refuses to compare the models.

I found this problem, because I wanted to understand why the LRT gave 
extremely different value than the Wald/t-stat index/statistics/test 
(choose one ;-) ), see below. It took me a while to realize that the 
number of observations was not the same in the two models (156 vs. 160). 
It happens when an explanatory variable is present in only one of the 
two models and contains NA's, which in not an uncommon situation. The 
MLdeviances then are very different but not comparable. I'm afraid it 
might have unnoticed by other users. I have reproduced the calls in lmer 
(R) and lme (S+).

Olivier

R version 2.6.1 (2007-11-26)
Package:       lme4
Version:       0.99875-9
 > (Fu.mod0.lmer <- lmer(Y ~ 1 + ( 1|subjec), data = Fu, 
na.action=na.omit, method="ML") )
Linear mixed-effects model fit by maximum likelihood
Formula: Y ~ 1 + (1 | subjec)
   Data: Fu
   AIC   BIC logLik MLdeviance REMLdeviance
 654.7 660.9 -325.4      650.7        651.4
Random effects:
 Groups   Name        Variance Std.Dev.
 subjec   (Intercept) 0.93595  0.96745
 Residual             2.96389  1.72160
number of obs: 160, groups: subjec, 16

Fixed effects:
            Estimate Std. Error t value
(Intercept)   3.4875     0.2775   12.57

 > (Fu.mod1.lmer <- lmer(Y ~ expvariable + ( 1|subjec), data = Fu, 
na.action=na.omit, method="ML") )
Linear mixed-effects model fit by maximum likelihood
Formula: Y ~ expvariable + (1 | subjec)
   Data: Fu
   AIC   BIC logLik MLdeviance REMLdeviance
 637.8 646.9 -315.9      631.8        633.1
Random effects:
 Groups   Name        Variance Std.Dev.
 subjec   (Intercept) 1.0357   1.0177 
 Residual             2.8796   1.6969 
number of obs: 156, groups: subjec, 16

Fixed effects:
            Estimate Std. Error t value
(Intercept)   3.2673     0.3145   10.39
expvariable   0.4140     0.2856    1.45

Correlation of Fixed Effects:
            (Intr)
expvariable -0.398

 > anova(Fu.mod0.lmer, Fu.mod1.lmer)
Data: Fu
Models:
Fu.mod0.lmer: Y ~ 1 + (1 | subjec)
Fu.mod1.lmer: Y ~ expvariable + (1 | subjec)
             Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)   
Fu.mod0.lmer  2  654.70  660.85 -325.35                            
Fu.mod1.lmer  3  637.78  646.93 -315.89 18.917      1  1.365e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



######## In S+8.0

 > Fu.mod0.lme <- lme(Y ~ 1, random=~1|subjec, data = Fu, 
na.action=na.omit, method="ML")
 > summary(Fu.mod0.lme)
Linear mixed-effects model fit by maximum likelihood
 Data: Fu
       AIC      BIC    logLik
  656.7007 665.9262 -325.3503

Random effects:
 Formula:  ~ 1 | subjec
        (Intercept) Residual
StdDev:   0.9674476 1.721595

Fixed effects: Y ~ 1
             Value Std.Error  DF  t-value p-value
(Intercept) 3.4875 0.2783988 144 12.52699  <.0001

Standardized Within-Group Residuals:
       Min         Q1        Med        Q3      Max
 -2.244476 -0.7371921 -0.2445665 0.8289158 2.608248

Number of Observations: 160
Number of Groups: 16
 > Fu.mod1.lme <- lme(Y ~ expvariable, random=~1|subjec, data = Fu, 
na.action=na.omit, method="ML")
 > summary(Fu.mod1.lme)
Linear mixed-effects model fit by maximum likelihood
 Data: Fu
       AIC     BIC    logLik
  639.7835 651.983 -315.8918

Random effects:
 Formula:  ~ 1 | subjec
        (Intercept) Residual
StdDev:    1.017715 1.696945

Fixed effects: Y ~ expvariable
               Value Std.Error  DF  t-value p-value
(Intercept) 3.267305 0.3164881 139 10.32363  <.0001
expvariable 0.413960 0.2874080 139  1.44032   0.152

Standardized Within-Group Residuals:
     Min         Q1        Med        Q3      Max
 -2.4028 -0.6934754 -0.2075537 0.7346619 2.768138

Number of Observations: 156
Number of Groups: 16
 > anova(Fu.mod0.lme, Fu.mod1.lme)
Problem in anova.lme(Fu.mod0.lme, Fu.mod1.lme): All fitted objects must 
use the same number of observations
Use traceback() to see the call stack

-- 
Olivier.Renaud at pse.unige.ch             http://www.unige.ch/~renaud/
Methodology & Data Analysis - Psychology Dept - University of Geneva
UniMail, Office 4142  -  40, Bd du Pont d'Arve   -  CH-1211 Geneva 4