Das nicht so seltene Auftreten bei komplexen maximal gemischten Modellen (Schätzung aller möglichen zufälligen Effekte für bestimmte Daten und Modelle) ist eine perfekte (+1 oder -1) oder nahezu perfekte Korrelation zwischen einigen zufälligen Effekten. Betrachten wir zum Zweck der Diskussion das folgende Modell und die folgende Modellzusammenfassung
Model: Y ~ X*Cond + (X*Cond|subj)
# Y = logit variable
# X = continuous variable
# Condition = values A and B, dummy coded; the design is repeated
# so all participants go through both Conditions
# subject = random effects for different subjects
Random effects:
Groups Name Variance Std.Dev. Corr
subject (Intercept) 0.85052 0.9222
X 0.08427 0.2903 -1.00
CondB 0.54367 0.7373 -0.37 0.37
X:CondB 0.14812 0.3849 0.26 -0.26 -0.56
Number of obs: 39401, groups: subject, 219
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.49686 0.06909 36.14 < 2e-16 ***
X -1.03854 0.03812 -27.24 < 2e-16 ***
CondB -0.19707 0.06382 -3.09 0.00202 **
X:CondB 0.22809 0.05356 4.26 2.06e-05 ***
Der vermeintliche Grund für diese perfekten Korrelationen ist, dass wir ein Modell erstellt haben, das für die uns vorliegenden Daten zu komplex ist. Der in diesen Situationen häufig gegebene Rat ist (z. B. Matuschek et al., 2017; Papier ), die überparametrisierten Koeffizienten auf 0 zu setzen, da solche entarteten Modelle dazu neigen, die Leistung zu verringern. Wenn wir in einem reduzierten Modell eine deutliche Änderung der festen Effekte beobachten, sollten wir diese akzeptieren. Wenn sich nichts ändert, ist es kein Problem, das Original zu akzeptieren.
Nehmen wir jedoch an, dass wir nicht nur an festen Effekten interessiert sind, die für RE gesteuert werden (zufällige Effekte), sondern auch an der RE-Struktur. In dem gegebenen Fall wäre es theoretisch vernünftig anzunehmen, dass Intercept
und die Steigung X
eine negative Korrelation ungleich Null aufweisen. Es folgen mehrere Fragen:
Was tun in solchen Situationen? Sollten wir die perfekte Korrelation angeben und sagen, dass unsere Daten nicht "gut genug" sind, um die "echte" Korrelation abzuschätzen? Oder sollten wir das 0-Korrelationsmodell angeben? Oder sollten wir vielleicht versuchen, eine andere Korrelation auf 0 zu setzen, in der Hoffnung, dass die "wichtige" nicht mehr perfekt ist? Ich glaube nicht, dass es hier 100% richtige Antworten gibt, ich würde gerne Ihre Meinung hören.
Wie schreibe ich den Code, der die Korrelation von 2 spezifischen Zufallseffekten auf 0 festlegt, ohne die Korrelationen zwischen anderen Parametern zu beeinflussen?
quelle
blme
,MCMCglmm
,rstanarm
,brms
...)Antworten:
Singuläre Kovarianzmatrizen mit zufälligen Effekten
Das Erhalten einer zufälligen Effektkorrelationsschätzung von +1 oder -1 bedeutet, dass der Optimierungsalgorithmus "eine Grenze" erreicht: Korrelationen können nicht höher als +1 oder niedriger als -1 sein. Selbst wenn keine expliziten Konvergenzfehler oder Warnungen vorliegen, weist dies möglicherweise auf einige Probleme mit der Konvergenz hin, da wir nicht erwarten, dass echte Korrelationen an der Grenze liegen. Wie Sie sagten, bedeutet dies normalerweise, dass nicht genügend Daten vorhanden sind, um alle Parameter zuverlässig zu schätzen. Matuschek et al. 2017 sagen, dass in dieser Situation die Leistung beeinträchtigt werden kann.
Eine andere Möglichkeit, eine Grenze zu erreichen, besteht darin, eine Varianzschätzung von 0 zu erhalten: Warum erhalte ich in meinem gemischten Modell trotz einiger Abweichungen in den Daten eine Nullvarianz eines zufälligen Effekts?
In beiden Situationen kann eine entartete Kovarianzmatrix mit zufälligen Effekten erhalten werden (in Ihrem Beispiel beträgt die Ausgabekovarianzmatrix ). Eine Nullvarianz oder eine perfekte Korrelation bedeutet, dass die Kovarianzmatrix nicht den vollen Rang hat und [mindestens] einer ihrer Eigenwerte Null ist. Diese Beobachtung legt sofort nahe, dass es andere , komplexere Wege gibt, um eine entartete Kovarianzmatrix zu erhalten: Man kann eine Kovarianzmatrix ohne Nullen oder perfekte Korrelationen haben, aber dennoch rangdefizient (Singular). Bates et al. 2015 sparsame gemischte Modelle4 × 44 × 4 4 × 4 (unveröffentlichter Preprint) empfehlen die Verwendung der Hauptkomponentenanalyse (PCA), um zu überprüfen, ob die erhaltene Kovarianzmatrix singulär ist. Wenn dies der Fall ist, schlagen sie vor, diese Situation genauso zu behandeln wie die oben genannten singulären Situationen.
Was tun?
Wenn nicht genügend Daten vorhanden sind, um alle Parameter eines Modells zuverlässig abzuschätzen, sollten wir eine Vereinfachung des Modells in Betracht ziehen. Anhand Ihres Beispielmodells
X*Cond + (X*Cond|subj)
gibt es verschiedene Möglichkeiten, es zu vereinfachen:Entfernen Sie einen der zufälligen Effekte, normalerweise die Korrelation höchster Ordnung:
Entfernen Sie alle Korrelationsparameter:
Update: Wie @Henrik bemerkt, entfernt die
||
Syntax Korrelationen nur, wenn alle Variablen links davon numerisch sind. Wenn kategoriale Variablen (wieCond
) beteiligt sind, sollte man lieber sein praktischesafex
Paket (oder umständliche manuelle Problemumgehungen) verwenden. Siehe seine Antwort für weitere Details.Entfernen Sie einige der Korrelationsparameter, indem Sie den Begriff in mehrere aufteilen, z.
lme4
, um dies zu erreichen. Siehe @ BenBolker Antwort auf SO für eine Demonstration, wie dies über einige Smart - Hacking zu erreichen.Im Gegensatz zu Ihren Aussagen glaube ich nicht, dass Matuschek et al. 2017 empfehlen speziell # 4. Der Kern von Matuschek et al. 2017 und Bates et al. 2015 scheint es so zu sein, dass man mit dem Maximalmodell a la Barr et al. 2013 und verringert dann die Komplexität, bis die Kovarianzmatrix den vollen Rang hat. (Darüber hinaus empfehlen sie häufig, die Komplexität noch weiter zu reduzieren, um die Leistung zu erhöhen.) Update: Im Gegensatz dazu haben Barr et al. empfehlen, die Komplexität NUR zu reduzieren, wenn das Modell nicht konvergiert; Sie sind bereit, singuläre Kovarianzmatrizen zu tolerieren. Siehe @ Henriks Antwort.
Wenn man Bates / Matuschek zustimmt, dann denke ich, ist es in Ordnung, verschiedene Möglichkeiten auszuprobieren, um die Komplexität zu verringern, um denjenigen zu finden, der die Arbeit erledigt, während er "den geringsten Schaden" anrichtet. In meiner obigen Liste enthält die ursprüngliche Kovarianzmatrix 10 Parameter. # 1 hat 6 Parameter, # 2 hat 4 Parameter, # 3 hat 7 Parameter. Welches Modell die perfekten Korrelationen beseitigt, kann man nicht sagen, ohne sie anzupassen.
Aber was ist, wenn Sie an diesem Parameter interessiert sind?
In der obigen Diskussion wird die Kovarianzmatrix mit zufälligen Effekten als Störparameter behandelt. Sie werfen eine interessante Frage auf, was zu tun ist, wenn Sie speziell an einem Korrelationsparameter interessiert sind, den Sie "aufgeben" müssen, um eine aussagekräftige Lösung mit vollem Rang zu erhalten.
Beachten Sie, dass das Festlegen des Korrelationsparameters auf Null nicht unbedingt zu BLUPs (
ranef
) führt, die nicht korreliert sind. Tatsächlich sind sie möglicherweise gar nicht so stark betroffen (siehe @ Placidias Antwort für eine Demonstration ). Eine Möglichkeit wäre also, die Korrelationen von BLUPs zu betrachten und dies zu melden.Eine andere, vielleicht weniger attraktive Option wäre, die Behandlung
subject
als festen Effekt zu verwendenY~X*cond*subj
, die Schätzungen für jedes Subjekt zu erhalten und die Korrelation zwischen ihnen zu berechnen. Dies entspricht dem Ausführen separaterY~X*cond
Regressionen für jedes Subjekt und dem Abrufen der Korrelationsschätzungen.Siehe auch den Abschnitt über einzelne Modelle in Ben Bolkers FAQ zu gemischten Modellen:
quelle
(Machine||Worker)
lmer
Schätzungen eine Varianz mehr als für(Machine|Worker)
. Waslmer
also||
mit Faktoren zu tun hat, kann nicht beschrieben werden durch "Dies beseitigt nur Korrelationen zwischen Faktoren, aber nicht zwischen Ebenen eines kategorialen Faktors." Es verändert die Struktur der Zufallseffekte auf etwas seltsame Weise (es erweitert sich(Machine||Worker)
auf(1|Worker) + (0+Machine|Worker)
, daher die zusätzliche Varianz). Fühlen Sie sich frei, meine Bearbeitung zu ändern. Mein Hauptpunkt ist, dass in diesen Aussagen die Unterscheidung zwischen numerischen und kategorialen Kovariaten klargestellt werden muss.machines2 <- subset(Machines, Machine %in% c("A", "B")); summary(lmer(score ~ Machine + (Machine || Worker), data=machines2))
. Es funktioniert im Allgemeinen nicht mit Faktoren aufgrund dieser Erweiterung und der Art und Weise, wieR
mit Faktoren umgegangen wirdmodel.matrix
.ranef
Werten zuzuwenden, um die Korrelation zwischen zufälligen Effekten zu untersuchen. Ich bin nicht zu tief in dieses Thema vertieft, aber ich weiß, dass es normalerweise nicht empfohlen wird, mit extrahierten Werten von zu arbeitenranef
, sondern mit geschätzten Korrelationen und Varianzen. Was ist deine Meinung dazu? Außerdem weiß ich nicht, wie man dem Prüfer erklären würde, dass die Korrelation im Modell nicht postuliert wurde, aber wir berechnen trotzdem die Korrelation der extrahierten Werte. Das macht keinen SinnIch stimme mit allem überein, was in Amöbas Antwort gesagt wurde, die eine großartige Zusammenfassung der aktuellen Diskussion zu diesem Thema bietet. Ich werde versuchen, ein paar zusätzliche Punkte hinzuzufügen und ansonsten auf das Handout meines letzten gemischten Modellkurses verweisen, in dem auch diese Punkte zusammengefasst sind.
Das Unterdrücken der Korrelationsparameter (Optionen 2 und 3 in Amöbens Antwort) über
||
funktioniert nur für numerische Kovariaten inlmer
und nicht für Faktoren. Dies wird mit Code von Reinhold Kliegl ausführlich besprochen .Mein
afex
Paket bietet jedoch die Funktionalität, um die Korrelation auch zwischen Faktoren zu unterdrücken, wenn Argumenteexpand_re = TRUE
im Aufruf vonmixed()
(siehe auch Funktionlmer_alt()
). Dies geschieht im Wesentlichen durch Implementierung des von Reinhold Kliegl diskutierten Ansatzes (dh Umwandlung der Faktoren in numerische Kovariaten und Angabe der Zufallseffektstruktur auf diesen).Ein einfaches Beispiel:
Für diejenigen, die es nicht wissen
afex
, besteht die Hauptfunktionalität für gemischte Modelle darin, p-Werte für die festen Effekte bereitzustellen, z.Dale Barr von Barr et al. (2013) empfehlen vorsichtiger, die Struktur der Zufallseffekte zu reduzieren, als in der Antwort von Amöbe dargestellt. In einem kürzlichen Twitter-Austausch schrieb er:
Daher ist Vorsicht geboten.
Als einer der Gutachter kann ich auch einen Einblick geben, warum wir, Bates et al. (2015) blieb unveröffentlicht. Ich und die beiden anderen Gutachter (die unterschrieben haben, aber hier nicht genannt werden) hatten einige Kritik am PCA-Ansatz (scheint nicht prinzipiell zu sein und es gibt keine Beweise dafür, dass er in Bezug auf die Macht überlegen ist). Darüber hinaus glaube ich, dass alle drei kritisiert haben, dass sich das Papier nicht auf die Frage der Spezifizierung der Zufallseffektstruktur konzentriert, sondern auch versucht, GAMMs einzuführen. So verwandelte sich das Papier von Bates et al. (2015) in das Papier von Matuschek et al. (2017) , das sich mit dem Problem der Zufallseffektstruktur mit Simulationen befasst, und das von Baayen et al. (2017) Papier zur Einführung von GAMMs.
Meine vollständige Übersicht über Bates et al. Entwurf finden Sie hier . IIRC, die anderen Bewertungen hatten ähnliche Hauptpunkte.
quelle
lmer_alt
funktioniert im Grunde genau wielmer
(oder sogarglmer
) mit dem einzigen Unterschied, dass es die||
Syntax erlaubt . Ich bin mir also nicht sicher, warum Sie dies unbedingt vermeiden möchtenafex
. Es sollte sogar ohne Anbringen funktionieren (dhafex::lmer_alt(...)
).cbind
zu den Daten. Dann wird der Zufallseffektterm in der Formel durch einen neuen ersetzt, in dem jede der neu erstellten Spalten mit einem + verknüpft ist. Siehe Zeilen 690 bis 730 in github.com/singmann/afex/blob/master/R/mixed.R||
ist dies ein wirklich wichtiger Punkt, danke, dass Sie ihn angesprochen und mir erklärt haben (ich habe meine Antwort bearbeitet, um sie wiederzugeben). Ich mag diese Funktionalität vonlmer_alt
inafex
. Der Vollständigkeit halberlmer
möchte ich hier nur erwähnen, dass man z(1+dummy(Machine,'B')+dummy(Machine,'C') || Worker)
. B. angeben kann , um die gleiche Ausgabe mit Vanilla Call ohne zusätzliche Vorverarbeitung zu erhalten . Dies wird eindeutig sehr umständlich, wenn kategoriale Variablen viele Ebenen haben.dummy()
nur mit den Standardbehandlungskontrasten funktioniert und nicht, wenn die Zufallseffekte Summen-Null-Kontraste verwenden (die verwendet werden sollten, wenn das Modell Interaktionen aufweist). Sie können dies beispielsweise sehen, wenn Sie die Varianzkomponenten im obigen Beispiel für denlmer_alt
Anruf mit demmixed
Anruf vergleichen.Auch ich hatte dieses Problem bei der Verwendung der Maximum-Likelihood-Schätzung - nur ich verwende den Goldstein IGLS-Algorithmus, wie er durch die MLwiN-Software implementiert wurde, und nicht LME4 in R. In jedem Fall hat sich das Problem jedoch gelöst, als ich mit derselben zur MCMC-Schätzung gewechselt bin Software. Ich hatte sogar eine Korrelation von mehr als 3, die sich auflöste, als ich die Schätzung änderte. Unter Verwendung von IGLS wird die Korrelation nach der Schätzung als Kovarianz geteilt durch das Produkt der Quadratwurzel des Produkts der zugehörigen Varianzen berechnet - und dies berücksichtigt nicht die Unsicherheit in jeder der konstituierenden Schätzungen.
Die IGLS-Software "weiß" nicht, dass die Kovarianz eine Korrelation impliziert, und berechnet lediglich Schätzungen einer konstanten, linearen, quadratischen usw. Varianzfunktion. Im Gegensatz dazu basiert der MCMC-Ansatz auf der Annahme von Stichproben aus einer multivariaten Normalverteilung, die Varianzen und Kovarianzen mit guten Eigenschaften und vollständiger Fehlerausbreitung entspricht, so dass die Unsicherheit bei der Schätzung der Kovarianzen bei der Schätzung der Varianzen berücksichtigt wird und umgekehrt.
MLwin ist die MCMC-Schätzkette mit IGLS-Schätzungen, und die nicht negative Kovarianzmatrix mit definitiver Varianz muss möglicherweise geändert werden, indem die Kovarianz zu Beginn vor Beginn der Stichprobe auf Null geändert wird.
Ein Beispiel finden Sie unter
Entwicklung von Mehrebenenmodellen zur Analyse von Kontextualität, Heterogenität und Veränderung unter Verwendung von MLwiN 3, Band 1 (aktualisiert September 2017); Band 2 ist auch auf RGate
https://www.researchgate.net/publication/320197425_Vol1Training_manualRevisedSept2017
Anhang zu Kapitel 10
quelle