AIC, ошибка anova: модели не все соответствуют одному и тому же количеству наблюдений, модели не все соответствуют одному и тому же размеру набора данных

У меня есть такие модели:

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)

Сейчас я пытаюсь оценить, должен ли случайный эффект присутствовать в модели. Поэтому я сравниваю модели с использованием AIC или anova и получаю следующую ошибку:

> 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

Как видите, в обоих случаях я использую один и тот же набор данных. Я нашел два лекарства, но я не считаю их удовлетворительными:

1. Добавление method = "ML"в вызов lme () - не уверен, стоит ли менять метод.
2. Использование lmer()вместо. Удивительно, но это работает, несмотря на то, что lmer () использует метод REML. Однако мне не нравится это решение, потому что lmer()оно не показывает p-значения для коэффициентов - lme()вместо этого я предпочитаю использовать более старые .

У вас есть идея, если это ошибка или нет, и как мы можем обойти это?

Быстрый поиск показывает, что это возможно (хотя я должен признать, что я думал, что это не так) и что это не ошибка ... просто еще один случай, когда методы в R скрыты и приводят к тому, что кажется «неожиданным» ', но толпа RTFM говорит: «Это в документации». Во всяком случае ... ваше решение заключается в том, чтобы использовать anovaв lmeкачестве первого аргумента аргументы, а lmмодели - в качестве второго (и третьего, если хотите) аргумента (ов). Если это кажется странным, это потому, что это немного странно. Причина заключается в том, что при вызове anova, то anova.lmeвызывается метод , только если первый аргумент является lmeобъектом. В противном случае он вызывает anova.lm(что, в свою очередь, вызывает anova.lmlist; если вы покопаетесь anova.lm, вы поймете почему). Подробнее о том, как вы хотите звонитьanovaв этом случае откройте справку для anova.lme. Вы увидите, что можете сравнивать другие модели с lmeмоделями, но они должны находиться в позиции, отличной от первого аргумента. По-видимому, также возможно использовать anovaна моделях подгонку, используя glsфункцию, не слишком заботясь о порядке аргументов модели. Но я не знаю достаточно деталей, чтобы определить, является ли это хорошей идеей или нет, или что именно она подразумевает (кажется, это хорошо, но ваш звонок). От этой связи по сравнению lmс по- lmeвидимому, хорошо документированы и привел в качестве метода, так что я бы заблуждаться в этом направлении, будь я вас.

Удачи.

m2m3$ML$REML$y$$k$$kX= 0$method="ML"

Сказав это, давайте посмотрим под капотом:

 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")
     }     
 }

где в вашем случае вы можете увидеть, что:

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

и, очевидно logLik, делает что-то (может быть?) неожиданное ...? нет, не очень, если вы посмотрите на документацию logLik, ?logLikвы увидите это явно указано:

 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’)

что возвращает нас к исходной точке, которую вы должны использовать ML.

Чтобы использовать общее высказывание в CS: «Это не ошибка; это (реальная) особенность!»

РЕДАКТИРОВАТЬ : (просто чтобы ответить на ваш комментарий :) Предположим, вы подходите другой модели, используя lmerэто время:

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

и вы делаете следующее:

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

Что кажется явным несоответствием между ними, но на самом деле это не так, как объяснил Гэвин. Просто чтобы заявить очевидное:

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

Есть веская причина, почему это происходит с точки зрения методологии, я думаю. lmeдействительно пытается понять регрессию LME для вас, в то время как lmerпри сравнении моделей она сразу же возвращается к результатам ML. Я думаю, что нет никаких ошибок в этом вопросе lmeи lmerпросто разные причины для каждого пакета.

См. Также комментарий Гэвина Симпосона о более глубоком объяснении того, что anova()происходило (практически то же самое происходит с AIC)