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
)