AIC, anova error: Models are not all fitted to the same number of observations, models were not all fitted to the same size of dataset

Question

I have models like this:

require(nlme)

set.seed(123)
n <- 100
k <- 5
cat <- as.factor(rep(1:k, n))
cat_i <- 1:k # intercept per kategorie
x <- rep(1:n, each = k)
sigma <- 0.2
alpha <- 0.001
y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma)
plot(x, y)

m1 <- lm(y ~ x)
summary(m1)

m2 <- lm(y ~ cat + x)
summary(m2)

m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit)
summary(m3)

Now I am trying to assess whether the random effect should be present in the model. So I compare the models using AIC or anova, and I get the following error:

> AIC(m1, m2, m3)
   df       AIC
m1  3 1771.4696
m2  7 -209.1825
m3  4 -154.0245
Warning message:
In AIC.default(m1, m2, m3) :
  models are not all fitted to the same number of observations  
> anova(m2, m3)
Error in anova.lmlist(object, ...) : 
  models were not all fitted to the same size of dataset

As you can see, in both cases I use the same dataset. I have found two remedies, but I don't consider them satisfying:

1. Adding method = "ML" to the lme() call - not sure if it is good idea to change the method.
2. Using lmer() instead. Surprisingly, this works, despite the fact that lmer() uses REML method. However I dont like this solution because the lmer() doesn't show p-values for coefficients - I like to use older lme() instead.

Do you have any idea if this is a bug or not and how can we go around that?

Answer 1

A quick search shows that it is possible (although I have to admit that I thought it wasn't) and that it isn't a bug...just another case where methods in R are hidden and result in things that seem 'unexpected', but the RTFM crowd say, "It is in the documentation." Anyway...your solution is to do anova with the lme as the first argument and the lm models as the second (and third if you like) argument(s). If this seems odd, it is because it is a little odd. The reason is that when you call anova, the anova.lme method is called only if the first argument is an lme object. Otherwise, it calls anova.lm (which in turn calls anova.lmlist; if you dig into anova.lm, you'll see why). For details about how you want to be calling anova in this case, pull up help for anova.lme. You'll see that you can compare other models with lme models, but they have to be in a position other than the first argument. Apparently it is also possible to use anova on models fit using the gls function without caring too much about the order of the model arguments. But I don't know enough of the details to determine whether that is a good idea or not, or what exactly it implies (it seems probably fine, but your call). From that link comparing lm to lme appears to be well documented and cited as a method, so I'd err in that direction, were I you.

Good luck.

Answer 2

This is peculiar definitely. As a first thought: when doing model comparison where models are having different fixed effects structures (m2 and m3 for example), it is best to us $ML$ as REML will "change" $y$. (It will multiply it with $k$, where $kX= 0$) That is interesting it that it works using method="ML" which makes me believe it might not be a bug. It seems almost like it enforces "good practice".

Having said that, let's look under the hood:

 methods(AIC)  
 getAnywhere('AIC.default')

 A single object matching ‘AIC.default’ was found
 It was found in the following places
   registered S3 method for AIC from namespace stats
   namespace:stats with value

 function (object, ..., k = 2) 
 {
     ll <- if ("stats4" %in% loadedNamespaces()) 
         stats4:::logLik
     else logLik
     if (!missing(...)) {
         lls <- lapply(list(object, ...), ll)
         vals <- sapply(lls, function(el) {
             no <- attr(el, "nobs") #THIS IS THE ISSUE!
             c(as.numeric(el), attr(el, "df"), if (is.null(no)) NA_integer_ else no)
         })
         val <- data.frame(df = vals[2L, ], ll = vals[1L, ])
         nos <- na.omit(vals[3L, ])
         if (length(nos) && any(nos != nos[1L])) 
             warning("models are not all fitted to the same number of observations")
         val <- data.frame(df = val$df, AIC = -2 * val$ll + k * val$df)
             Call <- match.call()
             Call$k <- NULL
         row.names(val) <- as.character(Call[-1L])
         val
     }
     else {
         lls <- ll(object)
         -2 * as.numeric(lls) + k * attr(lls, "df")
     }     
 }

where in your case you can see that :

  lls <- lapply(list(m2,m3), stats4::logLik)
  attr(lls[[1]], "nobs")
  #[1] 500
  attr(lls[[2]], "nobs")
  #[1] 498

and obviously logLik is doing something (maybe?) unexpected...? no, not really, if you look at the documentation of logLik, ?logLik, you'll see it is explicitly stated:

 There may be other attributes depending on the method used: see
 the appropriate documentation.  One that is used by several
 methods is ‘"nobs"’, the number of observations used in estimation
 (after the restrictions if ‘REML = TRUE’)

which brings us back to our original point, you should be using ML.

To use a common saying in CS: "It's not a bug; it's an (real) feature!"

EDIT: (Just to address your comment:) Assume you fit another model using lmer this time:

m3lmer <- lmer(y ~ x + 1|cat)

and you do the following:

lls <- lapply(list(m2,m3, m3lmer), stats4::logLik)
attr(lls[[3]], "nobs")
#[1] 500
 attr(lls[[2]], "nobs")
#[1] 498

Which seems like a clear discrepancy between the two but it really isn't as Gavin explained. Just to state the obvious:

 attr( logLik(lme(y ~ x, random = ~ 1|cat, na.action = na.omit, method="ML")),
 "nobs")
#[1] 500

There is a good reason why this happens in terms of methodology I think. lme does try to make sense of the LME regression for you while lmer when doing model comparisons it falls back immediately to it's ML results. I think there are no bugs on this matter in lme and lmer just different rationales behind each package.

See also Gavin Simposon's comment on a more insightful explanation of what went on with anova() (The same thing practically happens with AIC)