Was tun mit einer Korrelation mit zufälligen Effekten, die gleich 1 oder -1 ist?

9

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 Interceptund die Steigung Xeine negative Korrelation ungleich Null aufweisen. Es folgen mehrere Fragen:

  1. 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.

  2. Wie schreibe ich den Code, der die Korrelation von 2 spezifischen Zufallseffekten auf 0 festlegt, ohne die Korrelationen zwischen anderen Parametern zu beeinflussen?

User33268
quelle
Mit Paket nlme können Sie die Varianz-Kovarianz-Matrix zufälliger Effekte genau steuern. Ich habe das selbst nie wirklich gebraucht, aber ich würde Mixed-Effects-Modelle in S und S-PLUS (Pinheiro und Bates, 2000) noch einmal lesen , wenn ich das tun würde.
Roland
3
Eine radikale Alternative ist , um das Modell regularisieren, dh ein Bayes - Modell mit etwas informativ priors auf den zufälligen Effekten Strukturen passen (zB über blme, MCMCglmm, rstanarm, brms...)
Ben Bolker
@ BenBolker Ben. Ich bin mir nicht sicher, ob dies eine radikale Idee ist, da das Anpassen eines unregelmäßigen Modells der radikalere Weg sein könnte, um ein Modell
anzupassen
Ich danke euch allen für die tollen Antworten ... Leider war ich ein paar Tage offline, aber ich bin zurück.
User33268

Antworten:

13

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×44×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:

  1. Entfernen Sie einen der zufälligen Effekte, normalerweise die Korrelation höchster Ordnung:

    X*Cond + (X+Cond|subj)
  2. Entfernen Sie alle Korrelationsparameter:

    X*Cond + (X*Cond||subj)

    Update: Wie @Henrik bemerkt, entfernt die ||Syntax Korrelationen nur, wenn alle Variablen links davon numerisch sind. Wenn kategoriale Variablen (wie Cond) beteiligt sind, sollte man lieber sein praktisches afexPaket (oder umständliche manuelle Problemumgehungen) verwenden. Siehe seine Antwort für weitere Details.

  3. Entfernen Sie einige der Korrelationsparameter, indem Sie den Begriff in mehrere aufteilen, z.

    X*Cond + (X+Cond|subj) + (0+X:Cond|subj)
  4. Beschränken Sie die Kovarianzmatrix auf eine bestimmte Weise, z. B. indem Sie eine bestimmte Korrelation (diejenige, die die Grenze erreicht hat) auf Null setzen, wie Sie vorschlagen. Es gibt keinen eingebauten Weg 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 subjectals festen Effekt zu verwenden Y~X*cond*subj, die Schätzungen für jedes Subjekt zu erhalten und die Korrelation zwischen ihnen zu berechnen. Dies entspricht dem Ausführen separater Y~X*condRegressionen für jedes Subjekt und dem Abrufen der Korrelationsschätzungen.


Siehe auch den Abschnitt über einzelne Modelle in Ben Bolkers FAQ zu gemischten Modellen:

Es ist sehr häufig, dass überangepasste gemischte Modelle zu singulären Anpassungen führen. Technisch bedeutet Singularität, dass einige der -Parameter (Varianz-Kovarianz-Cholesky-Zerlegung), die diagonalen Elementen des Cholesky-Faktors entsprechen, genau Null sind, was der Rand des realisierbaren Raums ist, oder äquivalent, dass die Varianz-Kovarianz-Matrix etwas Null hat Eigenwerte (dh eher positiv semidefinit als positiv definit) oder (fast gleichwertig), dass einige der Varianzen als Null oder einige der Korrelationen als +/- 1 geschätzt werden.θ

Amöbe
quelle
1
Was mein Beispiel zeigt, ist, dass für (Machine||Worker) lmerSchätzungen eine Varianz mehr als für (Machine|Worker). Was lmeralso ||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.
Henrik
1
Nein, funktioniert auch nicht mit binären Variablen, überzeugen Sie sich selbst : 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, wie Rmit Faktoren umgegangen wird model.matrix.
Henrik
@amoeba: Ich denke, Sie haben einen interessanten Punkt hervorgehoben, indem Sie vorgeschlagen haben, sich den ranefWerten 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 arbeiten ranef, 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 Sinn
User33268
1
@RockyRaccoon Ja, ich denke, es ist besser, den geschätzten Korrelationsparameter zu verwenden / zu melden, aber hier sprechen wir über die Situation, in der wir ihn wohl nicht schätzen können, weil er gegen 1 konvergiert. Das würde ich in einem Artikel schreiben: "Das vollständige Modell konvergierte Um mit corr = 1 zu lösen, verwendeten wir gemäß den Empfehlungen in [Zitierungen] ein reduziertes Modell [Details]. Die Korrelation zwischen BLUPs mit zufälligen Effekten in diesem Modell betrug 0,9. " Wenn Sie die Korrelation nicht einbeziehen, beschränken Sie das Modell nicht darauf, sie als unkorreliert zu behandeln! Sie modellieren diese Korrelation einfach nicht explizit.
Amöbe
Ich habe noch eine Frage: Bedeuten Varianzen nahe Null und perfekte und nahezu perfekte Korrelationen von Zufallseffekten etwas über den tatsächlichen Wert der Parameter? Bedeuten -1-Korrelationen beispielsweise, dass die reale Korrelation mindestens negativ und / oder mindestens ungleich Null ist? Konkreter: Wenn wir versuchen, die Korrelation zu schätzen, die in Wirklichkeit 0 ist, ist es dann möglich, dass wir eine Schätzung von -1 erhalten?
User33268
9

