9
$\begingroup$

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?

$\endgroup$
6
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ Oh, and user11852's answer in regards to AIC with Gavin's addendum stands, there is no special AIC.lme or anything to address that issue and the whole thing starts to slip beyond my pay grade $\endgroup$ – russellpierce Jul 25 '13 at 1:42
3
$\begingroup$

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)

$\endgroup$
  • $\begingroup$ "you should be using ML" - but how can you explain that lmer is using REML (see the model summary) and works fine in AIC? So there are two possibilities: 1) the error message is *a feature, not bug, and the fact that it works for lmer is a bug. Or 2) the error message is a bug, not feature. $\endgroup$ – Curious Jul 24 '13 at 13:13
  • $\begingroup$ See updated post (I had to include some code). I had notice your valid point myself when writing your original response but I originally opted to keep it out so the rationale behind my answer is strictly computation-based. $\endgroup$ – usεr11852 says Reinstate Monic Jul 24 '13 at 13:51
  • 3
    $\begingroup$ @Tomas lmer() doesn't use REML when you ask it do comparisons. IIRC they included some fancy sugar in lmer() so you didn't have to refit the model with ML just to compare fits when you want REML on individual fits to get best estimates of the variance parameters. Look at ?lmer, run the first LME example up to and including the anova(fm1, fm2) call. Look at the log likelihoods reported by anova() and those reported earlier in the printed output. The anova() is getting ML estimates for you. $\endgroup$ – Reinstate Monica - G. Simpson Jul 24 '13 at 15:46
  • $\begingroup$ Good point Gavin, I forget that lmer got both at the same time (it uses PLS so it goes around estimating only one at the time). I forgot to the check what you mentioned. $\endgroup$ – usεr11852 says Reinstate Monic Jul 24 '13 at 16:09
  • 2
    $\begingroup$ @rpierce: The AIC reported within anova() is the one based on ML. The AIC reported just AIC() is the one based on REML. $\endgroup$ – usεr11852 says Reinstate Monic Jul 26 '13 at 1:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.