REML oder ML, um zwei Modelle mit gemischten Effekten mit unterschiedlichen festen Effekten zu vergleichen, aber mit demselben zufälligen Effekt?

18

Hintergrund: Hinweis: Mein Datensatz und R-Code sind unter dem Text enthalten

Ich möchte AIC verwenden, um zwei Modelle mit gemischten Effekten zu vergleichen, die mit dem lme4-Paket in R erstellt wurden. Jedes Modell hat einen festen und einen zufälligen Effekt. Der festgelegte Effekt unterscheidet sich zwischen den Modellen, der Zufallseffekt bleibt jedoch zwischen den Modellen gleich. Ich habe festgestellt, dass, wenn ich REML = T verwende, Modell2 die niedrigere AIC-Bewertung hat, aber wenn ich REML = F verwende, hat Modell1 die niedrigere AIC-Bewertung.

Unterstützung für die Verwendung von ML:

Zuur et al. (2009; SEITE 122) schlagen vor: "Um Modelle mit verschachtelten festen Effekten (aber mit derselben zufälligen Struktur) zu vergleichen, muss die ML-Schätzung und nicht REML verwendet werden." Dies bedeutet für mich, dass ich ML verwenden sollte, da meine zufälligen Effekte in beiden Modellen gleich sind, meine festen Effekte sich jedoch unterscheiden. [Zuur et al. 2009. Mixed Effect Models und Extensions in Ecology mit R. Springer.]

Unterstützung für die Verwendung von REML:

Ich stelle jedoch fest, dass die mit den Zufallseffekten verbundene Restvarianz bei der Verwendung von ML zwischen den beiden Modellen unterschiedlich ist (Modell1 = 136,3; Modell2 = 112,9). Bei Verwendung von REML ist sie jedoch bei den Modellen identisch (Modell1 = Modell2 =) 151,5). Dies impliziert für mich, dass ich stattdessen REML verwenden sollte, damit die zufällige Restvarianz zwischen Modellen mit derselben Zufallsvariablen gleich bleibt.

Frage:

Ist es nicht sinnvoller, REML als ML für den Vergleich von Modellen zu verwenden, bei denen sich die festen Effekte ändern und die zufälligen Effekte gleich bleiben? Wenn nicht, können Sie mir erklären, warum oder auf andere Literatur verweisen, die mehr erklärt?

# Model2 "wins" if REML=T:
REMLmodel1 = lmer(Response ~ Fixed1 + (1|Random1),data,REML = T)
REMLmodel2 = lmer(Response ~ Fixed2 + (1|Random1),data,REML = T)
AIC(REMLmodel1,REMLmodel2)
summary(REMLmodel1)
summary(REMLmodel2)

# Model1 "wins" if REML=F:
MLmodel1 = lmer(Response ~ Fixed1 + (1|Random1),data,REML = F)
MLmodel2 = lmer(Response ~ Fixed2 + (1|Random1),data,REML = F)
AIC(MLmodel1,MLmodel2)
summary(MLmodel1)
summary(MLmodel2)

Datensatz:

Response    Fixed1  Fixed2  Random1
5.20    A   A   1
32.50   A   A   1
6.57    A   A   2
24.77   A   B   3
41.69   A   B   3
34.29   A   B   4
1.80    A   B   4
10.00   A   B   5
15.56   A   B   5
4.44    A   C   6
21.65   A   C   6
9.20    A   C   7
4.11    A   C   7
12.52   B   D   8
0.25    B   D   8
27.34   B   D   9
11.54   B   E   10
0.86    B   E   10
0.68    B   E   11
4.00    B   E   11
Es zeigt
quelle
2
Faraway's (2006) Erweiterung des linearen Modells um R (S. 156): "Der Grund dafür ist, dass REML die Zufallseffekte durch Berücksichtigung linearer Kombinationen der Daten abschätzt, die die festen Effekte entfernen. Wenn diese festen Effekte geändert werden, werden die Wahrscheinlichkeiten von zwei Modelle werden nicht direkt vergleichbar sein. "
Jvh_ch
Obwohl AIC nach meinem besten Wissen auf Wahrscheinlichkeit basiert, wurde es zum Zweck der Vorhersage entwickelt. Wie würde man ein gemischtes Modell für die Vorhersage genau anwenden?
AdamO
@AdamO, könnten Sie präziser sein? Ein angepasstes gemischtes Modell kann für die Vorhersage verwendet werden, entweder auf Populationsebene (Vorhersage der Antworten für eine nicht spezifizierte / unbekannte Einheit durch Setzen der bedingten Modi / BLUPs auf Null) oder auf individueller Ebene (Bedingungsvorhersage für die Schätzungen der bedingten Modi / BLUPs) ). Wenn Sie genauer sein können, ist dies möglicherweise eine gute neue Frage zum Lebenslauf.
Ben Bolker
Mir war nur unklar, wie Sie dieses Modell anwenden wollten. Nichts in dem Problem deutete darauf hin, welche Art von Vorhersage, wenn überhaupt, durchgeführt wurde oder ob sie notwendig war und wenn ja, zu welchem ​​Zweck.
AdamO

