R help - Mar 2012 - AIC models are not all fitted to the same number of observation

Hi,

Using lme from the package nlme  3.1-103, I meet a strange warning. I am 
trying to compare to models with:

library(nlme)
lmez6=lme(lepus~vulpes,random=~1|troncon/an,data=ika_z6_test)
lmez60=lme(lepus~1,random=~1|troncon/an,data=ika_z6_test)

Both have the same number of observations and groups:

lmez6
Linear mixed-effects model fit by REML
   Data: ika_z6_test
   Log-restricted-likelihood: -2267.756
   Fixed: lepus ~ vulpes
(Intercept)      vulpes
  1.35017117  0.04722338

Random effects:
  Formula: ~1 | troncon
         (Intercept)
StdDev:   0.8080261

  Formula: ~1 | an %in% troncon
         (Intercept)  Residual
StdDev:    1.086611 0.4440076

Number of Observations: 1350
Number of Groups:
         troncon an %in% troncon
             169            1350


 > lmez60
Linear mixed-effects model fit by REML
   Data: ika_z6_test
   Log-restricted-likelihood: -2266.869
   Fixed: lepus ~ 1
(Intercept)
    1.435569

Random effects:
  Formula: ~1 | troncon
         (Intercept)
StdDev:   0.8139646

  Formula: ~1 | an %in% troncon
         (Intercept)  Residual
StdDev:    1.086843 0.4445815

Number of Observations: 1350
Number of Groups:
         troncon an %in% troncon
             169            1350

...but when I want to compare their AIC, I get:

AIC(lmez6,lmez60)
        df      AIC
lmez6   5 4545.511
lmez60  4 4541.737
Warning message:
In AIC.default(lmez6, lmez60) :
   models are not all fitted to the same number of observations


Has anybody an explanation about this strange warning ? To what extent 
this warning may limit the conclusions that could be drawn from AIC 
comparison ?

Thanks in advance,

Patrick

Le 21/03/2012 10:56, Patrick Giraudoux a ?crit :
> Hi,
>
> Using lme from the package nlme  3.1-103, I meet a strange warning. I 
> am trying to compare to models with:
>
> library(nlme)
> lmez6=lme(lepus~vulpes,random=~1|troncon/an,data=ika_z6_test)
> lmez60=lme(lepus~1,random=~1|troncon/an,data=ika_z6_test)
>
> Both have the same number of observations and groups:
>
> lmez6
> Linear mixed-effects model fit by REML
>   Data: ika_z6_test
>   Log-restricted-likelihood: -2267.756
>   Fixed: lepus ~ vulpes
> (Intercept)      vulpes
>  1.35017117  0.04722338
>
> Random effects:
>  Formula: ~1 | troncon
>         (Intercept)
> StdDev:   0.8080261
>
>  Formula: ~1 | an %in% troncon
>         (Intercept)  Residual
> StdDev:    1.086611 0.4440076
>
> Number of Observations: 1350
> Number of Groups:
>         troncon an %in% troncon
>             169            1350
>
>
> > lmez60
> Linear mixed-effects model fit by REML
>   Data: ika_z6_test
>   Log-restricted-likelihood: -2266.869
>   Fixed: lepus ~ 1
> (Intercept)
>    1.435569
>
> Random effects:
>  Formula: ~1 | troncon
>         (Intercept)
> StdDev:   0.8139646
>
>  Formula: ~1 | an %in% troncon
>         (Intercept)  Residual
> StdDev:    1.086843 0.4445815
>
> Number of Observations: 1350
> Number of Groups:
>         troncon an %in% troncon
>             169            1350
>
> ...but when I want to compare their AIC, I get:
>
> AIC(lmez6,lmez60)
>        df      AIC
> lmez6   5 4545.511
> lmez60  4 4541.737
> Warning message:
> In AIC.default(lmez6, lmez60) :
>   models are not all fitted to the same number of observations
>
>
> Has anybody an explanation about this strange warning ? To what extent 
> this warning may limit the conclusions that could be drawn from AIC 
> comparison ?
>
> Thanks in advance,
>
> Patrick
>
>
>
Sorry to go on on the thread, I have created, but the trouble I meet is 
above my level in stats... Actually, not using AIC but an anova 
approach, I get a more informative message:

anova(lmez6, lmez60)
        Model df      AIC      BIC    logLik   Test  L.Ratio p-value
lmez6      1  5 4545.511 4571.543 -2267.756
lmez60     2  4 4541.737 4562.566 -2266.869 1 vs 2 1.774036  0.1829
Warning message:
In anova.lme(lmez6, lmez60) :
   Fitted objects with different fixed effects. REML comparisons are not 
meaningful.

And fubbling a bit more, I disclosed that this was an effect of fitting 
the model using REML. If fitted using ML, things are going (apparently) 
smoothly:

lmez6=lme(lepus~vulpes,random=~1|troncon/an,data=ika_z6_test,method="ML")
 >
lmez60=lme(lepus~1,random=~1|troncon/an,data=ika_z6_test,method="ML")
 > anova(lmez6, lmez60)
        Model df      AIC      BIC    logLik   Test  L.Ratio p-value
lmez6      1  5 4536.406 4562.445 -2263.203
lmez60     2  4 4538.262 4559.093 -2265.131 1 vs 2 3.856102  0.0496

 > AIC(lmez6,lmez60)
        df      AIC
lmez6   5 4536.406
lmez60  4 4538.262

Now I have the following problem. What I understood from Pinheiro and 
Bates's book and some forums, is that ML estimations are biased to some 
extent tending to underestimate variance parameters. So probably not to 
recommend however results looks consistent here.

Thus, I am lost. The two models looks to me clearly embedded (one is 
just a null model with the only intercept to estimate and the other with 
intercept + one independent variable (numeric), both have the same 
random effects, the same response variable and the same number of 
observations). Warnings, from this point of view sounds inconsistent. 
They are probably not, but beyond my understanding...

Any idea ?

Patrick