Ich 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 in lmerund nicht für Faktoren. Dies wird mit Code von Reinhold Kliegl ausführlich besprochen .

Mein afexPaket bietet jedoch die Funktionalität, um die Korrelation auch zwischen Faktoren zu unterdrücken, wenn Argumente expand_re = TRUEim Aufruf von mixed()(siehe auch Funktion lmer_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:

library("afex")
data("Machines", package = "MEMSS") # same data as in Kliegl code

# with correlation:
summary(lmer(score ~ Machine + (Machine  | Worker), data=Machines))
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr       
#  Worker   (Intercept) 16.6405  4.0793              
#           MachineB    34.5467  5.8776    0.48      
#           MachineC    13.6150  3.6899   -0.37  0.30
#  Residual              0.9246  0.9616              
# Number of obs: 54, groups:  Worker, 6

## crazy results:
summary(lmer(score ~ Machine + (Machine  || Worker), data=Machines))
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr     
#  Worker   (Intercept)  0.2576  0.5076            
#  Worker.1 MachineA    16.3829  4.0476            
#           MachineB    74.1381  8.6103   0.80     
#           MachineC    19.0099  4.3600   0.62 0.77
#  Residual              0.9246  0.9616            
# Number of obs: 54, groups:  Worker, 6

## as expected:
summary(lmer_alt(score ~ Machine + (Machine  || Worker), data=Machines))
# Random effects:
#  Groups   Name         Variance Std.Dev.
#  Worker   (Intercept)  16.600   4.0743  
#  Worker.1 re1.MachineB 34.684   5.8894  
#  Worker.2 re1.MachineC 13.301   3.6471  
#  Residual               0.926   0.9623  
# Number of obs: 54, groups:  Worker, 6

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.

(m1 <- mixed(score ~ Machine + (Machine  || Worker), data=Machines, expand_re = TRUE))
# Mixed Model Anova Table (Type 3 tests, KR-method)
# 
# Model: score ~ Machine + (Machine || Worker)
# Data: Machines
#    Effect      df        F p.value
# 1 Machine 2, 5.98 20.96 **    .002
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

summary(m1)  
# [...]
# Random effects:
#  Groups   Name         Variance Std.Dev.
#  Worker   (Intercept)  27.4947  5.2435  
#  Worker.1 re1.Machine1  6.6794  2.5845  
#  Worker.2 re1.Machine2 13.8015  3.7150  
#  Residual               0.9265  0.9626  
# Number of obs: 54, groups:  Worker, 6
# [...]

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:

  • "Die Reduzierung des Modells birgt ein unbekanntes Risiko für Antikonservativität und sollte, wenn überhaupt, mit Vorsicht durchgeführt werden." und
  • "Mein Hauptanliegen ist, dass die Menschen die mit der Modellreduzierung verbundenen Risiken verstehen und dass die Minimierung dieses Risikos einen konservativeren Ansatz erfordert als allgemein üblich (z. B. jede bei 0,05 getestete Steigung)."

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.

Henrik
quelle
IN ORDNUNG. Dann könnte ich einige kleine Änderungen / Aktualisierungen einfügen, um einige der Punkte zu verdeutlichen, die Sie machen. In Bezug auf Bates Preprint kann es in verschiedener Hinsicht sehr suboptimal sein. Aber ich stimme Bates et al. dass singuläre Kovarianzmatrizen genau das gleiche Problem sind wie die Korrelationen von + 1 / -1. Mathematisch gibt es einfach keinen Unterschied. Also , wenn wir die perfekte Korrelationen Kompromiss Macht übernehmen, dann müssen wir sehr vorsichtig sein Singular cov. auch ohne explizite Simulationen. Ich bin nicht der Meinung, dass es "prinzipienlos" ist.
Amöbe
@amoeba lmer_altfunktioniert im Grunde genau wie lmer(oder sogar glmer) mit dem einzigen Unterschied, dass es die ||Syntax erlaubt . Ich bin mir also nicht sicher, warum Sie dies unbedingt vermeiden möchten afex. Es sollte sogar ohne Anbringen funktionieren (dh afex::lmer_alt(...)).
Henrik
@amoeba Was es tut, ist im Grunde der Ansatz, der im Code von Reinhold Kliegl beschrieben wird (dh die zufälligen Effekte erweitern). Für jeden Zufallseffektterm der Formel wird eine Modellmatrix erstellt (dh die Faktoren werden in numerische Kovariaten umgewandelt). Diese model.matrix ist dann cbindzu 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
Henrik
In Bezug auf kategoriale Variablen links von ||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 von lmer_altin afex. Der Vollständigkeit halber lmermö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.
Amöbe
2
@amoeba Es ist wichtig zu beachten, dass der verwendete Ansatz 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 den lmer_altAnruf mit dem mixedAnruf vergleichen.
Henrik
1

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

user55193
quelle