Wahlpflichtfach “Kognitions- und Klinisch-Psychologische Forschungsmethoden”
Wintersemester 2025/2026
Grundlagen von Linear Mixed-Models (LMMs)
Markus Janczyk und Valentin Koob
AG Psychologische Forschungsmethoden und Kognitive Psychologie
Institut für Psychologie, Universität Bremen
Willkommen bei Stats by Randolph. Hier geht es zur Hauptseite mit weiteren Informationen und Inhalten.
Autor:innen dieser Seite: An den Inhalten dieser Seite haben mitgearbeitet: Valentin Koob und Markus Janczyk. Der Inhalt dieser Seite wird in der Lehre in den Studiengängen Psychologie von der AG Psychologische Forschungsmethoden und Kognitive Psychologie an der Universität Bremen verwendet, steht aber allen Interessierten zur Verfügung. Rückmeldungen/Fehler/Vorschläge können gesendet werden an randolph@uni-bremen.
Versionshistory:
- v2.01: minimal korrigierte Version (29.9.2025)
- v2.0: zweite überarbeitete Version (5.9.2025)
- v1.01: kleinere Korrektur bei fixed-/random-factors bzw effects (30.12.2024)
- v1.0: erste vollständige Version (16.12.2024)
1 Einführung
In diesem Teil beschäftigen wir uns einführend mit sog. “linear mixed-effects models” (je nach Kontext auch: multilevel, hierarchical linear oder random-effects models; abgekürzt werden sie in der Literatur auch als LME oder LMM), einem fortgeschrittenen und durchaus komplexem Verfahren zur Erweiterung der Regression. Wir bezeichnen sie hier kurz als LMMs. Während wir hier konzeptuelle und mathematische Grundlagen sowie die Auswertung eines sehr einfachen Designs behandeln, findet sich hier ein Dokument mit zwei weiteren beispielhaften Auswertungen komplexerer Designs bzw. Fragestellungen.
Kernelement dieses Verfahrens ist, dass im Gegensatz zur Regression i.S. des Allgemeinen Linearen Modells (ALM) Datenstrukturen betrachtet werden können, die in gewisser Weise “gruppiert” sind bzw. “auf mehreren Ebenen” betrachtet werden können. Dadurch können Abhängigkeiten von Daten gezielt mit-modelliert werden. Was dies genau bedeuten kann, wird gleich an Beispielen erläutert. Insbesondere in den letzten 10 bis 15 Jahren sind diese Verfahren durch leichter zugängliche Softwarepakete und die Rechenleistung moderner Computer sehr populär geworden. Mitunter werden diese Verfahren auch den Standardverfahren (wie z.B. der within-subject ANOVA) als überlegen angesehen.
In Kapitel 2 werden wir mit einigen grundlegenden Begriffen im Kontext von LMMs beginnen. Dies umfasst v.a. den Aspekt der gruppierten Daten sowie fixed- vs. random-effects. Nach einigen eher mathematischen Details (u.a. zur Parameterschätzung) in Kapitel 3, widmet sich Kapitel 4 dann eher der praktischen Umsetzung anhand der Simon-Daten aus Teil 1. Kapitel 5 behandelt schließlich noch einige Ergänzungen, inklusive dem Schritt zu sog. generalized linear mixed-effects models, abgekürzt oft als GLMMs.
Viele der hier dargestellten Informationen über (G)LMMs basieren auf Fox (2015), Galecki und Burzykowski (2013), Singmann und Kellen (2019) und West et al. (2022), die auch geeignete Einführungen bzw. fortgeschrittene Lektüre darstellen.
Im Rahmen dieser Einführung arbeiten wir mit den folgenden Paketen:
# notwendige Pakete
library(car)
library(afex)
library(lme4)
library(pbkrtest)
library(lmerTest)
library(conflicted)Anmerkungen zu den verwendeten Paketen: Wir laden hier ein Paket namens conflicted. Dieses Paket bietet im eigentlichen Sinne keine Funktionen für das aktuelle Kapitel, aber es “trackt” ob es mehrere Funktionen gibt, die zwar gleich heißen, aber aus unterschiedlichen Paketen kommen. Der Grund dafür, warum wir es gerade hier nutzen ist, dass afex und lme4 beide die Funktion lmer() bereitstellen (wobei afex auf die lmer() Funktion des Pakets lmerTest() verweist). Um zu vermeiden, dass man fälschlicherweise Funktionen aus einem anderen Paket benutzt als beabsicht, produziert das conflicted Paket immer dann eine Fehlermeldung, wenn ein Funktionsaufruf zu mehreren geladenen Paketen passt. In diesen Fällen muss man mit Hilfe einer ::-Schreibweise das gewünschte Paket explizit adressieren (siehe weiter unten für Beispiele).
2 Grundlagen von LMMs
Wir beginnen dieses Kapitel mit einer Auflistung von Anwendungsbeispielen und verdeutlichen im Anschluss die grundlegende Idee von LMMs. Erst danach tauchen wir tiefer in die mathematischen Details ein.
2.1 Gruppierte Daten und interindividuelle Unterschiede
In Statistik 1 und 2 bzw. im Seminar zum ALM haben wir bereits ausführlich über die Regression gesprochen. Ziel der Regression ist es, den Zusammenhang zwischen einem oder mehreren Prädiktoren und einem Kriterium zu beschreiben, wobei der Beitrag der einzelnen Prädiktoren \(X_j\) durch die Regressionsgewichte \({\hat \beta}_j\) quantifiziert wird. Es ist wichtig zu verstehen, dass die Zusammenhänge, die durch \({\hat \beta}_j\) quantifiziert werden, für alle Beobachtungen (bei uns oft Personen) in gleicher Weise angenommen werden:
\[Y_i = {\hat \beta}_0 + \sum_{j=1}^p {\hat \beta}_j X_j + \varepsilon_i\,.\] In dieser Schreibweise ist \({\hat \beta}_0\) der Achsenabschnitt. Auch dieser beeinflusst natürlich im Zuge des ALM die Modellvorhersage für jede Beobachtung. Darüber hinaus geht das ALM davon aus, dass alle Fehler \(\varepsilon_i\) paarweise unkorreliert sind, so dass letztlich auch die Beobachtungen \(Y_i\) unkorreliert sind.
Im Bereich des ALM hatten wir darauf hingewiesen, dass diese Annahme sämtlichen weiteren Ableitungen von Verteilungen zugrunde liegt und dies insbesondere auch für die Verteilung gilt, die für die inferenzstatistischen Tests der Koeffizienten wichtig sind, so wie dies in Abbildung 2.1 zusammengefasst ist.
Abbildung 2.1: Illustration welche Verteilungen aus der Verteilungsannahme im ALM abgeleitet werden können.
Die klassische Regression i.S. des ALM stößt jedoch aus zwei Gründen schnell an ihre Grenzen, sobald die Daten in irgendeiner Art und Weise gruppiert sind:
- Wenn die Daten gruppiert sind, ist es plausibel anzunehmen, dass Zusammenhänge zwischen Prädiktoren und dem Kriterium in den verschiedenen Gruppen variieren. Beispielsweise kann ein Prädiktor, der für eine Gruppe einen starken Einfluss hat, für eine andere Gruppe einen schwachen Einfluss haben. Eine normale Regression würde darüber hinweg gehen und den Einfluss über die Gruppen hinweg quantifizieren. Eine Anmerkung dazu: Tatsächlich ist dies konzeptuell ähnlich zu einer Interaktion im ALM. Im Zuge von LMMs wird allerdings nicht pro Gruppe ein eigener Zusammenhang (isoliert) geschätzt. Vielmehr wird davon ausgegangen, die Zusammenhänge pro Gruppe würden aus einer Verteilung über alle Gruppen hinweg stammen.
- Wenn die Daten gruppiert sind, sind die Beobachtungen innerhalb einer Gruppe ähnlicher als zwischen den Gruppen. Dies führt zwangsläufig zu korrelierten Werten (relativ gesehen zu allen Datenpunkten), was eben genau gegen die Annahme des ALM verstößt.
An dieser Stelle setzen LMMs an: Sie erlauben es, eine Gruppierungsstruktur explizit zu modellieren, um die Einflüsse der Prädiktoren in Form von Regressionkoeffizienten (sog. “effects”) je Gruppe zu variieren und so Korrelationen zwischen den Beobachtungen angemessen zu berücksichtigen. Wichtig ist zu betonen, dass Gruppe hier nicht unbedingt als “zusammengesetzt aus mehreren Personen” zu verstehen ist; tatsächlich wären in within-subject Designs die Daten (verschiedener Bedingungen) nach jeder einzelnen Versuchsperson gruppiert.
2.1.1 Beispiele für gruppierte Daten
Es gibt eine Vielzahl von Fällen, in denen gruppierte Daten in der Psychologie auftreten können, und LMMs sind nicht auf eine bestimmte Art der Gruppierung beschränkt. Die folgenden Beispiele stellen daher nur einen Auszug aus einer Vielzahl möglicher Anwendungskontexte dar:
(Einfache) within-subject Designs: Insbesondere in der Kognitionspsychologie kommt es häufig vor, dass eine Vielzahl von Beobachtungen an jeder Versuchsperson erhoben werden. So kommt es z.B. durchaus vor, dass man im Labor mehrere hundert Durchgänge erhebt, wobei jede Beobachtung eine Reaktionszeit (RT) darstellt (z.B. im Zuge einer Simon-Aufgabe). Im Gegensatz zum Längsschnittdesign spielt es jedoch keine Rolle, wann genau eine Beobachtung erhoben wurde. Die Messwiederholung dient lediglich dazu, mehr Beobachtungen und damit genauere Schätzungen (z.B. für einen Mittelwert) zu erhalten. In diesem Fall ist also nur die Gruppierung der Daten innerhalb einer Versuchsperson relevant, um der Tatsache Rechnung zu tragen, dass sich die RTs (z.B. in ihrer Länge) von Versuchsperson zu Versuchsperson unterscheiden können. Ein Sonderfall, der hier jedoch nur erwähnt und nicht weiter betrachtet werden soll, ist, dass eine Gruppierung auch in Bezug auf das verwendete Stimulusmaterial erfolgen kann. Ein Beispiel wäre hier die Verwendung bestimmter Sätze im Rahmen einer sprachpsychologischen Untersuchung, die jemand für ein Experiment entworfen hat. Diese Sätze sind aber nur eine Auswahl aus vielen möglichen Sätzen, und nicht alle sind gleich schwierig und sie haben daher möglicherweise unterschiedlich starke Auswirkungen auf eine RT. In solchen Fällen ist es naheliegend, dass der Einfluss eines untersuchten Prädiktors von Satz zu Satz unterschiedlich ist, so dass letztlich eine Gruppierung der Beobachtungen sowohl nach Stimuli als auch Versuchspersonen erforderlich ist.
Längsschnittstudien: In der Entwicklungspsychologie aber auch in der Klinischen Psychologie oder Gesundheitspsychologie kommt es häufig vor, dass ein Merkmal im Längsschnitt über mehrere Zeitpunkte gemessen wird: Beispielsweise können bestimmte Patient:innengruppen hinsichtlich ihrer Symptomatik über mehrere Wochen vor und nach einem kritischen Ereignis untersucht werden. In diesem Fall gibt es zwei Gruppierungsebenen von Beobachtungen. Erstens liefert jede:r Patient:in mehrere Beobachtungen, so dass die Daten pro Patient:in gruppiert sind. Zweitens können sich Patient:innengruppen hinsichtlich der Schwere ihrer Symptome ähnlicher sein, so dass die Daten auch nach Patient:innengruppen gruppiert sind. Ähnliche Beispiele gibt es in der Gesundheitspsychologie. Beispielsweise könnte die Händehygiene in einem Krankenhaus zu verschiedenen Zeitpunkten vor und nach dem Ausbruch der Corona-Pandemie in Abhängigkeit von verschiedenen Berufsgruppen betrachtet werden.
Hierarchische Designs: Manchmal findet die Datenerhebung oder Stichprobenziehung auf verschiedenen “Ebenen” statt. Ein klassisches Beispiel hierfür ist, dass zunächst Schulen aus einer größeren Anzahl von Schulen nach dem Zufallsprinzip ausgewählt werden und anschließend eine Stichprobe von Schüler:innen pro Schule gezogen wird. In diesem Fall gibt es zwei “Ebenen” bzw. “Gruppierungsebenen”: Eine erste “übergeordnete” Gruppierung auf Schulebene (“Ebene 2”) und eine zweite Gruppierung innerhalb der Schulen auf Ebene der Schüler:innen (“Ebene 1”).
2.1.2 Zu den Begriffen fixed- vs. random-factor bzw. -effects
Um die nachfolgenden Erläuterungen zu erleichtern, führen wir hier bereits zwei wesentliche Begriffspaare (konzeptuell) ein, die später noch genauer erläutert werden: fixed- vs. random-factors und die mit ihnen assoziierten fixed- vs. random-effects:
Fixed-factors sind die typischen Faktoren, wie sie z.B. im Bereich der Varianzanalyse (ANOVA) verwendet werden: Nominalskalierte Variablen, bei denen die Stufen alle für die Studie/das Experiment relevanten Bedingungen umfassen. Im Bereich des ALM sind dies dann alle nominalskalierte Prädiktoren, die ggf. über Dummy- oder Effektkodierung in das Modell aufgenommen werden. Sogenannte fixed-effects sind dann die Regressionskoeffizienten, die für alle Beobachtungen (und damit auch für die gesamte dahinter stehende Grundgesamtheit) gelten sollen. In gewisser Weise sind die bereits bekannten Regressionsgewichte \({\hat \beta}_j\) im Rahmen des ALM daher fixed-effects, da sie für alle Beobachtungen identisch sind. Sie repräsentieren die aus der klassischen ANOVA bzw. Regression bekannten “Einflüsse”, “Haupteffekte”, “Interaktionen” etc. und ihre inhaltliche Interpretation hängt entsprechend von den betrachteten Prädiktor(-kombination)en ab.
Random-factors sind die diejenigen Faktoren, deren Stufen in einer bestimmten Studie durch eine Zufallsauswahl bestimmt worden sind. Der bekannteste random-effect sind oft die Versuchspersonen, die eine Zufallsauswahl aus einer interessierenden Population darstellen. Würde die Studie wiederholt werden, werden andere Versuchspersonen gezogen und daher ändern sich die Stufen dieses random-factors (während die fixed-effects gleich bleiben). Anders gesagt sind random-factors diejenigen, bzgl. derer man Inferenzen auf die dahinterliegende Populationen machen möchte. Im Zuge von LMMs, können die Einflüsse der Stufen der Prädiktoren allerdings unter Berücksichtigung der Gruppierungen variieren. Diese für jede Gruppe einzigartigen “Effekte” bzw. Koeffizienten werden relativ zu den fixed-effects gesehen, und als random-effects bezeichnet.
2.1.3 Interindividuelle Unterschiede
Im Folgenden betrachten wir ein einfaches Beispiel, in welchem Daten nach Versuchspersonen gruppiert werden; wir fokussieren also auf den Einfluss interindividueller Unterschiede in einem Experiment (das Beispiel ist angelehnt an Singmann & Kellen, 2019). Gegeben sei ein Experiment, in dem RTs auf Stimuli erhoben wurden, die entweder “kongruent” oder “inkongruent” sind (z.B. im Rahmen einer Simon-Aufgabe oder einer Flanker-Aufgabe). Insgesamt gibt es \(K = 3\) Versuchspersonen, wobei jede Versuchsperson \(n_j = 50\) Beobachtungen liefert. Die Daten befinden sich im Datensatz dummy_IV.csv (der hier heruntergeladen werden kann; sie werden aber erst in den folgende Abschnitten dann tatsächlich ausgewertet):
library(car)
data <- read.table("./Daten/dummy_IV.csv",
header = T,
sep = ";"
)
data$Cond <- as.factor(data$Cond)
contrasts(data$Cond) <- contr.sum(2)
some(data, n = 5) # aus carAbbildung 2.2 zeigt die Mittelwerte der drei Versuchspersonen über beide Bedingungen als schwarze Kästchen (in der Legende “beo.” für “beobachtet”) und visualisiert zusätzlich eine Reihe von Modellvorhersagen als graue Punkte (in der Legende “vor.” für “vorhergesagt”), um uns die Funktionsweise von LMMs näher zu bringen.
Abbildung 2.2: Illustration verschiedener Modelle mit random intercepts und random slopes.
Grundsätzlich lässt sich rein deskriptiv feststellen, dass die Versuchspersonen sich in ihrem allgemeinen RT-Niveau unterscheiden. Beispielsweise hat VP1 im Allgemeinen kürzere RTs als VP2. Alle Versuchspersonen zeigen kürzere RTs in kongruenten als in inkongruenten Bedingungen, aber die Versuchspersonen unterscheiden sich auch in der Differenz zwischen kongruenten und inkongruenten Durchgängen. Bei VP1 ist der Unterschied, also der Kongruenzeffekt, z.B. deutlich kleiner als bei VP3.
Das erste Modell (oben links) schätzt einen gemeinsamen Mittelwert (Grand Mean) für alle Versuchspersonen.
Das zweite Modell (oben rechts) schätzt den Gesamtmittelwert sowie eine Abweichung vom Gesamtmittelwert, um den Unterschied zwischen den beiden Bedingungen zu berücksichtigen (ein solches Modell wäre auch im Zuge des ALM möglich). Diese Schätzung ist allerdings für alle Versuchspersonen gleich, so dass es sich hierbei nur um den entsprechenden fixed-effect handelt, d.h. um den Teil der Regression, der für alle Versuchspersonen gleich ist und der damit interindividuelle Unterschiede ignoriert.
Im dritten Modell (unten links), welches nun ein LMM ist, kommen zusätzlich interindividuelle Unterschiede ins Spiel. Hier wird sowohl ein Gesamtmittelwert, eine allgemeine Differenz zwischen den Bedingungen, als auch ein individuelle Abweichung vom mittleren Niveau geschätzt (ein sog. “random intercept”).
Im letzten Modell (unten rechts) werden schließlich noch interindividuelle Unterschiede in der Differenz zwischen den beiden Bedingungen berücksichtigt (ein sog. random slope-Modell, wobei der slope hier die halbe Differenz zwischen den Bedingungen ist).
Man erkennt also, dass ein LMM über das ALM hinausgeht, welches im Rahmen einer Regression (hier ANOVA bzw. \(t\)-Test) keine interindividuellen Effekte bzw. Unterschiede modellieren könnte. Die Hinzunahme von random intercepts und random slopes verbessert letztlich auch die Vorhersagekraft (die über die vier Modelle steigt).
2.2 Grundlegende Aspekte von LMMs
2.2.1 Konzeptuelle Erklärung
In diesem Abschnitt wollen wir konzeptuell (daher auch etwas vereinfacht) erläutern, wie LMMs mit den interindividuellen Unterschieden im generellen RT-Level (random intercept) und der Größe des Effektes (random slopes) umgehen. Der linke Teil von Abbildung 2.3 ist dabei unser illustratives Datenbeispiel, bei dem wir noch ein paar mehr Versuchspersonen betrachten: Acht Versuchspersonen haben z.B. eine Simon-Aufgabe bearbeitet und es liegen daher (mittlere) RTs für die kongruente und die inkongruente Bedingung vor.
Zunächst betrachten wir nur den Intercept, also die generellen Unterschiede im RT-Level. Das generelle RT-Level kann dann z.B. als der Mittelwert der komgruenten und inkongruenten Bedingung aufgefasst werden. Dies sind die roten Punkte im linken Teil von Abbildung 2.3. Die horizontale rote Linie ist der Mittelwert der roten Punkte, also der fixed-effect des Intercepts. Entsprechend symbolisieren die vertikalen gestrichelt-roten Linie die Abweichungen der individuellen RT-Level vom Intercept (die random-effects des Intercepts). Ein erster Ansatz wäre es, diese Abweichungen für jede Versuchsperson als separaten Parameter zu behandeln. Allerdings würde dies dann zu einer sehr großen Anzahl an Parametern führen. Gleichzeitig würde ein weiterer Aspekt aus dem Fokus geraten: Die Versuchspersonen sind ja gerade zufällig gezogen worden und die individuellen Abweichungen sollten daher als (zufällige) Realisierungen einer Zufallsvariablen aufgefasst werden. Würden in einer Wiederholung des Experimentes andere Versuchspersonen gezogen werden, kommt es entsprechend zu anderen Realisierungen. Dies unterstreicht noch einmal, dass Versuchspersonen ein random-factor sind, über den man generalisieren möchte.
LMMs nehmen nun an, dass diese individuellen Abweichungen aus einer normalverteilten Zufallsvariable mit \(\mu_I = 0\) und einer unbekannten Varianz \(\sigma_I^2\) stammen.1 Dies ist im rechten Teil von Abbildung 2.3 illustriert. Die roten Punkte im unteren Bereich sind nun genau die gestrichelt-roten Linien im linken Teil dieser Abbildung (also die Abweichungen). Ein LMM fragt nun konzeptuell danach, welche Varianz die darüber eingezeichnete Normalverteilung haben muss, damit diese Abweichungen maximal wahrscheinlich sind. Dies geht beispielsweise durch eine Maximum-Likelihood Schätzung (siehe hier für eine interaktive Illustration), wie wir sie auch bereits in Teil 1 ja kennengelernt haben. Zusammengefasst wird also in diesem Fall nur ein einziger Parameter geschätzt, \(\sigma_I\) bzw. \(\sigma_I^2\), und dies ist unabhängig davon, wieviele Versuchspersonen tatsächlich zu den Daten beigetragen haben.
Abbildung 2.3: Illustration eines fiktiven Datenbeispiels (links) sowie der resultierenden Normalverteilung aus der die random intercepts stammen (rechts).
Eine ganz ähnlich Logik gilt nun für die random-slopes. Der linke Teil von Abbildung 2.4 stellt zunächst wieder die (mittleren) RTs der \(8\) Versuchspersonen dar. Die gestrichelt-blauen vertikalen Linien sind hier also die individuellen (Kongruenz-)Effekte der Versuchspersonen. Diese sind im mittleren Teil der Abbildung als blaue Punkte eingezeichnet. Die rote horizontale Linie ist dann der mittlere (Kongruenz-)Effekt, also der fixed-effect des Prädiktors Kongruenz, und die gestrichelt-roten vertikalen Linien sind dann die Abweichungen der individuellen Effekte vom fixed-effect Kongruenz. Der rechte Teil von Abbildung 2.4 stellt diese Abweichungen dann wieder als rote Punkte dar und visualisiert die Normalverteilung aus der diese Abweichungen stammen. Mit \(\mu_S=0\) ist hier also \(\sigma_S\) bzw. \(\sigma_S^2\) gesucht.
Abbildung 2.4: Illustration der random-slope Werte, für die die Varianz einer Normalverteilung bestimmt wird.
Schließlich kann es noch sein, dass es Kovariationen zwischen den Intercepts und Slopes der Versuchspersonen gibt. Dies ergibt sich dadurch, dass wir pro Versuchsperson sowohl eine Abweichung vom fixed-effect des Intercepts als auch vom fixed-effect der Kongruenz haben und beide Abweichungen können über die Versuchspersonen hinweg kovariieren: Beispielsweise könnte eine Versuchsperson mit generell eher längeren RTs auch einen größeren Kongruenzeffekt aufweisen als eine Versuchsperson mit generell eher kürzeren RTs. In diesem Fall würde es eine positive Kovarianz (bzw. Korrelation) \(\rho_{IS}\) geben.
Neben der Schätzung der beiden fixed-effects im aktuellen Beispiel, also \(\beta_0\) und \(\beta_1\) im klassischen ALM, sowie der Residualvarianz \(\sigma^2_{\varepsilon}\), kommen hier also drei (zu schätzende) Parameter im LMM dazu: \(\sigma^2_I\), \(\sigma_S^2\) und \(\rho_{IS}\). Wir werden später sehen, dass dies dadurch zum Ausdruck kommt, dass wir eine weitere multivariat-normalverteilte Zufallsvariablen einführen werden, deren Kovarianzmatrix die Varianz-Kovarianz-Struktur der random-effects modelliert. Im Allgemeinen wird für jeden random-effect ein (Varianz-)Parameter zusätzlich geschätzt. Allerdings nimmt die Zahl der zu schätzenden Kovarianzparameter schneller zu. Im Beispiel gab es bei zwei random-effects genau eine Kovarianz. Im Allgemeinen Fall gilt aber, dass es bei \(r\)-vielen random-effects \(\frac{r(r-1)}{2}\)-viele Kovarianzen gibt. Bei \(r=5\) random effects sind dies dann bereits 10 Kovarianzparameter zusätzlich.
2.2.2 Erklärung am Datenbeispiel
Um die folgende Syntax für LMMs besser zu verstehen, schauen wir uns zunächst noch einmal kurz den klassischen Aufruf im Rahmen des ALM mit Hilfe der Funktion lm() an, um deutlich zu machen, was die Koeffizienten bedeuten und woher sie kommen. Da wir hier jedoch grob vernachlässigen, dass eine Person mehrere RTs beiträgt, dient diese Vorgehensweise nur zu Demonstrationszwecken:2
# 1 + nicht nötig, aber zur Demonstration ausgeschrieben
lm_out <- lm(RT ~ 1 + Cond,
data = data
)
brief(lm_out)## (Intercept) Cond1
## Estimate 500.0 -48.5
## Std. Error 4.2 4.2
##
## Residual SD = 51.4 on 148 df, R-squared = 0.474Als Ergebnis erhalten wir hier die Koeffizienten \({\hat \beta}_0 = 500\) und \({\hat \beta}_1 = -48.5\). Da wir hier einen kategorialen Prädiktor betrachten, müssen wir um die Koeffizienten zu interpretieren uns erinnern, dass die einzelnen Stufen eines kategorialen Prädiktors in Form von Dummyvariablen kodiert werden (siehe kategoriale Prädiktoren, Statistik 2):
contrasts(data$Cond)## [,1]
## comp 1
## incomp -1Zu erkennen ist, dass hier eine sog. Effekt-Kodierung vorliegt (vgl. auch den constr.sum(2) beim Laden der Dateien). In diesem Fall entspricht \({\hat \beta}_0 = 500\) der mittleren RT aller Beobachtungen (= über alle Versuchspersonen hinweg) und stellt damit den Grand Mean dar. Der Steigungskoeffizient \({\hat \beta}_1 = -48.5\) ist die Differenz zwischen den mittleren RTs kongruenter Durchgängen und dem Grand Mean. Die Differenz zwischen kongruenten und inkongruenten Durchgängen, also der Kongruenzeffekt, ist das doppelte davon (also \(2\cdot {\hat \beta}_1 = -97\)). Da wir hier ein allgemeines RT-Niveau und eine Differenz zwischen den Bedingungen schätzen, aber eben über alle Versuchspersonen hinweg, entspricht diese Regression dem zweiten Modell (oben rechts) in Abbildung 2.2.
Um nun zusätzlich die interindividuellen Unterschiede im RT-Level und in der Größe des Kongruenzeffektes zu modellieren, müssen wir ein LMM formulieren und schätzen. In R gibt es u.a. hierfür die Funktion lme() aus dem nlme Paket bzw. die Funktion lmer() aus dem Paket lme4. Erstgenannte Funktion ist älter und in gewisser Weise aber auch mächtiger. Die zweite Funktion ist allerdings effizienter und leichter handzuhaben wenn es mehrere Gruppierungsvariablen gibt. Auch scheint uns diese Funktion in der Psychologie verbreiteter im Einsatz zu sein und aus diesem Grund werden wir uns i.W. auf die Funktion lmer() beschränken.
Der grundlegende Aufruf ist hierbei ähnlich zur lm() Funktion, mit dem Unterschied, dass wir explizit die Schätzung von random-effects in einer Klammerschreibweise angeben können. Die Syntax für ein Modell, welches zusätzlich individuelle RT-Level der Versuchsperson schätzt, sieht so aus (beachten Sie hier, dass durch die ::-Schreibweise explizit die Funktion lmer() aus aus dem Paket lme4 aufgerufen wird):
library(lme4)
# lme4:: um explizit die Funktion aus dem lme4 Paket zu nutzen
# man kann das Laden des Paketes dann auch weglassen, aber es
# werden im Folgenden noch weitere Funktionen daraus benötigt
lme_out <- lme4::lmer(RT ~ 1 + Cond + (1 | VP),
data = data
)Der erste Teil der Formelschreibweise ist identisch zum Aufruf von lm(). Neu hingegen ist die zusätzliche (additive) Angabe des Terms (1 | VP), welcher anzeigt, dass pro Versuchsperson Unterschiede im Intercept/Achsenabschnitt, also dem individuellen RT-Level berücksichtigt werden sollen; ein sog. random intercept Modell (siehe die Darstellung links unten in Abb. 2.2). Detailliertere Informationen zu den Ergebnissen fordern wir wie gewohnt mit der summary() Funktion an:
summary(lme_out)## Linear mixed model fit by REML ['lmerMod']
## Formula: RT ~ 1 + Cond + (1 | VP)
## Data: data
##
## REML criterion at convergence: 1372.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.08120 -0.79378 -0.04531 0.79521 2.73426
##
## Random effects:
## Groups Name Variance Std.Dev.
## VP (Intercept) 3117 55.83
## Residual 540 23.24
## Number of obs: 150, groups: VP, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 499.993 32.287 15.49
## Cond1 -48.513 1.897 -25.57
##
## Correlation of Fixed Effects:
## (Intr)
## Cond1 0.000Hier finden wir sehr viele Informationen die gegenüber der klassischen lm()-Funktion im Kontext des ALM neu sind. Wir werden im Laufe dieses Teils noch alle Details besprechen, wollen uns aber zunächst auf die Kerninformationen beschränken. Zwei wichtige Teile des Outputs stehen unter Random effects und Fixed effects:
- Der Teil
Fixed effectslistet die Koeffizienten, \({\hat \beta}_0\) und \({\hat \beta}_1\), der Regression auf, die für alle Versuchspersonen gleichermaßen gelten und diese decken sich (hier) mit den Koeffizienten aus dem Aufruf vonlm(). Man kann diese auch explizit mit der Funktionfixef()(aus dem Paketlme4) extrahieren:
fixef(lme_out)## (Intercept) Cond1
## 499.99333 -48.51333- Der Teil
Random effectsgibt Varianzen bzw. Standardabweichungen sowie die Anzahl der Gesamtbeobachtungen und die Anzahl der Gruppen pro Gruppierungsfaktor an (hier ist “Versuchsperson” [VP] der Grupperierungsfaktor für den wir random-effects zulassen). Die Varianzen beziehen sich dabei einerseits auf die Variabilität im Achsenabschnitt (Intercept, also dem allgemeinen RT-Level) als auch auf die Residualvarianz.
Um genauer zu klären, was ein random-intercept ist und wie dies mit der angegebenen Varianz zusammenhängt, schauen wir uns die berechneten random intercepts, \(\hat b_{i,0}\), mit der Funktion ranef() einmal an:
ranef(lme_out)$VPWas sind nun diese Werte? Diese ausgegebenen Werte sind individuelle Abweichungen im Achsenabschnitt (hier i.S. eines allgemeinen RT-Levels) relativ zum fixed-effect Achsenabschnitt, also dem Grand Mean von \({\hat \beta}_0 = 500\). Dies wird auch nochmal klar, wenn wir die Summe aus random intercept pro Versuchsperson und fixed intercept, \(\hat b_{i,0}\) + \({\hat \beta}_0\), einmal explizit ausrechnen und mit den RT-Mittelwerten pro Versuchsperson vergleichen:
# Mittelwert pro Versuchsperson
aggregate(RT ~ VP,
data = data,
FUN = mean
) ranef(lme_out)$VP + fixef(lme_out)[1] # random intercepts plus fixed effect interceptZu sehen ist, dass beide Varianten sehr ähnlich ausfallen, allerdings gibt es auch kleine Abweichungen: Die Summe aus random intercept und fixed intercept fällt etwas höher aus für Werte unter 500, und etwas niedriger für Werte über 500. In anderen Worten: Die Summe ist leicht “in Richtung 500”, also in Richtung des fixed intercept, \({\hat \beta}_0\), “gebiased”. Dieses Phänomen wird als “shrinkage” bezeichnet und ist ein Nebenprodukt der Schätzung der random-effects (mehr dazu später). Die angegebene Varianz des random intercepts im summary()-Output weiter oben ist dabei ein Varianzschätzer. Diese Varianz ist aber nicht die Varianz der random effects an sich, sondern die geschätzte Varianz einer Normalverteilung (mit Erwartungswert \(0\)) aus der die random intercepts stammen sollen, und im aktuellen Beispiel somit ein Varianzschätzer der individuellen RT-Levels unter Annahme derer Normalverteiltheit (mehr zu den Annahmen später).
Die Residualvarianz im summary() Output ist letztlich ähnlich zu der Residualvarianz im Zuge des ALM und beschreibt wie sehr die einzelnen Datenpunkte von der Modellvorhersage abweichen. Im Kontext von LMMs besteht die Modellvorhersage aber eben aus random- und fixed-effects und somit beschreibt die Residualvarianz die Varianz der Residuen innerhalb jeder Gruppierungsstufe (hier also innerhalb jeder Versuchsperson). Abbildung 2.5 zeigt die Rohdaten in Form von roten Punkten sowie die hieraus berechneten Mittelwerte pro Bedingung und Versuchsperson (in schwarz). Die Residuen pro Versuchsperson sind die Abweichungen der roten Punkte von den Modellvorhersagen (hier in grau).
Abbildung 2.5: Illustration der Beispieldaten und der Vorhersagen des random intercept Modells.
Wie wir auch schon im einführenden Beispiel gesehen haben, ist ein Modell mit einem random intercept, also mit individuell verschiedenen RT-Levels schon besser als ein Modell ohne interindividuelle Variabilität. Es ist allerdings auch noch nicht ganz ideal, da sich die einzelnen Versuchspersonen auch bezüglich der Größe des Kongruenzeffektes unterscheiden. Aus diesem Grund sind die Residuen pro Bedingung und Versuchsperson auch mal eindeutig negativ und mal eindeutig positiv. Um diese interindivduellen Unterschiede ebenfalls zu berücksichtigen, erweitern wir das Modell entsprechend und lassen nun auch noch einen random slope zu, wofür wir die Klammerschreibweise erweitern:
lme_out <- lme4::lmer(RT ~ 1 + Cond + (1 + Cond | VP),
data = data
)
summary(lme_out)## Linear mixed model fit by REML ['lmerMod']
## Formula: RT ~ 1 + Cond + (1 + Cond | VP)
## Data: data
##
## REML criterion at convergence: 1232.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.75339 -0.64083 -0.07181 0.65181 3.01885
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## VP (Intercept) 3123.5 55.89
## Cond1 503.9 22.45 -0.52
## Residual 194.9 13.96
## Number of obs: 150, groups: VP, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 499.99 32.29 15.486
## Cond1 -48.51 13.01 -3.729
##
## Correlation of Fixed Effects:
## (Intr)
## Cond1 -0.518Die geschätzten fixed-effects sind hierbei identisch zur vorherigen Ausgabe. Allerdings hat sich der Teil zu den random-effects erweitert. Hier steht nun auch, dass wir einen random slope für den Prädiktor bzw. die Dummyvariable Cond zugelassen haben, also interindividuelle Unterschiede in der Größe des Kongruenzeffektes. Als Konsequenz finden wir hier nun zusätzlich eine geschätzte Varianz der Normalverteilung aus der die random slopes stammen sowie eine geschätzte Korrelation zwischen den random slopes und den random intercepts (die hier negativ ausfällt; aber aufgrund der sehr wenigen Versuchspersonen wohl keine Bedeutung trägt). Hier beziehen sich die Varianzen und Korrelation nicht auf die reine Korrelation und Varianz der einzelnen random effects. Stattdessen stellen diese Parameter der Kovarianzmatrix einer multivariaten Normalverteilung mit Erwartungswertvektor \({\bf 0}\) dar, aus der die einzelnen individuellen Unterschiede am wahrscheinlichsten stammen.
Visualisieren wir die Mittelwerte und Residuen pro Versuchsperson erneut, sehen wir in Abbildung 2.6, dass es keine keine systematischen Unter- bzw. Überschätzungen pro Bedingung und Versuchsperson mehr gibt (vgl. auch die Darstellung rechts unten in Abb. 2.2).
Abbildung 2.6: Illustration der Beispieldaten und der Vorhersagen des random slopes Modells.
Einen Plot der Residuen an sich kann man auch “auf die Schnelle” mit Hilfe der generischen plot() Funktion anfordern (Abb. 2.7):
plot(lme_out)
Abbildung 2.7: Einfache Visualiserung der Residuen mit der plot()-Funktion.
3 Mathematische Grundlagen
Nachdem wir nun LMMs einführend und konzeptuell kennengelernt haben, wollen wir nun mehr auf die mathematischen Grundlagen eingehen. Dabei beschränken wir uns weitgehend auf das sog. klassische LMM mit nur einem Gruppierungsfaktor.
3.1 Das klassische LMM auf Populationsebene
Gegenüber der Schreibweise im Rahmen des ALM gibt es 3 wesentliche Neuerungen:
- Einführung von Gruppen bzw. Einzelschritten pro Gruppe, die entsprechend berücksichtigt werden müssen.
- Eine additive Trennung in fixed- und random-effects.
- Die Möglichkeit, eine bestimmte Kovarianzstruktur der Fehler zu spezifizieren (was nicht immer notwendig und im Paket
lme4auch nicht möglich ist).
Beginnen wir mit einer kurzen Wiederholung der Regressionsgleichung im Rahmen des ALM. Im Allgemeinen liegt für eine beliebige Versuchsperson \(i\) die folgende Regressionsgleichung vor:
\[ Y_i = \beta_1 x_{i,1} + \beta_2 x_{i,2} + \ldots + \beta_m x_{i,m} + \varepsilon_i\,. \] Die klassische Schreibweise mit einem Achsenabschnitt erhält man leicht, indem man \(x_{i,1} = 1\) setzt und die \(\beta\)s von \(0\) bis \(p\) durchnummeriert bzw. die restlichen Prädiktorwerte \(x_{i,2}\) bis \(x_{i,m}\) in \(x_{i,1}\) bis \(x_{i,p}\) umbenennt. Diese Notation ist in gewisser Weise allgemeiner und umfasst sowohl Modelle mit als auch ohne Achsenabschnitt.3
Die Variable \(\varepsilon_i\) ist das Residuum einer Versuchsperson und stellt im stochastischen Sinne eine Zufallsvariable dar, sodass letztlich auch \(Y_i\) eine (linear transformierte) Zufallsvariable ist. Mit Hilfe der Matrizenschreibweise kann dieselbe Information auch kompakt dargestellt werden als \[Y = {\bf X} \cdot {\bf \beta} + \varepsilon\,.\] Wenn \(n\) die Stichprobengröße und \(m\) die Anzahl der Prädiktoren ist, dann sind die Dimensionen der Vektoren bzw. Matrizen dabei wie folgt:
- \(Y\): \(n\times1\)
- \(\bf X\): \(n \times m\)
- \({\bf \beta}\): \(m \times 1\)
- \(\varepsilon\): \(n \times 1\)
Die zentrale Annahme ist schließlich, dass \(\varepsilon\) multivariat normalverteilt ist, wobei alle Varianzen identisch sind und alle Kovarianzen Null sind, also \[\varepsilon \sim N_n({\bf 0}, \sigma^2 {\bf I_n})\,.\]
Aufbauend auf dieser kurzen Wiederholung wollen wir im nächsten Schritt das ALM sukzessive erweitern zum klassischen LMM in der sog. Laird-Ware-Form (weil so von Laird & Ware, 1982, eingeführt), wobei wir nur einen Gruppierungsfaktor betrachten.
Als erstes müssen wir uns klar machen, dass sich alle \(n\) Beobachtungen unter Berücksichtigung der Gruppiertheit der Daten aufteilen lassen in \(K\)-viele Gruppen, jeweils vom Umfang \(n_i\) (bspw. liegen \(K\)-viele Versuchspersonen mit jeweils \(n_i\)-vielen Beobachtungen vor). In diesem Sinne können wir auch pro Stufe der Gruppierungsvariable einen Kriteriumsvektor \(Y_i\), eine Designmatrix \(\bf X_i\) sowie einen Messfehler \(\varepsilon_i\) definieren: \[Y_i = {\bf X_i} {\bf \beta} + \varepsilon_i \quad i\in\{1,..,K\}\,.\] Es ist wichtig, sich bereits jetzt daran zu gewöhnen, dass der Index \(i\) nicht für eine Beobachtung steht, sondern für eine Stufe bzw. Gruppe innerhalb des Gruppierungsfaktors (also z.B. für eine Versuchsperson). Daher liegen nun \(K\)-viele ALM-Gleichungen von oben vor, wobei \(Y_i\), \(\bf X_i\) und \(\varepsilon_i\) jeweils \(n_i\) viele Zeilen haben (ein Beispiel folgt weiter unten). Weiter ist zu beachten, dass \(\bf \beta\) keinen Index hat, da \(\bf \beta\) die fixed-effects darstellt, die für alle Versuchspersonen gleich sind.
Im nächsten Schritt führen wir nun \(q\)-viele random-effects additiv hinzu, wobei wir diese mit \(b_i\) bezeichnen (im Beispiel von oben hatten wir einen random intercept und einem random slope geschätzt; hier wäre also \(q=2\); siehe dazu auch den nächsten Abschnitt). Für jede Gruppe \(i\) lautet die Regressionsgleichung also nun \[\begin{align*} Y_i = {\bf X_i} {\bf \beta} + {\bf Z_i} b_i + \varepsilon_i \quad i\in\{1,\ldots,K\}\,. \end{align*}\] Die Matrix \(\bf Z_i\) hat hierbei die Dimension \(n_i \times q\) und sie enthält (konstante) Kovariaten für die \(i\)-te Gruppe und skaliert die unbekannten \(b_i\) (der Dimension \(q\times 1\)). Wie wir weiter unten sehen werden, gehören die \(q\)-vielen random-effects in \(b_i\) i.d.R. zu entsprechenden (fixed-)effects in \(\bf \beta\), sodass letztlich die Matrix \(\bf Z_i\) durch einen Auszug der Spalten von \(\bf X_i\) konstruiert wird. Auch wichtig ist zu betonen, dass \(b_i\) ein Zufallsvektor ist. Der entsprechende Vektor enthält also keine (!) Koeffizienten, die als zwar unbekannte aber konstante Werte auf Populationsebene angesehen werden können, sondern der Vektor enthält Zufallsvariablen.
Nun führen wir eine Reihe von Annahmen über die Verteilungen von \(b_i\) und \(\varepsilon_i\) ein, die sehr allgemein gehalten sind:
\(b_i\) ist ein mulitvariat-normalverteilter Zufallsvektor mit Erwartungswertvektor \(\bf 0\) und Kovarianzmatrix \(\bf {\mathcal D}\): \[b_i \sim N_q({\bf 0}, {\bf {\mathcal D}})\,.\] Da \(b_i\) ein Zufallsvektor ist werden die random-effects auch nicht im engeren Sinne geschätzt. Vielmehr wird ihre Kovarianzstruktur \({\bf {\mathcal D}}\) geschätzt, welche für jede Gruppe identisch ist. Hier stehen also die Varianzen bzw. Kovarianzen/Korrelationen der random-effects, welche wir in geschätzter Form im
summary()Output vom Anfangsbeispiel gesehen haben.
\(\varepsilon_i\) ist mulitvariat-normalverteilt mit Erwartungswertvektor \(\bf 0\) und Kovarianzmatrix \({\bf {\mathcal R}_i}\): \[\epsilon_i \sim N_{n_i}({\bf 0}, {\bf {\mathcal R}_i})\,.\] Hier besteht ein weiterer Unterschied zum klassichen ALM, bei dem die Fehler die gleiche Varianz besitzen und paarweise unkorreliert sind ( \(\bf {\mathcal R}_i = \sigma^2 {\bf I_{n_i}}\)). Das heißt auch, dass nun korrelierte Residuen mit unterschiedlicher Varianz modelliert werden können. Da in der allgemeinsten Form \(\bf {\mathcal R}_i\) zu viele unbekannte Varianzen und Kovarianzen hätte, wird \(\bf {\mathcal R}_i\) durch einen (kleineren) Satz an Parametern in Kombinationen mit Restriktionen konstruiert. Wir gehen hierauf allerdings nicht genau ein, da
lmer()stets von \(\bf {\mathcal R}_i = \sigma^2 {\bf I_{n_i}}\) ausgeht. Wir greifen das Thema kurz an bei der Besprechung von Längsschnittdaten hier auf und verweisen auf diese mögliche Erweiterung.Schließlich werden eine Reihe von Unabhängigkeitsannahmen formuliert (die durch Erweiterungen teilweise auch fallengelassen werden können/müssen):
\(b_i\) und \(\varepsilon_j\) (auch für \(i = j\)) sind stochastisch unabhängig
\(\varepsilon_i\) und \(\varepsilon_j\) (für \(i\neq j\)) sind stochastisch unabhängig
\(b_i\) und \(b_j\) (für \(i\neq j\)) sind stochastisch unabhängig.
Diese Annahmen bedeuten letztlich, dass die Residuen und random-effects zweier Gruppen jeweils unabhängig sind und auch die random-effects und Residuen innerhalb einer Gruppe unabhängig sind.
3.2 Beispiel: Parametrisierung für eine Bedingung mit random intercept und slope.
Um die Gestalt der Matrizen in der Gleichung \[\begin{align*} Y_i = {\bf X_i} {\bf \beta} + {\bf Z_i} b_i + \epsilon_i \quad i\in\{1,..,K\} \end{align*}\] klarer zu machen, wollen wir hier schematisch das Beispiel zum Kongruenzeffekt aufgreifen. Im entsprechenden Datensatz gibt es drei Versuchspersonen, d.h. es gilt \(K = 3\). Dabei liefert jede Versuchspersonen \(n_i = 50\) Beobachtungen (\(\forall i \in \{1, ..., K\}\)).
table(data$VP)##
## 1 2 3
## 50 50 50Betrachtet wird weiter ein kategorialer Prädiktor in Effektkodierung wobei jede Versuchsperson 25 Beobachtungen pro Stufe liefert:
table(data$VP, data$Cond) # Beobachtungen pro Zelle##
## comp incomp
## 1 25 25
## 2 25 25
## 3 25 25contrasts(data$Cond) # akuell: Effektkodierung## [,1]
## comp 1
## incomp -1Der Knackpunkt ist nun, dass für einen random slope und einen random intercept, \(b_{i,0}\) und \(b_{1,i}\), die Matrix \(\bf Z_i\) tatsächlich identisch aufgebaut ist wie \(\bf X_i\). Daher ergibt sich folgende schematische Darstellung für eine Versuchsperson \(i\):
\[\begin{align*} \underbrace{\begin{pmatrix} Y_{i,1} \\ \vdots \\ Y_{i,25} \\ Y_{i,26} \\ \vdots \\ Y_{i,50} \end{pmatrix}}_{Y_i} = \underbrace{\begin{pmatrix} 1 & 1 \\ \vdots & \vdots \\ 1 & 1 \\ 1 & -1 \\ \vdots & \vdots \\ 1 & -1 \end{pmatrix} \cdot \begin{pmatrix} \beta_0 \\ \beta_1 \end{pmatrix}}_{\bf X_i \bf \beta} + \underbrace{\begin{pmatrix} 1 & 1 \\ \vdots & \vdots \\ 1 & 1 \\ 1 & -1 \\ \vdots & \vdots \\ 1 & -1 \end{pmatrix} \cdot \begin{pmatrix} b_{i,0} \\ b_{i,1} \end{pmatrix}}_{{\bf Z_i} b_i} + \underbrace{\begin{pmatrix} \epsilon_{i,1} \\ \vdots \\ \epsilon_{i,25} \\ \epsilon_{i,26} \\ \vdots \\ \epsilon_{i,50} \end{pmatrix}}_{\epsilon_i} \end{align*}\]
Dabei gehören die ersten 25 Beobachtungen einer Versuchsperson zur Bedingung “kongruent” (vgl. Kodierung im letzten R Output) und die weiteren 25 Beobachtungen zur Bedingung “inkongruent”.
Formulieren wir einzelne Zeilen aus, dann ergibt sich also z.B. \[\begin{align*} Y_{i,1} &= \beta_0 + \beta_1 + b_{i,0} + b_{i,1} + \varepsilon_1\quad\text{für Zeile 1} \\ \text{und }Y_{i,26} &= \beta_0 - \beta_1 + b_{i,0} - b_{i,1} + \varepsilon_{26}\quad\text{für Zeile 26.} \\ \end{align*}\] Insofern ergibt sich die erste Beobachtung einer Person \(i\) in kongruenten Bedingungen durch den personen-unspezifischen Grand Mean, \(\beta_0\) (fixed intercept), die personen-unspezifische Abweichung für kongruente Durchgänge, \(\beta_1\) (fixed slope), die personen-spezifische Abweichung vom Grand Mean, \(b_{i,0}\) (random intercept), die personen-spezifische Abweichung für kongruente Durchgänge, \(b_{i,1}\), und schließlich ein Residuum für die erste Beobachtung, \(\varepsilon_1\). Ganz ähnlich berechnet sich die 26. Beobachtung, also die 1. Beobachtung der inkongruenten Bedingungen, wobei die beiden kongruenzabhängigen Abweichungen nun mit \(-1\) multipliziert werden und entsprechend \(-\beta_1\) und \(-b_{i,1}\) sind.
Man merke, dass hier auf Populationsebene \(b_{i,0}\), \(b_{i,1}\) und \(\varepsilon_1\) Zufallsvariablen sind. Typischerweise werden hierbei für die random-effects keine Restriktion bzgl. der Kovarianzstruktur gemacht, \[\begin{align*} \begin{pmatrix} b_{i,0} \\ b_{i,1} \end{pmatrix} \sim N_2({\bf 0}, {\bf {\mathcal D}}) \quad \text{mit} \quad {\bf {\mathcal D}} = \begin{pmatrix} \Psi_{11} & \Psi_{12} \\ \Psi_{21} & \Psi_{22} \end{pmatrix}\,, \end{align*}\] wohingegen die Messfehler pro Versuchsperson identisch verteilt sein können, oft mit gleicher Varianz und einer Kovarianz von 0, also4 \[\varepsilon_i \sim N_{50}({\bf 0}, \sigma^2 {\bf I_{50}} )\,.\]
3.3 Schätzung des klassischen LMMs
Die Schätzung der Parameter eines LMMs ist sowohl im rein mathematischen Sinn als auch numerisch alles andere als trivial und wir werden hier auch nur die Grundlagen grob klären können. Dennoch lohnt es sich diese einmal anzusehen, um einen Eindruck von der Komplexität zu erhalten sowie ein paar der Begriffe besser einordnen zu können, die einem in der Arbeit mit LMMs begegnen können.
Vorweg sei noch angemerkt, dass die Grundidee derjenigen entspricht die wir in Teil 1 zu mathematischen Modellen und Simulation kennengelernt haben. Das heißt, ein Algorithmus sucht auch hier i.W. diejenigen Parameterwerte, bei denen eine Kostenfunktion ein Minimum annimmt. Die Parameter sind hierbei nun die fixed-effects sowie die Varianz-Kovarianzparameter der random-effects Verteilung und der Residualverteilung.
Ein wenig formalere Ausführungen finden sich hier.
3.3.1 Maximum-Likelihood (ML) versus Restricted Maximum-Likelihood (REML) Schätzung
Das prinzipielle Herangehen entspricht einer Maximum-Likelihood (ML) Schätzung der Parameter (siehe hier für eine interaktive ShinyApp). Im Kontext von LMMs finden aber auch sog. Restricted Maximum-Likelihood (REML) Schätzungen Anwendung und diese Unterscheidung wird auch wichtig für die Wahl inferenzstatistischer Tests über die fixed- (und auch random-)effects (siehe dazu auch McNeish, 2017), die wir in Kapitel 4 behandeln werden.
Bei einer ML Schätzung werden die Varianzkomponenten der random-effects sowie die fixed-effects simultan geschätzt. Das dabei entstehende Problem ist, dass beide voneinander abhängen: Um Variablitäten zu schätzen, müssen wir wissen, um welchen Wert die Datenpunkte variieren sollen, und im Falle von LMMs sind dies eben die fixed-effects. In einem ersten Schritt werden daher die fixed-effects geschätzt, als gäbe es die random-effects nicht. Darauf basierend werden dann die Varianzkomponenten der random-effects für die gerade berechneten fixed-effects geschätzt. Nun werden die fixed-effects erneut geschätzt, dieses Mal allerdings unter Einbezug der random-effects, die im Anschluss auch “aktualisiert” werden. Diese Iterationen wiederholen sich, bis es kaum noch Veränderungen in den Schätzungen gibt. Da aber zur Schätzung der random-effects die fixed-effects in jeder Iteration “feststehen” wird (1) deren Variabilität ignoriert und (2) die verlorenen Freiheitsgrade für deren Schätzung werden nicht berücksichtigt. Genau wie \[ S_X^2=\frac{1}{n}\cdot\sum_{i=1}^n(x_i-M_X)^2 \] den Populationsparameter \(\sigma_X^2\) systematisch unterschätzt (da \(\mathbb{E}(\boldsymbol{S^2_X})=\frac{n-1}{n}\sigma_X^2\) ist; siehe hier)), da der Mittelwert nicht bekannt ist, sondern aus den Daten berechnet werden muss, werden die Varianzkomponenten der random-effects systematisch unterschätzt (insbesondere bei kleinen Stichprobengrößen ist dies der Fall). In die Berechnung der Standardfehler der fixed-effects gehen allerdings die Varianzkomponenten der random-effects mit ein und deren Unterschätzung kann daher zu liberalen Tests und einer erhöhten Rate an \(\alpha\)-Fehlern führen.
Eine REML Schätzung löst dieses Problem, indem die Schätzung der fixed-effects von der Schätzung der Varianzkomponenten der random-effects separiert wird. Dadurch kommt es i.W. zu besseren Schätzungen der Varianzkomponenten (in manchen Fällen zu erwartungstreuen Schätzungen). Etwas vereinfacht gesagt, wird eine lineare Transformation der Daten durchgeführt, in Folge derer die fixed-effects zunächst auf Null gesetzt werden. Dann werden die Varianzkomponenten der random-effects per ML geschätzt, dies geschieht aber nun eben unabhängig von den fixed-effects. Sind dann die Varianzkomponenten bestimmt, werden die fixed-effects durch eine Generalized Least Squares-Methode bestimmt. Diese kommt (wie die OLS-Methode) ohne Iterationen aus, benötigt aber bekannte Varianzkomponenten, wofür eben die Schätzungen aus dem ersten Schritt verwendet werden.
Wir merken uns an dieser Stelle schon einmal, dass zwei Modelle mit verschiedenen fixed-effects bzgl. ihrer REML Likelihood nicht vergleichbar sind, da sie nicht auf den gleichen Restriktionen (bzgl. der fixed effects) basieren und daher Likelihoods für verschiedene Daten sind. Dies gilt insbesondere auch dann, wenn das eine Modell ein Untermodell eines anderen (Ober-)Modells ist.
3.3.2 Schätzung mit dem Package lme4
Wir wollen uns mit den Details der numerischen Vorgehensweise und der Schätzung der Parameter eines LMMs nicht befassen und nur die wichtigsten Punkte hier festhalten:
Die Funktion
lmer()verwendet standardmäßig die REML Schätzung. Wollen wir, was manchmal nötig ist, dennoch eine ML verwenden, muss explizitREML = FALSEgesetzt werden.In der Praxis ist es so, dass die Berechnung eines LMMs oft mit Warnungen einhergeht. Dies liegt daran, dass die Algorithmen Vor- und Nachteile haben. Die Funktionen des
lme4Packages können aber eine ganze Reihe von Algorithmen verwenden und manchmal lohnt es sich dann, einen anderen Optimierungsalgorithmus zu verwenden. Wendet man die FunktionallFit()auf einlmer()-Objekt an, werden bspw. alle Optimierer ausprobiert und man bekommt einen Überblick, mit welchem Optimierer weiter gearbeitet werden kann. Für unser oben betrachtetes Beispiel sähe dies dann so aus (hier gibt es also keinerlei Probleme mit den Optimierern):
allFit(lme_out)## bobyqa : [OK]
## Nelder_Mead : [OK]
## nlminbwrap : [OK]
## nmkbw : [OK]
## nloptwrap.NLOPT_LN_NELDERMEAD : [OK]
## nloptwrap.NLOPT_LN_BOBYQA : [OK]## original model:
## RT ~ 1 + Cond + (1 + Cond | VP)
## data: data
## optimizers (6): bobyqa, Nelder_Mead, nlminbwrap, nmkbw, nloptwrap.NLOPT_LN_NELDERMEAD,nloptw...
## differences in negative log-likelihoods:
## max= 9.37e-08 ; std dev= 3.74e-083.3.2.1 Bestimmung der random-effects
Die bisherigen Ausführungen haben einen Eindruck davon gegeben, wie die zentralen Parameter eines klassischen LMMs bestimmt werden können. Dabei kam aber die Bestimmung der random-effects noch gar nicht vor. Dies wollen wir nun nachholen und grob skizzieren.
Grob gesagt, können die random-effects der Versuchspersonen aus der bedingten Verteilung von \(b_i\) gegeben die beobachteten Daten \(Y_i\) bestimmt werden, also aus \(f_{b_i|y_i}({\bf b_i}| {\bf y_i})\). Diese Verteilung kennen wir zwar nicht direkt, sie lässt sich allerdings durch den Satz von Bayes in Kombination mit Rechenregeln für bedingte Verteilungen herleiten, welche sich nicht nur auf einfache Wahrscheinlichkeiten, sondern auch auf ganze Verteilungen erstrecken: \[ f_{b_i|Y_i}({\bf b_i}| {\bf y_i}) = \frac{f_{Y_i|b_i}({\bf y_i} | {\bf b_i})\cdot f_{b_i}({\bf b_i})}{f_{y_i}({\bf y_i})} = \frac{f_{Y_i,b_i}({\bf y_i}, {\bf b_i})}{f_{y_i}({\bf y_i})}\,. \] Aus diesem Ausdruck sind tatsächlich alle Teile entweder bekannt oder können hergeleitet werden (siehe hier für mehr Details).
Unter wenigen Annahmen ist \(f_{b_i|Y_i}({\bf b_i}| {\bf y_i})\) eine multivariate Normalverteilung und die random-effects der Versuchspersonen werden durch den Schätzer ihres Erwartungswertsvektor bestimmt: \[ \bf \hat b_i = \bf \hat D Z_i'\hat V_i^{-1} (y_i - X_i\bf {\hat \beta})\,.\] Hierbei ist \(\bf \hat V_i\) eine Matrix, die sich auf Basis der geschätzten Kovarianzstruktur der Residuen und random effects ergibt: \[{\bf \hat V_i} = (Z_i \bf {\mathcal{\bf \hat D}} Z_i' + \bf {\mathcal{\hat R}_i})/\hat \sigma^2\]
Da die empirischen random-effects auf Basis der empirischen Schätzer für \(\theta\) und \(\bf \beta\) berechnet werden und im strengen Sinn nicht geschätzt, sondern vorhergesagt werden, nennt man sie auch EBLUPs als Abkürzung fürempirical best linear unbiased predictors.
Ein interessantes Nebenprodukt der Herleitung ist, dass die Varianz der vorhergesagten random-effects auf empirischer Ebene stets kleiner ausfällt als die Varianz der random-effects auf Populationsebene. Es gilt also \[V(\hat b_{i,j}) \leq V(b_{ij})\] bzw. \[\hat \Psi_{jj} \leq \Psi_{jj}\,.\]
Dies gilt nicht nur für jeden einzelnen random-effect, sondern auch für beliebige Linearkombinationen. Da sich die random-effects auf empirischer Ebene nahezu zu 0 addieren bzw. auf Populationsebene einen Erwartungswert von 0 haben, bedeutet dies, dass die empirischen random-effects näher um die \(0\) streuen als man es erwarten würde. In anderen Worten: Die empirischen random-effects werden ein Stück weit in Richtung \(0\) “hingezogen” bzw. die vorhergesagten Werte der Personen werden zu den Modellvorhersagen über alle Versuchspersonen hinweg “hingezogen”, was man als shrinkage bezeichnet. Dieser Effekt ist umso stärker, je größer die Varianz innerhalb jeder Gruppe ist und je kleiner die Stichprobengröße.
Um das Phänomen der shrinkage zu visualisieren, betrachten wir hier einen simulierten Datensatz, bei dem wir die zugrundeliegenden, “wahren” random-effects kennen. Der simulierte Datensatz und die wahren random effects sind in der Datei shrinkage.rds gegeben.
list_in <- readRDS(file = "./Daten/shrinkage.rds")
data <- list_in$data # auspacken der "Rohdaten"
true_ran_eff <- list_in$true_ran_effs # auspacken der wahren random effects
head(data)head(true_ran_eff)## (intercept) quality
## [1,] 1.9553440 1.31953469
## [2,] 0.2757388 0.37071823
## [3,] 0.6380666 3.17381335
## [4,] 0.6816115 1.57945543
## [5,] -1.4789450 -2.07896584
## [6,] 1.2071778 -0.07729409Der Datensatz stellt hierbei simulierte Daten für den Zusammenhang zwischen der Qualität einer sozialen Interaktion (quality) und dem anschließenden empfundenen “Glück” im Sinne des Wohlbefindens einer Person (happiness) zur Verfügung. Im Objekt true_ran_effs befinden sich die tatsächlichen random-effects in Form eines random intercepts und eines random slopes, welche der Simulation der Rohdaten zugrunde gelegt wurden. Wie verhalten sich nun diese “wahren” random-effects zu den vorhergesagten random-effects? Hierfür schätzen wir ein LMM mit einem random intercept und einem random slope und visualisieren im Anschluss die Ergebnisse als Scatterplot (siehe Abb. 3.1). Wir nutzen hier eine Funktion, die netterweise von einem Nutzer auf Github zur Verfügung gestellt wird, siehe hier.
lme_out <- lme4::lmer(happiness ~ 1 + quality + (1 + quality | vp),
data = data
)
# Empirische random effects
eblups <- ranef(lme_out)$vp
# Ursprüngliche random effects und empirische random effects zu einem
# Data.Frame basteln
true_ran_eff <- as.data.frame(true_ran_eff)
colnames(true_ran_eff) <- c("intercept", "slope")
true_ran_eff$orig <- "true"
eblups <- as.data.frame(eblups)
colnames(eblups) <- c("intercept", "slope")
eblups$orig <- "eblup"
joint <- rbind(eblups, true_ran_eff)
# Funktion von Github laden
source("./marginal_plot.R")
marginal_plot(
x = intercept,
y = slope,
group = orig,
data = joint,
bw = "nrd",
lm_formula = NULL,
pch = 15,
cex = 0.5
)## Warning: Paket 'scales' wurde unter R Version 4.4.3 erstellt
Abbildung 3.1: Illustration des shrinkage-Phänomens.
Zu sehen ist, dass die vorhergesagten random-effects auf Basis des geschätzten Modells (rot) deutlich weniger streuen und mehr in eine gemeinsame Richtung “zusammengezogen” sind als die ursprünglichen random-effects. Dies trifft hier vor allem auf den random intercept zu. Inhaltlich ist dies ähnlich wie die Regression zur Mitte, da eine zukünftig gezogene Versuchsperson wahrscheinlich weniger extreme Werte der zugrunde liegenden Normalverteilung der random-effects aufweisen würde.
4 Beispielhafte Auswertung und Inferenzstatistik bei LMMs
Nach den eher konzeptuellen und mathematischen Ausführungen in den Kapiteln 2 und 3 betrachten wir nun einen echten Datensatz, den wir zunächst in der einfachsten Form mit einem LMM auswerten. An diesem Beispiel führen wir dann auch die verschiedenen inferenzstatistischen Methoden ein, die zur Evaluation der fixed- und random-effects herangezogen werden können. Wir nutzen dazu den Simon-Datensatz, den wir auch bereits in Teil 1 verwendet haben.
4.1 Die Auswertung experimenteller Simon-Daten
Wir laden zunächst den Datensatz und führen die gleiche Filterung durch wie in Teil 1:
simon_daten <- read.table("./Daten/daten_Simon.txt",
sep = ";",
header = TRUE
)
# Filtern der Daten: Auschließen der Übungsblöcke, fehlerhafter Trials,
# und zu kurzer/langer RTs ("Ausreißer")
f_simon_daten <- subset(
simon_daten,
blocknr > 0 & error == 0 & errortype == 0 &
rt >= 200 & rt <= 1200
)Wir haben hier bereits festgestellt, dass dieser Datensatz einen Simon-Effekt aufweist, dort beruhte die Auswertung jedoch auf gemittelten Daten. Hier verwenden wir nun die nicht-aggregierten Rohdaten und wenden ein LMM auf diese an. Allerdings hatten wir auch bereits festgestellt, dass RTs nicht normalverteilt sind. Dies stellt zunächst ein Problem dar, da klassische LMMs normalverteilte Residuen annehmen (siehe Kap. 3). Zum Umgang damit gibt es i.d.R. 3 Strategien:
Man hofft auf die Robustheit von LMMs gegenüber leichten Verletzungen der Residualverteilungsannahme (vgl. bspw. Schielzeth et al., 2020).
Man transformiert das Kriterium, bspw. mit Hilfe des Logarithmus, um näher an eine Normalverteilung zu gelangen.
Man wechselt zu einem Framework, welches nicht-normalverteilte Residuen modellieren kann. Dies sind dann vor allem generalized linear mixed-models (GLMMs), die u.a. mit Hilfe der Funktion
glmer()aus demlme4-Paket ausgewertet werden können (siehe auch Kap. 5 für weiterführende Informationen dazu).
In einer idealen Welt wäre Strategie 3 die beste Variante. GLMMs sind jedoch noch komplexer als ihre klassischen LMM-Gegenstücke und stellen nicht selten zusätzliche Anforderungen an ihre Schätzung und Diagnose. Aus diesem Grund ist es nicht ungewöhnlich, dass in der Praxis häufig Strategie 1 oder 2 anzutreffen ist.5 Für unsere Beispieldaten haben wir uns entschieden, die übliche log-Transformation anzuwenden. Ein Histogramm der Original-RTs und der log-transformierten RTs zeigt, dass diese dann eher einer Normalverteilung folgen als die untransformierten RTs (siehe Abb. 4.1).
par(mfrow = c(1, 2))
hist(f_simon_daten$rt,
breaks = 20,
main = "RTs",
xlab = "RTs"
) # vor der Transformation
f_simon_daten$logrt <- log(f_simon_daten$rt)
hist(f_simon_daten$logrt,
breaks = 20,
main = "log-RTs",
xlab = "log-RTs"
) # nach der Transformation
Abbildung 4.1: Vergleich der Verteilungen Original-RTs (links) und der log-transformierten RTs (rechts).
Die weitere Auswertung folgt i.W. derjenigen, die wir in Kapitel 2 bereits an dem Minimalbeispiel mit 3 Versuchspersonen eingeführt haben. Zunächst faktorisieren wir Versuchspersonen vpid und den interessanten fixed-effect cong (also Simon-Kongruenz) und stellen eine Effektkodierung sicher:
f_simon_daten$vpid <- as.factor(f_simon_daten$vpid)
f_simon_daten$cong <- as.factor(f_simon_daten$cong)
contrasts(f_simon_daten$cong) <- contr.sum(2)Dann berechnen wir das random slope Modell, d.h., wir lassen sowohl random intercepts als auch random slopes zu und schauen uns den Output an:
library(lme4)
# lme4:: aber sinnvoll, um explizit die Funktion aus dem lme4 Paket zu nutzen
lme_out <- lme4::lmer(logrt ~ 1 + cong + (1 + cong | vpid),
data = f_simon_daten
)
summary(lme_out)## Linear mixed model fit by REML ['lmerMod']
## Formula: logrt ~ 1 + cong + (1 + cong | vpid)
## Data: f_simon_daten
##
## REML criterion at convergence: -3659.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5825 -0.6303 -0.0645 0.5288 5.0228
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## vpid (Intercept) 0.0072606 0.08521
## cong1 0.0004784 0.02187 0.50
## Residual 0.0341452 0.18478
## Number of obs: 7113, groups: vpid, 32
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 6.239220 0.015222 409.88
## cong1 -0.017519 0.004446 -3.94
##
## Correlation of Fixed Effects:
## (Intr)
## cong1 0.430Zunächst lesen wir REML criterion at convergence: -3659.3. Diese Wert gibt das REML Kriterium an, mit welchem der Algorithmus die Parameterbestimmung beendet hat und er ist so etwas wie eine log-Likelihood. Unter Fixed effect finden wir die Koeffizienten (und geschätzten Standardfehler) der fixed-effects, die also quasi für alle Versuchspersonen gleichermaßen gelten (auch abrufbar mit fixef(lme_out). Unter Random effects finden wir die die Varianzen (bzw. Standardabweichungen) der Normalverteilungen (mit Erwartungswert \(0\)) aus denen die individuellen Abweichungen im Intercept und Slope stammen (die individuellen Werte können mit ranef(lme_out) abgerufen werden), sowie die Korrelation der individuellen random intercepts und random slopes. Der positive Wert hierbei bedeutet also z.B., wenn eine Versuchspersonen mehr nach oben im allgemeinen RT-Level abweicht, dann wäre ihr Slope auch größer, und damit der (Simon-)Effekt also größer. Diese drei Parameter sind die Komponenten der Kovarianzmatrix einer multivariaten Normalverteilung (mit Erwartungswertvektor \(\bf{0}\)). Unter Residual fnden wir letztlich noch die Schätzung der Fehlervarianz. Um nur diese (Ko-)Varianzparameter auszugeben, können wir auch VarCorr(lme_out) verwenden.
Am Ende finden wir noch Correlation of Fixed Effects. Dies mag auf den ersten Blick seltsam erscheinen, da die fixed-effects ja für alle Versuchspersonen gleichermaßen gelten. Allerdings sehen wir ja auch geschätzte Standardfehler der fixed-effects. Würden wir das LMM an immer wieder auf Basis anderer Zufallsstichproben schätzen, würden wir daher immer verschiedene Werte für die fixed-effects erhalten und deren Korrelation wird hier abgebildet. Eine ausführlichere Darstellung findet sich hier.
Was den bisherigen Outputs fehlt, sind inferenzstatistische Tests der fixed- und random-effects.
4.2 Inferenzstatistik bei LMMs
Im vorherigen Abschnitt haben wir ein LMM mit Hilfe der Funktion lmer() geschätzt. Hier haben wir auch gesehen, dass keinerlei \(p\)-Werte für die fixed-effects ausgegeben worden sind, obwohl diese in der Forschungspraxis natürlich oft benötigt werden. Wir beginnen mit den inferenzstatistischen Tests für die fixed-effects und wenden uns dann auch noch den Tests für random-effects zu.
4.2.1 Inferenzstatistische Tests der fixed-effects
Wie können wir nun \(p\)-Werte für die fixed-effects erhalten? Zunächst einmal muss man festhalten, dass die Berechnung von \(p\)-Werten im Zuge von LMMs alles andere als trivial ist und es gibt auch eine ganze Reihe von Ansätzen zu ihrer Bestimmung. Allerdings scheint uns auch nach wie vor nicht ganz klar zu sein, ob es eine grundsätzlich “beste” Variante gibt oder ob sich dies für verschiedene Situationen unterscheiden kann. Etwas vereinfacht gesagt ist das Hauptproblem, dass die Bestimmung der Freiheitsgrade unklar ist. Daher wird zwar im Output von lmer() ein \(t\)-Wert angeben, aber daraus kein \(p\)-Wert mehr berechnet.
4.2.1.1 Wald-Tests
Wald-Tests nutzen die geschätzten fixed-effect Koeffizienten und den zugehörigen geschätzten Standardfehler um eine Teststatistik herzuleiten. Unterschieden werden hierbei \(t\)-, \(F-\), \(\chi^2\) und \(z-\)Wald-Tests, welche aber alle sehr ähnlich sind. Beim Wald-\(z\)- und \(t\)-Test wird der Quotient aus Schätzer und Standardfehler berechnet, also \[t_0 = z_0 = \frac{\hat \beta_j}{\sqrt{\widehat{V(\hat \beta_j)}}}\,.\] Der Unterschied liegt lediglich darin, ob die Teststatistik anhand der \(t\)-Verteilung oder anhand der Standardnormalverteilung geprüft wird. Der Wald \(F\)-Test ist letztlich das Quadrat des \(t\)-Tests, d.h. es würde also \(F_0 = t_0^2\) anhand einer \(F\)-Verteilung getestet. Schließlich ist die Teststatistik beim Wald \(\chi^2\)-Test eine direkte Funktion des \(F\)-Werts und beim Test eines einzelnen Koeffizienten identisch. Getestet wird hierbei allerdings anhand einer \(\chi^2\)-Verteilung (mehr hierzu findet sich bei West et al. 2022, S. 37).
In R können wir Wald \(\chi^2\)-Tests mit Hilfe der Funktion Anova() aus dem Paket car anfordern. Hierbei sollte aber das Argument type = 3 definiert werden, damit der Signifikanztest nur den jeweiligen einzelnen fixed-effect Koeffizienten testet (also das Modell ohne den jeweiligen fixed-effect Koeffizienten wird gegenüber einem Modell mit allen fixed-effect Koeffizienten getestet und damit letztlich die Hypothese \(H_0: \beta_j = 0\)):
library(car)
Anova(lme_out,
type = 3
)Wald-Tests werden allerdings im Allgemeinen nicht empfohlen, da sie das nominelle \(\alpha\)-Niveau nicht halten können. Wie bereits erwähnt, liegt dies primär daran, dass die Freiheitsgrade der Tests nicht korrekt sind und andererseits, dass die geschätzten Standardfehler der fixed-effects, \(\sqrt{\widehat{V(\hat \beta_j)}}\), substanziell verzerrt seien können.
4.2.1.2 “\(t\) as \(z\)”-Analyse
Ein sehr einfacher Ansatz, der auch gerade in den frühe(re)n Jahren der LMM-Verwendung zum Tragen kam, basiert auf den Wald \(t\)-Werten und darauf, schlichtweg die Standardnormalverteilung zur Beurteilung der Signifikanz heranzuziehen. Die Logik dabei ist, dass die \(t\)- und die Standardnormalverteilung sich umso ähnlicher werden, je höher die Freiheitsgrade der \(t\)-Verteilung sind (wenn \(\text{df}\rightarrow \infty\) dann gehen beide ineinander über).
Eine implizit darauf basierende und manchmal verwendete Daumenregel ist es, fixed-effects als signifikant einzustufen wenn \(|t|>1.96\) ist (bzw. manchmal auch, wenn er größer als \(2\) ist). Insgesamt scheint diese Heuristik aber zu erhöhten Type 1 Fehlerraten zu führen (siehe dazu auch Luke, 2017).
4.2.1.3 Satterthwaite-Approximation und die Kenward-Roger Methode
Als Lösung zu den Problemen beim Wald-Test gibt es zwei Korrekturen, die helfen, korrekte \(p\)-Werte zu erhalten. Zum einen hat Sattertwaite (1946) eine Methode zur Korrektur der Freiheitsgrade vorgeschlagen, die tatsächlich auch viel Verwendung in der LMM Literatur findet. Darauf aufbauend haben Kenward und Roger (1997) zusätzlich eine Methode vorgeschlagen, bei der die Standardfehler nach oben korrigiert werden. Dies verbessert die Validität der \(p\)-Werte dann noch einmal.
Für typische Datensätze in der Psychologie wird diese Methode oft empfohlen (siehe bspw. Singmann & Kellen, 2019). In R können wir sie ebenfalls über die Anova() Funktion anfordern. Da die letztlich berechnete Statistik einem \(F\)-Wert entspricht, müssen wir das Argument test.statistic = "F" definieren:
Anova(lme_out,
type = 3,
test.statistic = "F"
)Dass hierbei die Freiheitsgrade approximiert werden, erkennt man gut an den entsprechenden nicht-ganzzahligen Werten. Auch erkennt man, dass die berechneten \(F\)-Teststatistiken leicht kleiner ausfallen als die vorherigen \(\chi^2\)-Werte (ohne Korrektur wären diese beim Test eines einzelnen Koeffizienten identisch).
Da die Methode nach Kenward-Roger je nach Komplexität des Modells recht rechenintensiv sein kann, kann ggf. auch nur auf die Korrektur der Freiheitsgrade nach Satterthwaite ausgewichen werden. Eine Implementation dieser Methode bietet das Paket lmerTest (Kuznetsova et al., 2017). Im Prinzip überschreibt dieses Paket die lmer()-Funktion und liefert dann direkt die \(p\)-Werte in der Ausgabe. Beachten Sie, dass hierbei wichtig ist, explizit mit der ::-Schreibweise zu arbeiten, um R zu erklären, welche lmer()-Funktion benutzt werden soll (man kann dann i.Ü. auf das Laden eines Paketes mit library() auch gänzlich verzichten):
library(lmerTest)
lme_out_SAT <- lmerTest::lmer(logrt ~ 1 + cong + (1 + cong | vpid),
data = f_simon_daten
)
summary(lme_out_SAT)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logrt ~ 1 + cong + (1 + cong | vpid)
## Data: f_simon_daten
##
## REML criterion at convergence: -3659.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5825 -0.6303 -0.0645 0.5288 5.0228
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## vpid (Intercept) 0.0072606 0.08521
## cong1 0.0004784 0.02187 0.50
## Residual 0.0341452 0.18478
## Number of obs: 7113, groups: vpid, 32
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 6.239220 0.015222 30.954650 409.88 < 2e-16 ***
## cong1 -0.017519 0.004446 30.437180 -3.94 0.000441 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cong1 0.4304.2.1.4 Likelihood-Ratio Tests
Ein Likelihood-Ratio Test (LRT) basiert auf der (skalierten) Differenz der Likelihoods (bzw. der Deviance-Statistik) zweier genesteter Modelle (siehe auch bereits hier im Kontext der logistischen Regression). Konkret wird hierbei ein Obermodell mit allen Parametern mit einem Untermodell, bei dem die interessierenden Parameter weggelassen werden, auf Basis der (maximierten) Likelihood verglichen. Der Name “Ratio” rührt dabei daher, dass ein Bruch der beiden Likelihoods betrachtet wird, der aber Rechenregeln folgend oft in eine Differenz überführt wird. Die resultierende Teststatistik ist dabei \(\chi^2\)-verteilt mit so vielen Freiheitsgraden wie Parameter im Vergleich vom Ober- zum Untermodell weggelassen wurden (beim Test eines einzelnen fixed-effects also 1): \[\begin{align*} -2\log\left( \frac{\mathcal{L_{\text{Untermodell}}}}{\mathcal{L_{\text{Obermodell}}}} \right) &= [-2\log(\mathcal{L}_{\text{Untermodell}})]-[-2\log(\mathcal{L}_{\text{Obermodell}})]\\ &=2\cdot[\log(\mathcal{L}_{\text{Obermodell}})-\log(\mathcal{L}_{\text{Untermodell}})]\sim\chi^2_{df = \Delta p} \end{align*}\]
LRTs können in R durch eine Übergabe von Ober- und Untermodell der Funktion anova() erhalten werden. Hierbei ist jetzt allerdings wichtig, dass die Modelle per ML und nicht per REML geschätzt wurden. Dies liegt eben daran, dass sich bei einer REML Schätzung auch die Schätzungen der Varianzkomponenten der random-effects unterscheiden und die Likelihoods so nicht mehr vergleichbar sind. Außerdem sollten die Modelle beim Test der fixed-effects sich trivialerweise auch nur in der fixed-effects Struktur unterscheiden und nicht zusätzlich in der random-effects Struktur (ansonsten testet man ja mehr als nur den einen intendierten fixed-effect). Da aber nun lmer() als Default mit einer REML Schätzung arbeitet, muss hier explizit eine ML Schätzung mit REML = FALSE erzwungen werden. Am Beispiel des einzigen fixed-factors cong sähe dies bspw. so aus:
# Obermodell
lme_out_ober <- lme4::lmer(logrt ~ 1 + cong + (1 + cong | vpid),
data = f_simon_daten,
REML = FALSE # kein REML, sondern ML-Schätzung
)
# Untermodell ohne cong als fixed effect
lme_out_unter <- lme4::lmer(logrt ~ 1 + (1 + cong | vpid),
data = f_simon_daten,
REML = FALSE # kein REML, sondern ML-Schätzung
)
# Likelihood-Ratio-Test
anova(
lme_out_ober,
lme_out_unter
)Tatsächlich würde die Funktion anova() selber ein “Refitting” mit ML vornehmen, wenn ihr zwei mit REML geschätzte Modelle übergeben werden.
Die manuelle Berechnung der Test-Statistik (Chisq im Output) können wir wie oben geschildert nachvoll-ziehen (siehe auch das entsprechende Kapitel in Statistik 2). Da diese Statistik im engeren Sinne kein \(\chi^2\)-Wert ist und nur für sehr große Stichproben approximativ \(\chi^2\)-verteilt ist, bezeichnen wir diese hier mit \(\lambda\):
# Teststatistik aus den log-Likelhoods berechnen
lambda <- 2 * (logLik(lme_out_ober) - logLik(lme_out_unter))
lambda <- as.numeric(lambda)
lambda # Teststatistik## [1] 13.00237# dann den p-Wert anhand einer chi^2-Verteilung bestimmen
1 - pchisq(lambda,
df = 1
)## [1] 0.0003110964Da die R-internen Funktionen zum Bau einer Designmatrix sich manchmal vehement “weigern” eine dummykodierte Variable aus dem Modell zu entfernen (und in eine effektkodierte Variable zu verwandeln), lohnt es sich zur Anwendung von LRTs manchmal, auf die Funktion mixed() aus dem afex Paket auszuweichen. Wir werden diese Funktion später noch kurz einführen. Außerdem werden LRTs zum Testen von fixed-effects nicht generell empfohlen, da deren approximative \(\chi^2\)-Verteilung nur für sehr viele Gruppen pro Gruppierungsvariable (bspw. also Versuchspersonen) und Beobachtungen pro Gruppe gilt (Pinheiro & Bates, 2000). Singmann und Kellen (2019) empfehlen daher LRT Tests erst ab einer Anzahl von ca. 40 bis 50 Gruppen pro Gruppierungsvariable.
4.2.1.5 Parametric Bootstrap
Bootstrapping haben wir in Statistik 1 und 2 bisher noch nicht besprochen. Es handelt sich hierbei um eine recht mächtige und weit verbreitete Technik, die die Rechenleistung moderner Computer ausnutzt. Um fixed-effects auf Signifikanz zu testen, ist die Idee hierbei, dass man zuerst das Untermodell (ohne den interessierenden fixed-effect) schätzt und die so erhaltenen Koeffizienten einer Simulation zugrunde legt, bei der wiederholt neue Daten simuliert werden. Diese simulierten Daten stellen also mögliche Datensätze dar, unter der Annahme, es gäbe keinen Effekt, d.h. der im Obermodell zusätzlich enthaltene Effekt sei nicht vorhanden. Im Anschluss wird sowohl das Ober- als auch Untermodell auf die simulierten Daten gefittet und so pro simuliertem Datensatz eine “synthethische” LRT-Statistik berechnet. Die Verteilung aller LRT-Statistiken kann als simulierte Verteilung unter der \(H_0\) angesehen werden und dient als Referenzverteilung. In anderen Worten “bastelt” man sich eine simulierte Verteilung der Prüfgröße auf Basis der vorhandenen Daten und nutzt keine mathematisch-hergeleitete Verteilung (wie bspw. die \(\chi^2\)-Verteilung mit \(x\) Freiheitsgraden). Der Name “Bootstrap” kommt daher, dass man sich auf Basis der eigenen Daten eine Verteilung generiert und damit in Analogie sich selbst am Schopfe aus dem Schlammassel zieht.6 Ein Aufruf ist direkt über die Funktion PBmodcomp() aus dem Paket pbkrtest möglich (dabei können Warnungen bzgl. der Modellkonvergenz entstehen, die wir hier in der Ausgabe unterdrückt haben; solange diese selten sind, stellt dies kein Problem dar):
library(pbkrtest)
set.seed(1)
PBmodcomp(
largeModel = lme_out_ober,
smallModel = lme_out_unter,
nsim = 100, # Anzahl der Simulationen; sollte > 1000 sein
seed = 1 # für die Reproduzierbarkeit
) Bootstrap-Methoden werden im Kontext von LMMs generell empfohlen, können aber je nach Komplexität des Modells sehr rechenintensiv sein und sind daher aus praktischer Sicht nicht immer möglich.
4.2.1.6 Zusammenfassung und Bezug zum Beispiel
Im allgemeinen sind die aktuellen “go-to” Optionen entweder die Methode nach Kenward-Roger oder das Parametric Bootstrapping. Sollten beide Verfahren zu rechenintensiv sein, wird oft die Methode der Freiheitsgradanpassung nach Satterthwaite empfohlen. Nur als “letzter Ausweg” sollten (naive) Wald-Tests verwendet werden. Für eine genauere Betrachtung der \(\alpha\)-Fehlerraten, siehe z.B. Luke (2017).
Da wir uns in den letzten Abschnitten nur mit den Tests an sich beschäftigt haben, hier nun die entsprechenden Tests der fixed-effects nach Kenward-Roger:
Anova(lme_out,
type = 3,
test.statistic = "F"
)Vor allem sehen wir hier also, dass der Simon-Effekt (dargestellt durch den fixed-factor cong) signifikant ist.
4.2.2 Inferenzstatistische Tests der random-effects
In den meisten Fällen ist man an den fixed-effects interessiert und die Berücksichtigung von random-effects ist “nur” ein Mittel zum Zweck, um Abhängigkeiten bzw. Gruppierungen in den Daten akkurat zu berücksichtigen. Dennoch gibt es Fälle in denen man gerne testen möchte, ob die Variabilität in den random-effects signifikant ist; bspw. um eben genau diese Variabilität als Ausdruck individueller Unterschiede zu testen (siehe z.B. Janczyk & Miller, 2024, für eine Anwendung). Wichtig ist dabei allerdings nochmal zu betonen, dass nicht die eigentlichen random-effects getestet werden, sondern die Varianzen und Kovarianzen der random-effects innerhalb der Matrix \(\bf {\mathcal D}\) (bzw. die jeweiligen Parameter innerhalb des Parametervektors \(\theta\)).
Außerdem sei betont, dass es meist sinnvoll ist möglichst alle random-effects im Modell zu lassen, gerade dann, wenn diese nicht der primäre Fokus der Analyse sind (außer das Modell konvergiert nicht oder liefert unzulässige Kovarianzmatrizen der random-effects). Denn die Annahme, dass es keine Variabilität in einem Effekt zwischen den Gruppen gibt, ist oft aus inhaltlichen Gründen nicht haltbar. Außerdem führt die Nicht-Berücksichtigung eines random-effects in einem Modell, wenn es aber in Wahrheit tatsächlich Variabilität bzgl. diesem gibt, zu einem deutlichen Anstieg von \(\alpha\)-Fehlern beim Test der entsprechenden fixed-effects (vgl. bspw. Barr et al., 2013). Vor diesem Hintergrund ist es auch kritisch “nicht-signifikante” random-effects auszuschließen, um Tests über die fixed-effects signifikant zu bekommen.
Es gibt zwei verbreitete Varianten, um random-effects zu testen, die beide auf LRTs basieren. Da hierbei die Varianzkomponenten als Parameter im Fokus stehen, wird i.d.R. empfohlen, hier nun die REML Methode zur Schätzung der Modelle zu nutzen (Fox & Weisberg, 2019; West et al., 2022).
4.2.2.1 Likelihood Ratio Tests für random-effects
Eine erste direkte Möglichkeit die Varianzen bzw. Kovarianzen von random-effects zu testen bieten LRTs anhand einer \(\chi^2\)-Verteilung. Das Vorgehen ist dabei ganz identisch zu den LRTs für fixed-effects: Ein Ober- und ein Untermodell werden geschätzt und die Funktion anova() führt den LRT (approximativ) anhand einer \(\chi^2\)-Verteilung durch.
Als Beispiel wollen wir hier testen, ob die Berücksichtigung der Varianzkomponente für die random slopes das Modell signifikant verbessert. Dies würde darauf hindeuten, dass nicht vernachlässigbare interindividuelle Unterschiede tatsächlich vorliegen. Hierfür müssen wir ein Untermodell spezifizieren, bei dem die Varianz nicht zugelassen wird, d.h. wir schätzen als Untermodell das random intercept Modell:
# Schätzung der beiden Modelle, wobei das Untermodell keinen random slope
# mehr enhält
lme_out_ober <- lme4::lmer(logrt ~ 1 + cong + (1 + cong | vpid),
data = f_simon_daten
)
lme_out_unter <- lme4::lmer(logrt ~ 1 + cong + (1 | vpid),
data = f_simon_daten
)
VarCorr(lme_out_ober) ## Groups Name Std.Dev. Corr
## vpid (Intercept) 0.085209
## cong1 0.021872 0.500
## Residual 0.184784VarCorr(lme_out_unter) # nur noch Varianz des intercepts vorhanden## Groups Name Std.Dev.
## vpid (Intercept) 0.085449
## Residual 0.186001Nun können wir den LRT durchführen. Allerdings führt die Funktion anova() “sicherheitshalber” einen Refit mit der ML Methode durch. Da wir dies an dieser Stelle aber explizit verhindern wollen, müssen wir dies mit refit = FALSE unterdrücken:
# LRT-Test
anova(lme_out_unter,
lme_out_ober,
refit = FALSE # nicht neu mit ML fitten, sondern REML beibehalten
)Das Ergebnis zeigt, dass sich der Modellfit zwischen beiden Modellen also unterscheidet, der random slope also substanziell zur Varianzaufklärung beiträgt. Ein wichtiger Punkt zum Verständnis hier ist noch: Der Test hat offenbar \(2\) Freiheitsgrade. Dies liegt daran, dass man sich bewusst sein muss, dass tatsächlich auch zwei Komponenten entfernt wurden: zum einen die Varianz der random slope Verteilung, aber als Konsequenz dessen auch die Kovarianz zwischen random slope und random intercept.
Tatsächlich müssten wir außerdem eigentlich zwei Fälle unterscheiden, nämlich den wenn Parameter unter der Nullhypothese an der Grenze ihrer möglichen Ausprägungen liegen (man sagt dann, sie lägen auf dem Rand des Parameterraums) und wenn sie dies nicht tun. Testen wir z.B. auf einen Wert von \(0\), dann gilt der erste Fall für Varianzen während der zweite Fall für Kovarianzen/Korrelationen zutrifft.
Der direkte Vergleich des Ober- und des Untermodells mit einem LRT ist daher gewissermaßen eine Vereinfachung dieses Vorgehens, die allerdings konservativer ist und daher, im Fall eines signifikanten Unterschieds im Modellfit, nicht problematisch ist (und in der Praxis eher üblich scheint). Wir haben diese beiden Fälle hier im Detail behandelt und weitere Informationen gibt es auch bei West et al. (2022) und bei Fox und Weisberg (2019).
4.2.2.2 Parametric Bootstrapping
Analog zum Test derfixed-effects, kann auch zum Test der random-effects auf Parametric Bootstrapping zurückgegriffen werden, wobei der Ansatz grundsätzlich identisch ist. Die Umsetzung in R ist leider nicht ganz ideal, denn es gibt die Funktion PBmodecomp(), allerdings erzwingt diese eine ML-Schätzung, keine REML-Schätzung. Ein beispielhafter Aufruf um zu testen, ob der random slope im Beispieldatensatz weggelassen werden kann, würde wie folgt aussehen (die Modelle wurden im vorherigen Code-Chunk geschätzt):
PBmodcomp(
largeModel = lme_out_ober, # Obermodell
smallModel = lme_out_unter, # Untermodell
seed = 2,
nsim = 100 # auch hier nsim höher stellen
) Eine leicht nutzbare Funktion für einen parametric bootstrap mit Hilfe der REML-Schätzung ist uns nicht bekannt. Insofern man diese wirklich benötigt, muss selbst eine entsprechende Funktion auf Basis von simulate() und refit() geschrieben werden.
4.3 Auswertung mit der Funktion mixed() aus dem Paket afex()
Das Paket afex (Singmann et al., 2024; siehe auch Singmann & Kellen, 2019) stellt mit der Funktion mixed() eine Wrapper-Funktion zur Arbeit mit LMMs zur Verfügung, die den Aufruf und die korrekte Berechnung von p-Werten erleichtern soll (siehe auch hier für ein Beispiel). Der Aufruf erfordert zwar einerseits die korrekte Spezifizierung der Modellformel, andererseis kann die Art der p-Wert-Berechnung mit dem Argument method = "..." direkt angegeben werden. Hierbei stehen S (default), KR, LRT und PB für die Satterthwaite-Methode, die Kenward-Roger Korrektur, LRTs und Parametric Bootstrapping. Wollen wir in unserem Beispiel die Signifikanz des Faktor Kongruenz mit einem LRT evaluieren, würde der Aufruf wie folgt aussehen (das notwendige Untermodell erzeugt mixed() dann intern von selbst):
library(afex)
mixed(logrt ~ 1 + cong + (1 + cong | vpid),
data = f_simon_daten,
method = "LRT"
)## Mixed Model Anova Table (Type 3 tests, LRT-method)
##
## Model: logrt ~ 1 + cong + (1 + cong | vpid)
## Data: f_simon_daten
## Df full model: 6
## Effect df Chisq p.value
## 1 cong 1 13.00 *** <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1Zwei weitere Aspekte von mixed() wollen wir hier kurz ansprechen:
- Die Funktion
mixed()stellt sicher, dass die fixed-factors effektkodiert sind. In experimentellen ANOVA-ähnlichen Designs ist dies tatsächlich wichtig, spätestens dann, wenn Interaktionen betrachtet werden (daher haben wir die Kodierung hier auch jeweils direkt verwendet). Es gibt aber Situationen, in denen eine Dummykodierung sinnvoller wäre. Dann ist oft ein Rückgriff auf dielmer()Funktion sinnvoller. - Bei der Verwendung mancher fortgeschrittener Syntax kommt es bei
lmer()manchmal zu Problemen, wenn die Stufen der Faktoren semantische Bezeichnungen haben. Die Verwendung der Optionexpand_re = TRUEinmixed()sorgt dafür, dass diese Probleme nicht auftreten.
Insgesamt ist afex ein vielbenutztes Paket und wir greifen bei den Beispielen hier auch auf die Funktion mixed() weiter zurück.
4.4 Zusammenfassung und Ausblick
Wir haben hier nur eine einfache Version der Auswertung von Daten mit LMMs erläutert und Wert auf die verschiedenen inferenzstatistischen Verfahren gelegt. Im Kontext dieser Einführung soll dies reichen. Weitere Ausführungen bei komplexeren Modellen und anderen Fragestellungen (inkl. Fragen der Modelldiagnostik, … ) sind hier zu finden. Dort besprechen wir auch den (nicht seltenen) Fall von Konvergenzproblemen und dem möglichen Umgang damit.
5 Ausblick
Wir haben in den letzten Kapiteln nur eine Einführung in LMMs gegeben. In diesem Kapitel geben wir schließlich noch kurze Ausblicke auf weitere relevante Themen.
5.1 Die gleichzeitige Modellierung mehrerer random-factors (d.h. mehrere Gruppierungsvariablen)
Die bisherigen Betrachtungen und Gleichungen zum klassischen LMM haben immer nur einen Gruppierungsfaktor bzw. nur einen Term \({\bf Z_i} b_i\) beinhaltet und werden daher auch als “single-level” Modell bezeichnet (Pinheiro & Bates, 2000). Es ist aber grundsätzlich möglich, das Modell zu erweitern bzw. leicht abzuändern, um weitere Gruppierungsvariablen und deren Beziehung zueinander zu modellieren. Eine hierbei häufig gemachte Unterscheidung sind sog. crossed vs. nested random-effects. Diese zwei Fälle sind in Abbildung 5.1 illustriert (in Anlehnung an Singmann & Kellen, 2019) und werden im Folgenden kurz erläutert.
Abbildung 5.1: Illustration von crossed und nested random-effects.
5.1.1 Crossed random-effects
Bisher haben wir nur über Versuchspersonen hinweg generalisieren wollen und daher mit Versuchspersonen auch nur einen random-factor betrachtet und modelliert. In den entsprechenden Funtionsaufrufen von lmer() gab es auch nur einen random-effect Term wie etwa ... + (1 | subject). Eine weitere Stärke von LMMs ist es allerdings, dass auch mehrere random-factors simultan berücksichtigt werden können. Wenn Beobachtungen sowohl unter dem einen als auch unter dem anderen random-factor (Gruppierungsfaktor) zusammengefasst werden können, dann sprechen wir von crossed random-effects.
Um zu verdeutlichen, wo dies typischerweise in der Psychologie vorkommt, greifen wir zunächst auf ein empirisches Beispiel zurück, bei dem mit sprachlichen Stimuli gearbeitet wurde.7 In der Sprachpsychologie gibt es bspw. den Effekt, dass es einen Zusammenhang zwischen (deiktischer) Zeit und Räumlichkeit gibt. In einem typischen Experiment sollen Versuchspersonen dann auf einen sprachlichen Stimulus mit einer linken oder einer rechten Reaktion angeben, ob dieser Stimulus sich auf die Vergangenheit (“damals”, “gestern”, …) oder auf die Zukunft (“bald”, “übermorgen”, …) bezieht. In einer kongruenten Bedingung soll auf vergangenheitsbezogene Stimuli mit links und auf zukunftsbezogene Stimuli mit rechts reagiert werden; in einer inkongruenten Bedingung umgekehrt, also auf vergangenheitsbezogene Stimuli mit rechts und auf zukunftsbezogene Stimuli mit links. Typischerweise sind in solchen Experimenten die RTs in den kongruenten Bedingungen kürzer als in den inkongruenten Bedingungen (z.B. Janczyk & Ulrich, 2019; für eine Metaanalyse, siehe von Sobbe et al., 2019).
Hier gibt es nun zwei random-effects über die generalisiert werden soll:
- Zum einen gibt es wieder die Versuchsperson, die aus einer Population von Menschen stammen und die Befunde sollen für diese Population gelten. Dabei unterscheiden sich die Menschen in ihrem allgemeinen RT-Level (random intercept), aber auch in der Größe des Kongruenzeffektes (random slope).
- Zum anderen, und dies ist neu, sind auch die Stimuli als Stichprobe aus einer größeren Population möglicher sprachlicher Ausdrücke für die Experimente aufzufassen. Das Ergebnis der Studie soll aber nicht nur für diese speziellen Stimuli gelten, sondern für alle potenziell möglichen Stimuli. Dabei können, genau wie bei den Versuchspersonen auch, manche der Wörter generell zu eher schnellen, andere zu langsamen Reaktionen führen (random intercept) und, da jeder Stimulus auch in der kongruenten und der inkongruenten Bedingung verwendet wird, kann auch ein Kongruenzeffekt für die Stimuli berechnet werden, der für die verwendeten Wörter unterschiedlich groß ausfallen kann (random slope).
Ähnliche Situationen kommen auch in anderen Kontexten vor, z.B. wenn Versuchspersonen Emotionen von Gesichtern beurteilen sollen und die Gesichter als Stimuli nur eine Stichprobe möglicher Gesichter sind oder wenn in der Entscheidungspsychologie Entscheidungen auf Basis von Hinweisreizen (“Cues”) getroffen werden sollen und diese nur eine Stichprobe möglicher Szenarien beschreiben. In beiden Fällen sollen aber die Schlussfolgerungen auch generalisiert werden über die Stimuli, also die Gesichter bzw. die Szenarien. Typischerweise redet man daher oft von Subjects und Items als random- factors.
Wir illustrieren dies und die entsprechende Auswertung hier nun einmal mit simulierten Daten. Dies hat für uns hier mehrere Vorteile:
- Wir knüpfen an den Teil 1 an und üben das grundsätzliche Vorgehen einer Simulation noch einmal. Gleichzeitig erlaubt uns dieses Vorgehen einen eher praktischeren Blick auf die Zerlegung der Varianz einer abhängigen Variablen durch LMMs.
- Da wir die Daten simuliert haben und sie quasi ideal sind, kennen wir die ground truth. Wir können daher leicht erkennen, wo sich die entsprechenden Werte im Output von
lmer()finden lassen (vgl. dazu auch DeBruine & Barr, 2021).
Wir bauen die Simulation hier schrittweise auf und es empfiehlt sich, die einzelnen Schritte einmal nachzuvollziehen, um das Geschehen besser verstehen zu können. Ziel ist es, Daten für das folgende Modell zu generieren: \[Y_{ijk} = \beta_0 + \beta_1 \cdot x_{ijk} + b_{0,j} + b_{1,j} \cdot x_{ijk} + b_{0,k} + b_{1,k} \cdot x_{ijk}+\varepsilon_{ijjk}\] Hier ist \(Y_{ijk}\) eine Beobachtung \(i\) für Item \(j\) und Person \(k\). Die fixed-effects sind \(\beta_0\) und \(\beta_1\). Der Wert \(x_{ijk}\) ist ein Prädiktorwert (im Folgenden einer nominalskalierten Variable). \(b_{0,j}\) und \(b_{1,j}\) sind ein random intercept und slope für ein Item \(j\). Entsprechend sind \(b_{0,k}\) und \(b_{1,k}\) random intercepts und slopes für eine Person \(k\).
Wir beginnen mit einigen Parametern für die Simulation im Allgemeinen und der Festlegung der Werte für die Modellparameter:
set.seed(1) # um die Simulation reproduzierbar zu machen
n_subjects <- 50 # Wieviele Versuchspersonen?
n_items <- 50 # Wieviele Items?
conditions <- c(-0.5, 0.5) # Wir benennen die beiden Bedingungen mit -0.5 und 0.5 (x_ijk)
# damit der fixed effect nachher den generellen RT Unterschied
# zwischen beiden Bedingungen darstellt. Die beiden Bedingungen
# könnten auch, wie im Beispiel, kongruente und inkongruente
# Zuordnungen von Stimuli zu Tasten sein
sigma_residual <- 20 # SD der Residuen
# fixed effects
beta_0 <- 800 # Intercept
beta_1 <- 100 # Slope
# Varianzkomponenten subject random factor
SD_0s <- 30 # Intercept
SD_1s <- 200 # Slope
COV_01s <- 4000 # Kovarianz Intercept-Slope
# Varianzkomponenten item random factor
SD_0i <- 20 # Intercept
SD_1i <- 40 # Slope
COV_01i <- 200 # Kovarianz Intercept-SlopeFür die Simulation brauchen wir gleich die Kovarianzmatrizen für die beiden random-factors, die wir daher als nächstes zusammenstellen:
# Kovarianzmatrix für Subject random-factor
CovMatrixSubject <- matrix(
c(
SD_0s^2, COV_01s,
COV_01s, SD_1s^2
),
nrow = 2, ncol = 2,
byrow = TRUE
)
# Kovarianzmatrix für Item random-factor
CovMatrixItems <- matrix(
c(
SD_0i^2, COV_01i,
COV_01i, SD_1i^2
),
nrow = 2, ncol = 2,
byrow = TRUE
)Nun können wir für die Subjects und die Items die notwendigen Werte für Intercept und Slope simulieren, die wir benötigen, um die entsprechenden RTs gemäß dem zugrunde liegenden Modell zu berechnen:
# Wir beginnen mit Intercept und Slope für den Subject random-factor
# simulieren der Daten gemäß der oben spezifizirten Kovarianzmatrix und mit einem
# Erwartungswert 0 -> egibt alle b_0k und b_1k
randomEffects_Subjects <- MASS::mvrnorm(
n = n_subjects,
mu = c(0, 0),
Sigma = CovMatrixSubject,
empirical = TRUE
)
# Die simulierten Werte und die Labels für die Versuchspersonen in ein DataFrame...
randomEffects_Subjects <- data.frame(
Subject = c(1:n_subjects),
round(randomEffects_Subjects)
)
# ...und passende Namen geben
names(randomEffects_Subjects)[2:3] <- c("Intercept_Subject", "Slope_Subject")
# Das gleiche machen wir dann für den Item random-factor
# und erhalten so alle b_0j und b_1j
randomEffects_Items <- MASS::mvrnorm(
n = n_items,
mu = c(0, 0),
Sigma = CovMatrixItems,
empirical = TRUE
)
randomEffects_Items <- data.frame(
Item = c(1:n_items),
round(randomEffects_Items)
)
names(randomEffects_Items)[2:3] <- c("Intercept_Item", "Slope_Item")Nun beginnen wir mit der Zusammenstellung des eigentlichen DataFrames, der am Ende die simulierten RTs beinhaltet und den wir danach mit einem passenden LMM auswerten wollen:
# Zunächst legen wir den DataFrame und entsprechende Variablen für Subject,
# Condition und Item an.
data <- data.frame(
Subject = rep(c(1:n_subjects), each = 2 * n_items),
Condition = c(-0.5, 0.5),
Item = rep(rep(c(1:n_items), each = 2), times = n_subjects)
)
# Dann fügen wir die Werte für die fixed-effects hinzu, die ja für alle Versuchspersonen
# und Item gelten sollen.
data$beta_0 <- beta_0
data$beta_1 <- beta_1
# Nun fügen wir dem DataFrame die Intercepts und Slopes des Item random-factors hinzu...
data <- merge(data, randomEffects_Items,
by = "Item"
)
# ...und danach das Ganze noch für den Subject random-factor.
data <- merge(data, randomEffects_Subjects,
by = "Subject"
)
# Als letzter Schritt wird eine Variable mit den Residuen hinzugefügt.
data$residual <- round(rnorm(
n = 2 * n_subjects * n_items,
mean = 0,
sd = sigma_residual
))Nun können wir für jede einzelne Beobachtung die RT gemäß dem zugrunde liegenden LMM berechnen (siehe auch nochmal die Gleichung oben, vor dem ersten R-Code):
data$RT <- data$beta_0 + data$Condition * data$beta_1 +
data$Intercept_Item + data$Condition * data$Slope_Item +
data$Intercept_Subject + data$Condition * data$Slope_Subject +
data$residualInsgesamt sieht der Datensatz nun wie folgt aus (wegen des Umbruchs hier leider nicht so schön formatiert):
car::some(data) Nun stellen wir in Abbildung 5.2 als einen Check einmal zwei Darstellungen nebeneinander: Einmal aggegieren wir über die Items und erhalten die RTs für die Bedingungen pro Versuchspersonen (links) und einmal aggregieren wir über die Subjects und erhalten die RTs für die Bedingungen pro Item (rechts; in den Abbildungen dann jeweils nochmal über Subjects oder Items gemittelt). Ersichtlich ist, dass die “kongruente” Bedingung (also Condition = -0.5) zu kürzeren RTs führt als die inkongruente Bedingung, und dies sowohl für Subjects als auch für Items:
# wir kopieren uns den Datensatz in ein neues DataFrame, damit lmer() nachher ohne
# faktorisierte Variablen arbeiten kann
dataPlot <- data
dataPlot$Subject <- as.factor(dataPlot$Subject)
dataPlot$Item <- as.factor(dataPlot$Item)
dataPlot$Condition <- as.factor(dataPlot$Condition)
# beide Plots werden mit ezPlot() erstellt und gespeichert
plot1 <- ez::ezPlot(
data = dataPlot,
dv = RT,
wid = Subject,
within = Condition,
x = Condition
)
plot2 <- ez::ezPlot(
data = dataPlot,
dv = RT,
wid = Item,
within = Condition,
x = Condition
)
# ...und dann nebeneinander dargestellt und mit einem Titel versehen.
gridExtra::grid.arrange(plot1 + ggplot2::ggtitle("Pro Subject"),
plot2 + ggplot2::ggtitle("Pro Item"),
ncol = 2
)
Abbildung 5.2: Illustration des Kongruenzeffektes für die Subjects (links) und die Items (rechts).
rm(dataPlot)Nun nutzen wir lmer(), um ein entsprechendes LMM mit zwei random-factors zu schätzen. Die Formulierung ist dabei ganz ähnlich zu dem, was wir bisher genutzt haben, nur dass es eben zwei Terme mit random-effects gibt. Da wir einen \(p\)-Wert für den fixed-effect Condition erhalten möchten, benutzen wir hier das entsprechende Paket lmerTest, welches die Satterthwaite Approximation der Freiheitsgrade durchführt:
M1 <- lmerTest::lmer(RT ~ Condition + (Condition | Subject) + (Condition | Item),
data = data
)
summary(M1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RT ~ Condition + (Condition | Subject) + (Condition | Item)
## Data: data
##
## REML criterion at convergence: 45458.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4572 -0.6370 -0.0087 0.6635 3.7238
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 915.2 30.25
## Condition 39991.2 199.98 0.66
## Item (Intercept) 409.9 20.25
## Condition 1661.7 40.76 0.24
## Residual 420.4 20.50
## Number of obs: 5000, groups: Subject, 50; Item, 50
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 799.849 5.156 85.135 155.123 < 2e-16 ***
## Condition 98.755 28.869 53.040 3.421 0.00121 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Condition 0.565Wir können die einzelnen Ergebnisse nun einmal der Reihe nach bewerten. Dazu extrahieren wir diese aus dem Output und vergleichen Sie mit den Werten, die der Simulation zugrunde liegen. Wir beginnen mit den fixed-effects und stellen fest, dass diese sehr gut mit den Vorgaben der Simulation übereinstimmen:
# Vorgaben fixed effects
cat(paste0("Intercept: ", beta_0, "; Slope/Condition: ", beta_1))## Intercept: 800; Slope/Condition: 100# Schätzungen
fixef(M1)## (Intercept) Condition
## 799.8492 98.7552Also nächstes schauen wir uns die Varianzkomponenten an. Dazu extrahieren wir mit VarCorr() die Varianzkomponenten:
# Vorgaben Varianzkomponenten Subject
cat(paste0(
"SD Intercept: ", SD_0s, "; SD Slope/Condition: ", SD_1s,
"; Korrelation: ", round(COV_01s / (SD_0s * SD_1s), 2)
))## SD Intercept: 30; SD Slope/Condition: 200; Korrelation: 0.67# Vorgaben Varianzkomponenten Item
cat(paste0(
"SD Intercept: ", SD_0i, "; SD Slope/Condition: ", SD_1i,
"; Korrelation: ", round(COV_01i / (SD_0i * SD_1i), 2)
))## SD Intercept: 20; SD Slope/Condition: 40; Korrelation: 0.25# Schätzungen
VarCorr(M1)## Groups Name Std.Dev. Corr
## Subject (Intercept) 30.252
## Condition 199.978 0.663
## Item (Intercept) 20.247
## Condition 40.764 0.235
## Residual 20.504Schließlich können wir auch die Standardabweichung der Residuen extrahieren (oder sie der Ausgabe der Varianzkomponenten entnehmen) und mit dem vorgegebenen Wert vergleichen:
# Vorgabe
cat(paste0("SD Residuen: ", sigma_residual))## SD Residuen: 20# Schätzung
round(sigma(M1), 2)## [1] 20.5Es bleibt noch die Korrelation der fixed-effects. Wie vorher bereits erwähnt, ist dies die Korrelation zwischen den beiden fixed-effects, die entstehen würde, wenn wir wiederholt eine Zufallsstichprobe ziehen und dann das gleiche LMM schätzen würden. Dies ist daher auch kein Spezifikum eines LMMs, hat aber oft auch keine direkt wichtige Implikation für die Interpretation. Wir haben hier die obige Simulation ergänzt, um dieses Konzept zu illustrieren.
5.1.2 Nested random-effects
Nested random-effects treten dann auf, wenn Gruppen aus einem Gruppierungsfaktor wieder zu einer Gruppe eines zweiten Gruppierungsfaktors auf einem höheren Level zusammengefasst werden können. Ein klassisches Beispiel hierfür wäre eine Studie, bei der Schüler:innen als Beobachtungen betrachtet werden (Level 1). Diese gehören wiederum zufällig gewählten Klassen an (Level 2), wobei jede Klasse wiederum Teil einer zufällig gezogenen Schule (Level 3) ist. Da eine Klasse aber eindeutig (nur) zu einer Schule gehört, sind die Gruppierungen “geschachelt”, wobei der übliche Fachausdruck dafür genested ist. Während die abhängige Variable (das Kriterium) stets auf Level 1 gemessen wird, können die Variablen, deren Wirkung auf die abhängige Variable untersucht werden sollen, auf allen Levels definiert sein. So kann z.B. für die Klasse deren Größe erfasst werden und für die Schule der Schultyp.
Oft wird dann eine ganz ähnliche Syntax wie im vorangegangenen Beispiel verwendet. In Kapitel 4 von West et al. (2022) wird ein typisches Schüler:innen-Klassen-Schulen-Beispiel besprochen und eine Übersicht findet sich in Tabelle 15.2 bei Galecki und Burzykowski (2013).
5.2 Effektstärken und Power-Analyse
Die Frage nach standardisierten Effektstärken ist eine relativ ungeklärte Thematik. Es gibt, sehr vereinfacht, für sehr spezielle Situationen wie z.B. einen fixed-effect mit zwei Abstufungen Vorschläge, Effektstärken in der Art von Cohen’s \(d\) zu berechnen (siehe Westfall et al., 2014). Eine generelle Empfehlung gibt es allerdings nicht. Stattdessen werden oft die (unstandardisierten) Regressionskoeffizienten als Effekte berichtet (bzw. z.B. Odds-Ration, wenn die abhängige Variable dichotom ist; siehe weiter unten im Kontext GLMMs).
Power-Analysen haben wir konzeptuell bereits im Zuge des t-Tests und der ANOVA in Statistik 1 und 2 eingeführt. Ziel der Power-Analyse ist es zu evaluieren, wie häufig ein Test signifikant werden würde, unter der Annahme, ein bestimmter Effekt (bspw. ein fixed- oder random-effect) von einer bestimmten Größe läge vor. Die häufigste Power-Analyse ist die sog. a-priori Power-Analyse. Hier plant man die Stichprobengröße für ein intendiertes Analyse-Design derart, dass ein interessierender Effekt mit einer Wahrscheinlichkeit von z.B. \(1-\beta=0.8\) signifikant wird, unter der Annahme, dass dieser Effekt tatsächlich vorliegt. Die Power-Analyse ist also ein wichtiges Instrument in der Stichprobenplanung und essenziell für die Beurteilung, ob eine Studie überhaupt in der Lage ist einen vermuteten Effekt aufzudecken oder nicht.
Wie auch im Fall anderer Power-Analysen, wird bei LMMs zwischen simulativen und analytischen Ansätzen unterschieden. Bei analytischen Ansätzen wird die Power auf Basis geschlossener Gleichungen hergeleitet und ist damit konzeptuell ähnlich zur Power-Analyse für den \(t\)-Test bzw. einfache Fälle der ANOVA. Weitere Informationen hierzu findet sich z.B. bei Verbeke und Molenberghs (2000) und eine Variante ihrer Umsetzung in R wird in Galecki und Burzykowski (2013) gezeigt. Auch das Paket powerlmm (welches inzwischen allerdings nicht mehr auf CRAN zu finden ist), liefert einen Zugang zu diesen analytischen Lösungen (siehe auch Westfall et al., 2014). Allerdings sind diese analytischen Zugänge häufig auf bestimmte Designs beschränkt. In den nicht-abgedeckten Fällen wird dann auf einen simulativen Ansatz zurückgegriffen. Die grundsätzliche Logik lässt sich wie folgt zusammenfassen:
Modell festlegen: Das LMM wird mit bestimmten Parametern, die nachgewiesen werden sollen, mit plausiblen Werten exakt spezifiziert. Beispielsweise wird in einem Modell mit einem Prädiktor der fixed intercept, der fixed slope, die Varianz der entsprechenden random-effects sowie die Residualvarianz definiert.
Zufällige Beobachtungen generieren: Eine festgelegte Anzahl von Beobachtungen wird dann basierend auf dem Modell zufällig erzeugt. Zum Beispiel könnten 30 Versuchspersonen mit jeweils 20 Beobachtungen simuliert werden.
Modell anpassen und testen: Das Modell wird dann auf den generierten Datensatz gefittet und die intendierten Signifikanztests werden durchgeführt.
Die Schritte 1 bis 3 werden mit unterschiedlicher Anzahl von Versuchspersonen und Beobachtungen pro Versuchsperson wiederholt. Die Power wird dann über die relative Häufigkeit signifikanter Effekte geschätzt.
In R bieten die Pakte powerlmm (Magnusson, 2024), mixedpower (Kumle et al., 2021; auch nicht mehr auf CRAN verfügbar) und simr (Green & MacLeod, 2016) einfachere Zugänge zu simulativen Ansätzen und eine anwendungsfreundliche Einführung zu deren Umsetzung findet sich in Kumle et al. (2021).
Es sei noch erwähnt, dass die Power-Analyse für LMMs meist deutlich mehr Überlegungen benötigt als eine Power-Analyse im Zuge von t-Tests oder ANOVAs. Der Grund hierfür ist, dass mehrere Komponenten des LMMs Einfluss auf die Power haben (bspw. die Anzahl der Gruppen, Beobachtungen aber auch die Varianzkomponenten der random-effects) und es zudem bisher keine Einigung über gängige Effektstärken gibt. Insofern man nicht detaillierte Informationen über die intendierte Analyse auf Basis einer Vorstudie besitzt, muss man durch eigene Überlegungen realistische Annahmen über die einzelnen Komponenten des Modells machen, was naheliegender Weise nicht immer ganz trivial ist.
5.3 Generalized (Linear) Mixed-Models (GLMMs)
Genau wie die Verteilungsannahmen des ALMs aufgelockert werden können, wenn zum “Verallgemeinerten Linearen Modell” (GLM für Generalized Linear Model) übergegangen wird (siehe hier für eine kurze Betrachtung in Statistik 2), gibt es Generalized (Linear) Mixed-Models (GLMMs), mit deren Hilfe nicht-intervallskalierte und nicht-normalverteilte abhängige Variablen analysiert werden können. Ein typisches Beispiel sind auch hierbei (binomalverteilte) dichotome Variablen (wie z.B. Fehler in einem typischen RT-Experiment).
Während grundlegend nun LMMs mit der entsprechenden Variante im GLM zusammengebracht werden müssen, gibt es dennoch vieles zu beachten, wenn mit GLMMs gearbeitet wird (siehe z.B. Bolker et al., 2008, oder Stroup, 2013). Das lme4-Paket bietet z.B. die Funktion glmer() zur Analyse an.
Eine für unsere Bereiche typische Situation ist die Analyse von Fehler in Experimenten. Um dies kurz zu demonstrieren, nutzen wir den Simon-Datensatz aus Kapitel 4, allerdings benutzen wir nun alle Trials, unabhängig davon ob sie korrekt (0) oder fehlerhaft (1) waren (Variable error):
f_simon_daten <- subset(
simon_daten,
blocknr > 0 & errortype == 0 &
rt >= 200 & rt <= 1200
)
f_simon_daten$vpid <- as.factor(f_simon_daten$vpid)
f_simon_daten$cong <- as.factor(f_simon_daten$cong)
contrasts(f_simon_daten$cong) <- contr.sum(2)Nun können wir die Funktion glmer() nutzen. Während die Formel und die random-effect Struktur quasi identisch ist zum Aufruf weiter oben, müssen wir nun noch die Verteilungsannahmen und die Link-Funktion spezifieren (siehe auch logistische Regression):
glme_out <- lme4::glmer(error ~ cong + (1 + cong | vpid),
data = f_simon_daten,
family = binomial(link = "logit")
)
summary(glme_out)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: error ~ cong + (1 + cong | vpid)
## Data: f_simon_daten
##
## AIC BIC logLik -2*log(L) df.resid
## 3306.9 3341.5 -1648.4 3296.9 7565
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.5970 -0.2682 -0.2117 -0.1682 6.9071
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## vpid (Intercept) 0.3447 0.5871
## cong1 0.1335 0.3654 -0.04
## Number of obs: 7570, groups: vpid, 32
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.97244 0.11931 -24.913 <2e-16 ***
## cong1 -0.19980 0.08734 -2.288 0.0222 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cong1 0.042Der Ausgabe von summary() können wir viele Komponenten entnehmen, die wir nun bereits kennen. Neu ist, dass hier direkt ein Signifikanztest (in Form eines Wald-Tests) ausgegeben wird. In diesem Fall spräche dessen Ergebnis für einen signifikanten Kongruenzeffekt auch auf die Fehlerrate. Im Kontext von GLMMs sind Wald-Tests allerdings mit Vorsicht zu genießen und häufig wird auf entsprechende LRTs zurückgegriffen, die nach dem gleichen Schema wie bei LMMs durchgeführt werden. Wenn wir glme_out bereits als das Obermodell auffassen, wäre eine Durchführung wie folgt:
# Untermodell ohne fixed-effect schätzen
glme_out_unter <- lme4::glmer(error ~ 1 + (1 + cong | vpid),
data = f_simon_daten,
family = binomial(link = "logit")
)
# LRT durchführen
anova(glme_out, glme_out_unter)In diesem Fall kommen beide Wege zur gleichen Entscheidung.
5.4 Schlussanmerkungen
Zur Vorsicht wollen wir darauf hinweisen, dass es im Gebiet der LMMs nach wie vor viel Unklarheit über den Umgang mit speziellen Situationen und speziellen Problemen gibt. Dies rührt auch daher, dass viele der intensiv mit diesen Modellen befassten Personen aus dem Kreis der mathematischen Statistik stammen, während die Nutzung dieser Modelle z.B. in der Psychologie noch in den Anfängen steckt. Da es zudem viele “handwerkliche Freiheitsgrade” gibt (z.B. wie die random-effect Struktur reduziert werden kann, wie \(p\)-Werte bestimmt werden, …), halten wir für es essentiell, in entsprechenden Arbeiten, Publikationen genau Auskunft über das Vorgehen zu geben, um die Ergebnisse (und ggf. auch die Probleme damit) als Leser:in besser einschätzen zu können.
6 Referenzen
Baayen, H. Davidson, D.J. & Bates, D.M. (2008). Mixed-effects modeling with crossed random efect for subjects and items. Journal of Memory and Language, 59, 390-412.
Barr, D. J., Levy, R., Scheepers, C. & Tily, H. J. (2013). Random effects structure for confirmatory hypothesis testing: Keep it maximal. Journal of Memory and Language, 68(3), 255–278.
Bolker, B.M., Brooks, M.E, Clark, C.J., Geange, S.W., Poulsen, J.R., Stevens, M.H.H. & White, J.-S.S. (2008). Generalized linear mixed models: A practical guide for ecology and evolution. Trends in Ecology and Evolution, 24, 127-135.
DeBruine, L.M. & Barr, D.J. (2021). Understanding mixed-effect models through data simulation. Advances in Methods and Practices in Psychological Science, 4, 1-15.
Fox, J. (2015). Applied regression analysis and generalized linear models. Sage Publications.
Fox, J. & Weisberg, S. (2019). An R companion to applied regression. Sage.
Galecki, A. & Burzykowski, T. (2013). Linear mixed-effects model using R. Springer New York.
Green, P., & MacLeod, C. J. (2016). “simr: an R package for power analysis of generalised linear mixed models by simulation.” Methods in Ecology and Evolution, 7(4), 493-498.
Janczyk, M. & Miller, J. (2024). Generalization of unpredictable action effect features: Large individual differences with little on-average effect. The Quarterly Journal of Experimental Psychology, 77, 898-908.
Janczyk, M., & Ulrich, R. (2019). Action consequences affect the space-time congruency effect on reaction time. Acta Psychologica, 198, 102850.
Kenward, M.G. & Roger, J. (1997). Small sample inference for fixed effects from restricted maximum likelihood. Biometrics, 53, 983-997.
Kumle, L., Võ, M. L., & Draschkow, D. (2021). Estimating power in (generalized) linear mixed models: An open introduction and tutorial in R. Behavior Research Methods, 53, 2528–2543.
Kuznetsova, A., Brockhoff, P.B. & Christensen, R.H.B. (2017). lmerTest package: Tests in linear mixed effects models. Journal of Statistical Software, 82, 1-26.
Laird, N. M. & Ware, J. H. (1982). Random-effects models for longitudinal data. Biometrics, 38(4), 963-974.
Luke S.G. (2017). Evaluating significance in linear mixed-effect model in R. Behavior Research Methods, 49, 1494-1502.
Magnusson, K. (2024). powerlmm: Power analysis for longitudinal multilevel models. R package version 0.4.0.
McNeish, D. (2017). Small sample methods for multilevel modeling: A colloquial elucidation of REML and the Kenward-Roger correction. Multivariate Behavioral Research, 52, 661-670.
Pinheiro, J. & Bates, D. (2000). Mixed-effects Models in S and S-PLUS. Springer.
Satterthwaite, F.E. (1946). An approximate distribution of estimates of variance components. Biometrics Bulletin, 2, 110-114.
Schielzeth, H., Dingemanse, N. J., Nakagawa, S., Westneat, D. F., Allegue, H., Teplitsky, C., … & Araya‐Ajoy, Y. G. (2020). Robustness of linear mixed‐effects models to violations of distributional assumptions. Methods in Ecology and Evolution, 11(9), 1141-1152.
Singmann, H., Bolker, B., Westfall, J., Aust, F., & Ben-Shachar, M. (2024). afex: Analysis of factorial experiments. R package version 1.4-1, https://CRAN.R-project.org/package=afex
Singmann, H. & Kellen, D. (2019). An introduction to mixed models for experimental psychology. In D. Spieler & E. Schumacher (Eds.), New methods in cognitive psychology (pp. 4-31). Routledge.
Stroup, W.W. (2013). Generalized linear mixed models: Modern concepts, methods and applications. CRC Press.
Verbeke, G. & Molenberghs, G. (2000). Linear mixed models for longitudinal data. Springer.
von Sobbe, L., Scheifele, E., Maienborn, C., & Ulrich, R. (2019). The space–time congruency effect: A meta-analysis. Cognitive Science, 43, 1-23.
West, B. T., Welch, K. B. & Galecki, A. T. (2022). Linear mixed models: A practical guide using statistical software. Chapman and Hall.
Westfall, J., Kenny, D. A., & Judd, C. M. (2014). Statistical power and optimal design in experiments in which samples of participants respond to samples of stimuli. Journal of Experimental Psychology: General, 143(5), 2020-2045.
Wir bleiben hier bei der bekannten Notation von Varianzen, führen aber später eine etwas generellere Schreibweise ein.↩︎
Wir können im Prinzip auch schreiben. Wir führen hier die , die den Intercept repräsentiert, aber explizit an, um die Vergleichbarkeit mit den folgenden Schreibweisen nachvollziehbarer zu machen.↩︎
Allerdings werden Modelle ohne Achsenabschnitt tatsächlich äußerst selten verwendet.↩︎
Wie oben schon erwähnt, sind LMMs generischer und die Varianz-Kovarianzmatrix der Residuen könnte auch unterschiedliche Varianzen bzw. Kovarianzen ungleich 0 enthalten. Auch wäre es denkbar, dass man individuelle Residualverteilungen pro Versuchsperson schätzt. Bei der Funktion ist dies allerdings nicht möglich. Hier muss für jede Versuchsperson die gleiche Residualverteilung, mit identischer Fehlervarianz und unabhängigen Fehlern angenommen werden.↩︎
Wenn die Auswertung von RTs mit einer herkömmlichen ANOVA erfolgt, gehen die gemittelten RTs in die Auswwertung ein. Wegen des Zentralen Grenzwertsatzes gehen diese aber gegen eine Normalverteilung, insbesondere bei großen Stichproben.↩︎
Da der Baron von Münchhausen u.a. bekundete, sich selber am Schopfe aus dem Sumpf gezogen zu haben, wird manchmal (aber echt selten), auch der Begriff “Münchhausenmethode” benutzt.↩︎
Tatsächlich geht die heutige Verbreitung von LMMs in der psychologischen Forschung auf sprachpsychologische Fragestellungen zurück (z.B. Baayen et al., 2008).↩︎