AIC, Anova-Fehler: Modelle sind nicht alle an die gleiche Anzahl von Beobachtungen angepasst, Modelle wurden nicht alle an die gleiche Größe des Datensatzes angepasst

Ich habe solche Modelle:

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)

Jetzt versuche ich zu beurteilen, ob der zufällige Effekt im Modell vorhanden sein sollte. Ich vergleiche die Modelle mit AIC oder Anova und erhalte folgende Fehlermeldung:

> 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

Wie Sie sehen können, verwende ich in beiden Fällen denselben Datensatz. Ich habe zwei Mittel gefunden, aber ich halte sie nicht für befriedigend:

1. Hinzufügen method = "ML"zum Aufruf von lme () - nicht sicher, ob es sinnvoll ist, die Methode zu ändern.
2. Verwenden Sie lmer()stattdessen. Überraschenderweise funktioniert dies, obwohl lmer () die REML-Methode verwendet. Diese Lösung gefällt mir lmer()jedoch nicht, da die p-Werte für Koeffizienten nicht angezeigt lme()werden. Stattdessen verwende ich gerne ältere .

Haben Sie eine Idee, ob dies ein Fehler ist oder nicht und wie können wir das umgehen?

Eine schnelle Suche zeigt, dass es möglich ist (obwohl ich zugeben muss, dass ich dachte, dass es nicht so ist) und dass es kein Fehler ist ... nur ein weiterer Fall, in dem Methoden in R versteckt sind und zu Dingen führen, die unerwartet erscheinen ', aber die RTFM-Menge sagt: "Es ist in der Dokumentation." Wie auch immer ... Ihre Lösung besteht darin, anovadas lmeals erstes Argument und die lmModelle als zweites (und drittes, wenn Sie möchten) Argument (e) zu verwenden. Wenn dies seltsam erscheint, liegt es daran, dass es etwas seltsam ist. Der Grund ist, dass beim Aufrufen anovadie anova.lmeMethode nur aufgerufen wird, wenn das erste Argument ein lmeObjekt ist. Andernfalls ruft es an anova.lm(was wiederum anruft anova.lmlist; wenn Sie sich vertiefen anova.lm, werden Sie sehen, warum). Einzelheiten dazu, wie Sie anrufen möchtenanovaIn diesem Fall rufen Sie die Hilfe für auf anova.lme. Sie werden sehen, dass Sie andere Modelle mit lmeModellen vergleichen können , diese müssen sich jedoch in einer anderen Position als dem ersten Argument befinden. Anscheinend ist es auch möglich, anovaModelle zu verwenden, die mit der glsFunktion passen, ohne sich zu sehr um die Reihenfolge der Modellargumente zu kümmern. Aber ich weiß nicht genug über die Details, um festzustellen, ob dies eine gute Idee ist oder nicht oder was genau dies impliziert (es scheint wahrscheinlich in Ordnung zu sein, aber Ihr Anruf). Von diesem Link aus scheint der Vergleich lmmit lmegut dokumentiert und als Methode zitiert zu sein, also würde ich mich in diese Richtung irren, wenn ich Sie wäre.

Viel Glück.

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

Nachdem wir das gesagt haben, schauen wir unter die Haube:

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

Wo in Ihrem Fall können Sie das sehen:

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

und macht offensichtlich logLiketwas (vielleicht?) Unerwartetes ...? nein, nicht wirklich, wenn Sie sich die Dokumentation von ansehen logLik, ?logLikwerden Sie sehen, dass es ausdrücklich angegeben ist:

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

Das bringt uns zurück zu unserem ursprünglichen Punkt, den Sie verwenden sollten ML.

Um ein allgemeines Sprichwort in CS zu verwenden: "Es ist kein Fehler, es ist eine (echte) Funktion!"

BEARBEITEN : (Nur um Ihren Kommentar anzusprechen :) Angenommen, Sie passen in lmerdieser Zeit zu einem anderen Modell :

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

und du machst folgendes:

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

Das scheint eine klare Diskrepanz zwischen den beiden zu sein, aber es ist wirklich nicht so, wie Gavin es erklärt hat. Nur um das Offensichtliche zu sagen:

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

Es gibt einen guten Grund, warum dies in Bezug auf die Methodik geschieht, denke ich. lmeversucht, die LME-Regression für Sie zu verstehen, während sie lmerbei Modellvergleichen sofort auf die ML-Ergebnisse zurückgreift. Ich denke, es gibt keine Fehler in dieser Angelegenheit lmeund lmernur unterschiedliche Gründe für jedes Paket.

Siehe auch Gavin Simposons Kommentar zu einer aufschlussreicheren Erklärung dessen, was vor sich ging anova()(das gleiche passiert praktisch mit AIC)