AIC, błąd anova: nie wszystkie modele są dopasowane do tej samej liczby obserwacji, nie wszystkie modele są dopasowane do tego samego rozmiaru zbioru danych

Mam takie modele:

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)

Teraz próbuję ocenić, czy efekt losowy powinien być obecny w modelu. Porównuję więc modele za pomocą AIC lub anova i otrzymuję następujący błąd:

> 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

Jak widać, w obu przypadkach korzystam z tego samego zestawu danych. Znalazłem dwa lekarstwa, ale nie uważam ich za zadowalające:

1. Dodanie method = "ML"do wywołania lme ​​() - nie jestem pewien, czy warto zmienić metodę.
2. Używanie lmer()zamiast. O dziwo, działa to pomimo faktu, że lmer () używa metody REML. Jednak nie podoba mi się to rozwiązanie, ponieważ lmer()nie pokazuje wartości p dla współczynników - lme()zamiast tego wolę używać starszych .

Czy masz pojęcie, czy to błąd, czy nie i jak możemy to obejść?

Szybkie wyszukiwanie pokazuje, że jest to możliwe (chociaż muszę przyznać, że myślałem, że tak nie jest) i że nie jest to błąd ... to po prostu kolejny przypadek, w którym metody w R są ukryte i powodują rzeczy, które wydają się „nieoczekiwane” ”, ale tłum RTFM mówi:„ Jest w dokumentacji ”. W każdym razie ... rozwiązanie jest zrobić anovaz lmejako pierwszy argument i lmmodeli jako drugi (i trzeci, jeśli lubisz) argumentu (ów). Jeśli wydaje się to dziwne, to dlatego, że jest trochę dziwne. Powodem jest to, że podczas rozmowy anovaThe anova.lmemetoda jest wywoływana tylko jeśli pierwszy argument jest lmeobiektem. W przeciwnym razie wywołuje anova.lm(co z kolei wywołuje anova.lmlist; jeśli się zagłębisz anova.lm, zobaczysz dlaczego). Szczegółowe informacje o tym, jak chcesz dzwonićanovaw takim przypadku poproś o pomoc anova.lme. Zobaczysz, że możesz porównać inne modele z lmemodelami, ale muszą one znajdować się w pozycji innej niż pierwszy argument. Najwyraźniej możliwe jest również użycie anovaw modelach dopasowanych za pomocą glsfunkcji bez zbytniego dbania o kolejność argumentów modelu. Ale nie znam wystarczających szczegółów, aby ustalić, czy to dobry pomysł, czy nie, lub co dokładnie implikuje (wydaje się prawdopodobnie w porządku, ale twój telefon). Od tego linku lmdo porównania lmewydaje się być dobrze udokumentowane i cytowane jako metoda, więc popełniłbym błąd w tym kierunku, gdybym był tobą.

Powodzenia.

To jest zdecydowanie szczególne. Pierwsza myśl: gdy porównujemy modele, w których modele mają różne struktury efektów stałych ( m2i m3na przykład), najlepiej dla nas$ML$jak REMLsię zmieni$y$. (Pomnoży to$k$, gdzie $kX= 0$) To ciekawe, że działa, method="ML"co sprawia, że ​​uważam, że to nie błąd. Wygląda na to, że wymusza „dobrą praktykę”.

Powiedziawszy to, spójrzmy pod maską:

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

gdzie w twoim przypadku możesz zobaczyć, że:

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

i oczywiście logLikrobi coś (może?) nieoczekiwanego ...? nie, nie bardzo, jeśli spojrzysz na dokumentację logLik, ?logLikzobaczysz, że jest to wyraźnie określone:

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

co przywraca nas do pierwotnego punktu, którego powinieneś używać ML.

Aby użyć wspólnego powiedzenia w CS: „To nie jest błąd; to (prawdziwa) funkcja!”

EDYCJA : (tylko w celu skomentowania komentarza :) Załóżmy, że pasujesz do innego modelu, używając lmertego czasu:

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

i wykonujesz następujące czynności:

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

Co wydaje się wyraźną rozbieżnością między nimi, ale tak naprawdę nie jest tak, jak wyjaśnił Gavin. Wystarczy powiedzieć oczywiste:

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

Myślę, że istnieje dobry powód, dla którego dzieje się tak pod względem metodologii. lmestara się zrozumieć dla ciebie regresję LME, podczas lmergdy podczas porównywania modeli natychmiast wraca do wyników ML. Myślę, że nie ma żadnych błędów w tej sprawie lmei lmerpo prostu inne uzasadnienia dla każdego pakietu.

Zobacz także komentarz Gavina Simposona na bardziej wnikliwe wyjaśnienie tego, co się działo anova()(To samo dzieje się praktycznie z AIC)