Antworten:

22

Zuur et al. Und Faraway (aus @ janhoves Kommentar oben) haben recht; Die Verwendung von Wahrscheinlichkeitsmethoden (einschließlich AIC) zum Vergleich von zwei Modellen mit unterschiedlichen festen Effekten, die durch REML angepasst werden, führt im Allgemeinen zu Unsinn.

Ben Bolker
quelle
4
Danke @janhove, AdamO und Ben Bolker. Ich fand auch diesen Link von Aaron hilfreich bei der Beantwortung dieser Frage. Darin heißt es: "Die REML-Wahrscheinlichkeit hängt davon ab, welche festen Effekte im Modell enthalten sind, und ist daher nicht vergleichbar, wenn sich die festen Effekte ändern. Da REML im Allgemeinen bessere Schätzungen für die zufälligen Effekte liefert, sollten die üblichen Empfehlungen eingehalten werden Ihr bestes Modell mit REML für Ihre endgültige Schlussfolgerung und Berichterstattung. "
Es Zahlen
11

Ich werde ein Beispiel geben, um zu veranschaulichen, warum die REML-Wahrscheinlichkeit nicht für Dinge wie AIC-Vergleiche verwendet werden kann. Stellen Sie sich vor, wir hätten ein normales Mischeffektmodell. Lassen bezeichnen die Designmatrix und davon ausgehen , dass diese Matrix vollen Rang hat. Wir können eine Reparametrisierung des Mittelwertraums finden, die durch die Matrix . Die beiden Matrizen überspannen denselben linearen Unterraum von . Somit können die Spalten von als lineare Kombinationen der Spalten von . Daher können wir eine quadratische Matrix , so dassXX~RnX~XB

X~=XB .

Darüber hinaus hat den vollen Rang (dies kann durch die Annahme bewiesen werden, dass dies nicht der Fall ist; dann wäre auch kein Widerspruch). Dies bedeutet, dass invertierbar ist.BXB

Wenn wir mit der zweiten Parametrisierung des Mittelwertraums beginnen und eine Kovarianzmatrix sein lassen, betrachten wir das REML-Kriterium, das wir maximieren sollten (ich lasse eine Konstante weg).V

|V|1/2|X~V1X~|1/2exp((yX~β~)V1(yX~β~)/2) ,

über den Parametersatz, wobei . Anhand der Tatsache, dass , können wir erkennen, dass dies umgeschrieben werden kann alsβ=(X~V1X~)1yX=X~B

|B||V|1/2||XV1X|1/2|exp((yXβ¯)V1(yXβ¯)/2) ,

wobei . Dies ist die REML-Wahrscheinlichkeit für die anderen Parametrisierungszeiten der Determinante von.| B|β¯=(XV1X)1y|B|

Wir haben daher ein Beispiel für zwei verschiedene Parametrisierungen desselben Modells, die unterschiedliche Wahrscheinlichkeitswerte ergeben, unter der Annahme, dass (eine solche Matrix kann leicht gefunden werden). Der gleiche Parameterwert maximiert in beiden Fällen das Kriterium, aber der Wert der Wahrscheinlichkeit ist unterschiedlich. Dies zeigt, dass der Wahrscheinlichkeitswert ein willkürliches Element enthält, und veranschaulicht daher, warum der Wert der Wahrscheinlichkeit nicht für den Vergleich zwischen Modellen mit unterschiedlichen festen Effekten verwendet werden kann: Sie könnten die Ergebnisse ändern, indem Sie einfach die Mittelwertraumparametrierung in ändern eines der Modelle.|B|1

Dies ist ein Beispiel dafür, warum REML beim Vergleich von Modellen mit unterschiedlichen festen Effekten nicht verwendet werden sollte. REML schätzt jedoch die Zufallseffektparameter häufig besser. Daher wird manchmal empfohlen, ML für Vergleiche und REML für die Schätzung eines einzelnen (möglicherweise endgültigen) Modells zu verwenden.

swmo
quelle