Wahlpflichtfach “Kognitions- und Klinisch-Psychologische Forschungsmethoden”
Wintersemester 2025/2026

Fortgeschrittene Verwendung von Linear Mixed-Models (LMMs)

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:

• v1.0: erste vollständige Version (6.9.2025)

1 Einführung

Nachdem wir hier die Grundlagen von Linear Mixed-Models (LMMs) gelegt haben und erste Auswertung an einem einfachen Beispiel durchgeführt haben, bauen wir nun darauf auf und behandeln exemplarisch zwei anderen, komplexere Fragestellungen. In diesem Zuge diskutieren wir auch Themen wie Modelldiagnostik etc. etwas ausführlicher.

Das erste Beispiel (Kap. 2) behandelt dabei eine längsschnittliche Untersuchung, während wir im zweiten Beispiel (Kap. 3) erneut auf die Simon-Daten zurückgreifen, diese aber mehrfaktoriell auswerten. Daran erläutern wir dann auch, wie mit Warnungen bzgl. der Konvergenz und Validität der Parameter im Zuge der Modellschätzung umgegangen werden kann.
Insgesamt ist dieser Überblick natürlich selektiv, er wird bei Notwendigkeit um weitere Ausführungen und ggf. auch um weitere Auswertungssituationen ergänzt.

Im Rahmen dieses Textes arbeiten wir mit den folgenden Paketen:

# notwendige Pakete
library(car)
library(afex)
library(lme4)
library(pbkrtest)
library(lmerTest)
library(influence.ME)
library(ez)
library(schoRsch)
library(conflicted)

2 Auswertung einer Längsschnittstudie

Als erstes Beispiel betrachten wir eine Längsschnittstudie mit einem recht einfachen Design und umgehen damit zunächst einen Teil der Probleme komplexer LMMs. Der Datensatz kann über das Paket car direkt in R geladen werden. In der (etwas älteren) Studie von Davis et al. (2005) ging es darum, den Zusammenhang zwischen sportlicher Aktivität und Anorexia nervosa (AN) zu untersuchen, da AN häufig mit erhöhter körperlicher Aktivität assoziiert ist.

Untersucht wurden 138 Kinder bzw. junge Frauen im Alter von 11 bis 18 Jahren, die im Rahmen eines Krankenhausaufenthaltes rekrutiert wurden, sowie eine entsprechende Kontrollgruppe von 98 Kindern bzw. jungen Frauen in vergleichbarem Alter. Da wir hier mehrere hundert Teilnehmerinnen betrachten, ist die Nutzung eines LMMs grundsätzlich möglich und sinnvoll, da typischerweise eie Empfehlung ist, dass es mindestens 6 Gruppen pro Gruppierungsvariable (also bspw. 6 Versuchspersonen) gibt. Ansonsten sind die geschätzten Varianz-Kovarianz-Parameter der random-effects mehr als ungenau (Bolker, 2015).

Die zentrale abhängige Variable ist eine leitfadengestützte retrospektive Schätzung der eigenen sportlichen Aktivität im Alter von 8, 10, 12 und 14 Jahren sowie zum aktuellen Zeitpunkt der Befragung. Da das Alter der Teilnehmerinnen variiert, variiert auch die Anzahl der retrospektiven Schätzungen. Außerdem variiert das Alter bei der Einschätzung der aktuellen sportlichen Aktivität. Im Durchschnitt wurde die Diagnose AN in der Patientinnengruppe ca. 1 Jahr vor dem Klinikaufenthalt gestellt. Die Daten beinhalten also Einschätzungen der sportlichen Aktivität vor und nach der offiziellen Diagnose.

Wir laden zunächst die Daten und log-transformieren die geschätzte Sportaktivität leicht, damit Extremwerte weniger ins Gewicht fallen:

library(car)
data <- Blackmore
data$exercise <- log(data$exercise + 5 / 60, base = 2)
head(data)

Die Spalte subject gibt dabei die Versuchspersonennummer an, age das Alter zur aktuellen bzw. retrospektiven Schätzung, exercise die (nun log-transformierte) Schätzung der sportlichen Aktivität (“log(Stunden pro Woche)”) und group die Zugehörigkeit zur Patientinnen- bzw. Kontrollgruppe.

2.1 Eine erste Datenexploration

Vor jeder Analyse sollten die Daten grob exploriert werden, um eine Vorstellung davon zu erhalten, wie sie sinnvoll modelliert werden könnten. Dazu schauen wir uns die Verteilung der abhängigen Variablen pro Gruppe sowie einige ausgewählte Regressionen zwischen age und exercise an.

Zunächst plotten wir ein Histogramm der Variablen exercise separat für jede Gruppe (siehe Abb. 2.1). Eine nützliche Funktion hierfür ist histogram() aus dem Paket lattice (aber natürlich gibt es etliche Alternativen):

histogram(~ exercise | group,
  data = data
)

Abbildung 2.1: Histogramm der Variablen Exercise separat für beide Gruppen.

Zu erkennen ist, dass beide Verteilungen am unteren Ende eine klare Häufung haben. Ein Teil der Teilnehmerinnen gab also an keinerlei Sport gemacht zu haben, zumindest zu einem gewissen Zeitpunkt. Dies ist bereits ein Hinweis auf ein mögliches Problem der Daten, wie wir später im Rahmen der Modelldiagnostik auch sehen werden. Da wir den Datensatz aus primär didaktischen Gründen betrachten, soll uns dies aber erst einmal nicht stören.

Unser Ziel ist letztlich die Schätzung eines LMMs, mit dessen Hilfe das Alter und die Gruppenzugehörigkeit einer Versuchsperson als Prädiktoren und die sportliche Aktivität als Kriterium betrachtet werden. Da es sich um ein Längsschnittdesign handelt und wir mehrere Werte pro Versuchsperson vorliegen haben, möchten wir mögliche interindividuelle Unterschiede berücksichtigen.

Es lohnt sich, zunächst eine Reihe einfacher, linearer Regressionen für einen Teil der Teilnehmerinnen zu betrachten, um die Angemessenheit des vorgeschlagenen LMMs besser beurteilen zu können. Zu diesem Zweck plotten wir die Daten von 20 zufällig ausgewählten Versuchspersonen pro Gruppe (da alle 236 schwer zu visualisieren sind):

# Versuchspersonen(nummern) zufällig ziehen
set.seed(1)
# Patientinnengruppe
samp_patients <- unique(data$subject[data$group == "patient"])
samp_patients <- sample(samp_patients,
  size = 20,
  replace = FALSE
)
# Kontrollgruppe
samp_control <- unique(data$subject[data$group == "control"])
samp_control <- sample(samp_control,
  size = 20,
  replace = FALSE
)

Eine nützliche Funktion zur Visualisierung solcher Daten ist xyplot(), ebenfalls aus dem Paket lattice. Abbildung 2.2 visualisiert die Daten der 20 Patientinnen und Abbildung 2.3 die Daten der 20 Kontrollpersonen.

# Plot für Patientinnengruppe
xyplot(exercise ~ age | subject,
  data = data[data$subject %in% samp_patients, ],
  main = "Patientinnengruppe",
  xlab = "Alter",
  ylab = "sp. Aktivität",
  layout = c(10, 2),
  type = c("p", "r"),
  col = "black",
  col.line = "red"
)

Abbildung 2.2: Visualisierung der Daten von 20 zufällig ausgewählten Patientinnen und der Regression von sportlicher Aktivität auf Alter.

# Plot für Kontrollgruppe
xyplot(exercise ~ age | subject,
  data = data[data$subject %in% samp_control, ],
  main = "Kontrollgruppe",
  xlab = "Alter",
  ylab = "sp. Aktivität",
  layout = c(10, 2),
  type = c("p", "r"),
  col = "black",
  col.line = "red"
)

Abbildung 2.3: Visualisierung der Daten von 20 zufällig ausgewählten Kontrollpersonen und der Regression von sportlicher Aktivität auf Alter.

Zu erkennen sind leichte Aufwärtstrends, d.h. die sportliche Aktivität scheint zumindest leicht mit dem Alter zuzuehmen, und dies ist ein wenig offensichtlicher in der Patientinnen- als in der Kontrollgruppe. Da von jeder Versuchsperson allerdings nur wenige Datenpunkte vorliegen, sollten wir nicht mehr als einen linearen Zusammenhang annehmen. Allerdings sind die (linearen) Zusammenhänge nicht gleich für jede Versuchsperson: Sie variieren sowohl im Achsenabschnitt/Intercept (= dem allgemeinen Level) als auch in der Steigung.

Wir können uns diese Sichtweise auch einmal etwas klarer machen, indem wir tatsächlich für jede Versuchsperson eine Regression berechnen. Dazu können wir sehr einfach die Funktion lmList() aus dem lme4 Paket verwenden und wir tun dies auch getrennt für die Kontroll- und die Patientinnengruppe:

# zunächst für alle Patientinnen
perSub_patients <- lmList(exercise ~ age | subject,
  data = data,
  subset = subject %in% samp_patients
)
# dann für die Kontrollgruppe
perSub_control <- lmList(exercise ~ age | subject,
  data = data,
  subset = subject %in% samp_control
)

Eine ausführliche Darstellung der Ergebnisse würden wir bekommen, indem wir summary() auf die Ergebnisse anwenden. Wir sind hier aber zunächst an den Koeffizienten interessiert (also Intercept und Slope), die wir sehr einfach mit coef() extrahieren können:

# einmal als Beispiel für die ersten sechs VPen der Patientinnengruppe
head(coef(perSub_patients))

# dann für beide Gruppen speichern:
coefficients_patients <- coef(perSub_patients)
coefficients_control <- coef(perSub_control)

Schließlich lassen wir uns die mittleren Slopes für beide Gruppen ausgeben und sehen, dass die Steigung in der Patientinnengruppe tatsächlich größer als in der Kontrollgruppe ausfällt:

mean(coefficients_patients[, 2]) # 2 = Slope, 1 = Intercept

## [1] 0.3020855

mean(coefficients_control[, 2])

## [1] 0.05069

2.2 Auswertung mit einem LMM

Nun werten wir die Daten mit einem LMM aus, bei dem age und group Prädiktoren sind und exercise das Kriterium ist. Da der allgemeine Zusammenhang zwischen age und exercise pro Stufe von group variieren könnte, lassen wir eine Interaktion zu. Weiter schätzen wir sowohl einen random intercept als auch einen random slope für die Variable age, da wir keinerlei Einschränkungen für die random-effects annehmen können.¹ Zuvor subtrahieren wir noch den Wert 8 vom Alter, welches dem niedrigsten Wert des Prädiktors age entspricht. Auf diese Art und Weise können wir den Achsensabschnitt/Intercept besser interpretieren, denn er reflektiert so die sportliche Aktivität im Alter von 8 Jahren:

data$age <- data$age - 8
lme_out <- lme4::lmer(exercise ~ group + age + group * age + (age | subject),
  data = data
)

Eine Zusammenfassung der Ergebnisse erhalten wir dann mit summary():

summary(lme_out)

## Linear mixed model fit by REML ['lmerMod']
## Formula: exercise ~ group + age + group * age + (age | subject)
##    Data: data
## 
## REML criterion at convergence: 3614.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7349 -0.4245  0.1228  0.5280  2.6362 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  subject  (Intercept) 2.08385  1.4436        
##           age         0.02716  0.1648   -0.28
##  Residual             1.54777  1.2441        
## Number of obs: 945, groups:  subject, 231
## 
## Fixed effects:
##                  Estimate Std. Error t value
## (Intercept)      -0.27602    0.18237  -1.514
## grouppatient     -0.35399    0.23529  -1.504
## age               0.06402    0.03136   2.041
## grouppatient:age  0.23986    0.03941   6.087
## 
## Correlation of Fixed Effects:
##             (Intr) grpptn age   
## grouppatint -0.775              
## age         -0.489  0.379       
## groupptnt:g  0.389 -0.489 -0.796

Schauen wir zunächst auf die fixed-effects, da diese für die zu testenden Hypothesen meistens am relevantesten sind. Um diese korrekt interpretieren zu können, müssen wir uns noch ansehen, wie die kategoriale Variable group kodiert ist:

contrasts(data$group)

##         patient
## control       0
## patient       1

Die Variable ist dummykodiert und die Kontrollgruppe ist die Referenzgruppe. Somit ergibt sich die folgende Interpretation (vgl. auch Teil 8.3.4 aus Statistik 2 sowie unsere Interpretationen aus dem ALM-Teil):

• Der Achsenabschnitt \(\hat \beta_0 = \hat \alpha = -0.28\) ist die geschätzte sportliche Aktivität im Alter von 8 Jahren in der Kontrollgruppe (da nach der Subtraktion von 8, ein kodiertes Alter von 0, ein wahres Alter von 8 Jahren impliziert).

• Der Steigungskoeffizient für group, \(\hat \beta_1 = \hat \gamma = -0.35\), ist die Differenz im Achsenabschnitt zwischen der Patientinnengruppe und der Kontrollgruppe. Im Mittel zeigt die Patientinnengruppe also eine geringere sportliche Aktivität als die Kontrollgruppe im Alter von 8 Jahren.

• Der Steigungskoeffizient für age, \(\hat \beta_3 = \hat \beta = 0.06\), ist die Steigung der Regression bzgl. age und exercise, innerhalb der Kontrollgruppe. Das heißt, erhöht man den Prädiktor age um eine Einheit, so steigt die sportliche Aktivität (in der log-Skala) in der Kontrollgruppe im Schnitt um \(0.06\) Einheiten .

• Der Interaktionsterm, \(\hat \beta_4 =\hat \delta = 0.24\), ist schließlich die Differenz in der Steigung der Regression bzgl. age und exercise, zwischen der Patientinnen- und der Kontrollgruppe. Das heißt, erhöht man den Prädiktor age um eine Einheit, so steigt die sportliche Aktivität (in der log-Skala) in der Patientinnengruppe im Schnitt um \(0.06 + 0.24 = 0.3\) Einheiten.

Mit dem Wissen aus dem Grundlagenteil können wir nun auch die restlichen Teile des summary()-Outputs klären.

• Im oberen Teil des Outputs sieht man den Aufruf und die Information, dass die Schätzung auf Basis der maximierten log-restricted-Likelihood (REML) erfolgt ist, gefolgt von der Eingabe-Formel beim Aufruf von lmer().²

• Im Anschluss gibt es eine kurze Zusammenfassung der skalierten/standardisierten Residuen . Hierbei werden die Residuen durch die geschätzte Standardabweichung der Residuen dargestellt (bei lmer() also \(\hat \sigma\), siehe auch nächster Abschnitt). Wir können diese zum Verständnis auch selbst ausrechnen:

residuen <- resid(lme_out)
sigma <- sigma(lme_out) # sigma dach
summary(residuen / sigma)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -2.7349 -0.4245  0.1228  0.0000  0.5280  2.6362

Im Wesentlichen sollten diese Werte (so wie hier) symmetrisch um Null herum sein. Unterscheiden sich die Quantile im Betrag allzu sehr, ist dies ein Hinweis auf eine asymmetrische Verteilung.

• Dann folgt die Kovarianzstruktur der random-effects mit deren Varianzen, Standardabweichungen und Korrelation. Konzeptuell ergeben diese Werte die Matrix \(\bf \mathcal D\) aus dem Abschnitt zu den mathematischen Grundlagen (wobei natürlich die Korrelation in eine Kovarianz umgerechnet werden muss). Im gleichen Block finden wir auch nochmal eine Information über die Gesamtanzahl der Beobachtungen \(N\) sowie die Anzahl der Gruppen \(K\) in der jeweiligen Gruppierungsvariable (hier also die Anzahl der Versuchspersonen). Müsste im Verlauf der weiteren Analyse die random-effect Struktur reduziert werden, dann können hier Hinweise gefunden werden, welche Parameter aus dem Modell entfernt werden könnten.

• Den Block zu den fixed-effects haben wir teilweise schon betrachtet. Hier stehen neben den eigentlichen Schätzern \(\bf {\hat \beta}\) auch die geschätzten Standardfehler der Schätzer, also \(\widehat{V(\hat \beta_j)}\) sowie \(t\)-Werte. Mögliche inferenzstatistische Tests haben wir hier eingeführt und werden darauf gleich zurückkommen.

• Schließlich gibt der letzte Abschnitt Auskunft über die Korrelationsstruktur der fixed-effects. Konkret ist hiermit die Korrelationsmatrix der fixed-effects Schätzer \(\hat \beta_j\) gemeint. Da wir die entsprechende Formel bis hierhin noch nicht besprochen haben, reichen wir diese für die Kovarianzmatrix nach: \[ \widehat{ V(\bf {\hat \beta})} = \hat \sigma \cdot \left(\sum_{i=1}^K \bf X_i' \bf{\hat V_i^{-1}} \bf X_i \right)^{-1}\,. \] Die Korrelationsmatrix kann auf Basis der Kovarianzmatrix leicht errechnet werden (bspw. mit der R-Funktion cov2cor()). Wichtig ist zu betonen, dass es sich hierbei um geschätzte Kovarianz- bzw. Korrelationsmatrix der Schätzer handelt. Das heißt, wenn wir immer und immer wieder eine Stichprobe vom gleichen Umfang ziehen würden, dann würden wir auch immer und immer wieder leicht unterschiedliche geschätzte fixed-effects erhalten. Die Kovarianz und Varianz dieser fixed-effects ist das was hier geschätzt wird. Diese Korrelation liefern z.B. Anhaltspunkte dafür, ob höhere Werte auf einem Prädiktor tendenziell auch mit höheren Werten auf einem anderen Prädiktor einhergehen (siehe hier für mehr Informationen.

Zum Schluss dieses Abschnitts wollen wir noch die Ergebnisse visualisieren (siehe Abb. @rtef(fig:ergebnis-beispiel-1). Hierfür erstellen wir zuerst einen Scatter-Plot aller Daten und zeichnen anschließend die vorhergesagten Werte mit Hilfe der Funktion predict() ein. Das Argument re.form = NA in diesem Zusammenhang gibt an, dass wir an der Modellvorhersage über alle Versuchspersonen hinweg interessiert sind; ansonsten müsste man auch noch eine konkrete subject-Nummer übergeben:

# Farben definieren, macht den nachfolgenden Code etwas leserlicher
red <- rgb(
  red = 1, green = 0, blue = 0, # rot in RGB Koordinaten
  alpha = 0.2
) # Transparenz
blue <- rgb(
  red = 0, green = 0, blue = 1, # blau in RGB Koordinaten
  alpha = 0.2
) # Transparenz

# Grundlegender Scatter-Plot mit Gruppe farbig kodiert
plot(
  data$exercise ~ jitter(data$age,
    amount = 0.1
  ),
  col = ifelse(data$group == "patient", red, blue),
  pch = 19,
  main = "Zusammenhang Alter und sp. Aktivität \n für beide Gruppen",
  ylab = "Aktivität (log-Skala)",
  xlab = "Alter",
  ylim = c(-3, 7)
)

# Legende dazu
legend(
  x = "topright",
  legend = c("AN-Gruppe", "Kontrollgruppe"),
  col = c(red, blue),
  pch = 19
)

# Vorhersage des Modells
# ... benötigt einen Datensatz mit Prädiktorwerten
new_data <- data.frame(
  age = rep(seq(0, 9, 1), 2),
  group = rep(c("control", "patient"),
    each = 10
  )
)
# ... dann die Vorhersage anfordern über alle Personen hinweg
new_data$pred_exercise <- predict(lme_out,
  new_data,
  re.form = NA
)

Ein approximatives Konfidenzintervall für die Vorhersage kann auf Basis des Standardfehlers der Vorhersage hergeleitet werden. Den Standardfehler erhalten wir über das Argument se.fit der predict() Funktion. Allerdings sollte hierbei beachtet werden, dass die Schätzunsicherheit in den random-effects hierbei nicht berücksichtigt wird.

conf_interval <- predict(lme_out,
  new_data,
  re.form = NA,
  se.fit = TRUE
)
new_data$pred_exercise <- conf_interval$fit # Vorhersage
new_data$lower <- conf_interval$fit - conf_interval$se.fit * 1.96 # untere Grenze
new_data$upper <- conf_interval$fit + conf_interval$se.fit * 1.96 # obere Grenze

Nun können wir die Regressionsgeraden und die approximativen Konfidenzintervalle im Scatterplot ergänzen:

# Patientinnengruppe
points(pred_exercise ~ age, # Regressionsgerade
  data = subset(
    new_data,
    group == "patient"
  ),
  type = "l",
  col = "red"
)
points(lower ~ age, # untere Grenze des KIs
  data = subset(
    new_data,
    group == "patient"
  ),
  type = "l",
  lty = 2,
  col = "red"
)
points(upper ~ age, # obere Grenze des KIs
  data = subset(
    new_data,
    group == "patient"
  ),
  type = "l",
  lty = 2,
  col = "red"
)

# Kontrollgruppe
points(pred_exercise ~ age, # Regressionsgerade
  data = subset(
    new_data,
    group == "control"
  ),
  type = "l",
  col = "blue"
)
points(lower ~ age, # untere Grenze des KIs
  data = subset(
    new_data,
    group == "control"
  ),
  type = "l",
  lty = 2,
  col = "blue"
)
points(upper ~ age, # obere Grenze des KIs
  data = subset(
    new_data,
    group == "control"
  ),
  type = "l",
  lty = 2,
  col = "blue"
)

Abbildung 2.4: Visualisierung der Ergebnisse der LMM-Auswertung der Längsschnittdaten.

2.3 Inferenzstatistische Auswertung der fixed-effects

Beispielhaft wenden wir nun die Korrektur nach Kenward und Roger (1997) an, um p-Werte zu den fixed-effects zu erhalten:

Anova(lme_out,
  type = 3,
  test.statistic = "F"
)

Mit Bezug auf die Daten des Beispiels müssten wir also schließen:

1. (Intercept): Es gibt nicht genügend Evidenz, um anzunehmen, dass die Kontrollgruppe im Alter von 8 Jahren eine von 0 verschiedene körperliche Aktivität (auf der Log-Skala!) aufweist (der Intercept).

2. Prädiktor group: Es gibt nicht genug Evidenz um davon auszugehen, dass sich Kontroll- und Patientinnengruppe im Alter von 8 Jahren in ihrer sportlichen Aktivität unterscheiden.

3. Prädiktor age: Es gibt einen signifikanten Zusammenhang zwischen dem Alter einer Person und der sportlichen Aktiviät in der Kontrollgruppe.

4. Interaktion group:age: AN-Patientinnen haben einen stärkeren Anstieg in der sportlichen Aktivität über das Jugendalter hinweg als die Personen der Kontrollgruppe. Ggf. könnte also eine überdurchschnittlich starke Zunahme in der sportlichen Aktivität im Jugendalter ein Kennzeichen einer später ausgebildeten AN sein (so zumindest diskutieren die Autor:innen der Orginalstudie ihre Daten).

2.4 Modelldiagnostik

Eine korrekte Inferenz über random- und fixed-effects setzt die Angemessenheit der statistischen Annahmen voraus. Zur Beurteilung ist es unerlässlich sich einige Plots anzusehen, anhand derer diese Angemessenheit beurteilt werden kann. Im Gegensatz zu klassischen Regressionen im Zuge des ALMs, ist das Vorgehen bei LMMs aufgrund derer Komplexität nicht wirklich standardisiert und es gibt verschiedene Empfehlungen für notwendige/sinnvolle diagnostische Plots. Auch gibt es immer noch Forschung und Diskussion darüber, wie Modelldiagnostik sinnvoll vollzogen werden kann. Wir werden daher hier bei grundlegenden Checks bleiben. Für mehr Informationen siehe bspw. Claeskens (2013) oder West et al. (2022).

2.4.1 Residuendiagnostik

Sog. unskalierte (englisch: “raw”) “conditional” bzw. bedingte Residuen beziehen sich auf die Residuen unter Berücksichtigung der random-effects. Sie beschreiben die Abweichungen der Beobachtungen relativ zu den vorhergesagten Werten auf Basis der fixed- und random-effects: \[ \bf e_{(c),i} = \bf y_i - X_i \bf {\hat \beta} - Z_i \hat b_i\,. \] Diese können dazu genutzt werden um die Heterogenität der (bedingten) Residuen (= \(\varepsilon_i\)) zu prüfen. Die rohen Residuen eignen sich allerdings u.U. nicht, um die Normalverteiltheit von \(\varepsilon_i\) oder die Anwesenheit von “Ausreißern” zu prüfen, da diese ggf. korreliert sind oder unterschiedliche Varianzen besitzen (je nachdem wie das LMM bzgl. \(V(\varepsilon_i) = \bf {\mathcal R_i}\) parametrisiert ist; kein Problem bei lmer(), da lmer() stets unkorrelierte Fehler mit gleicher Varianz annimmt). Sollten sich die Varianzen der Residuen pro Gruppe unterscheiden, bietet es sich an, die Residuen zu standardisieren (man unterscheidet dann zwischen Pearson und studentisierten Residuen, wobei bei Pearson Residuen die Variabilität in \(\bf {\hat \beta}\) ignoriert wird). Gibt es zusätzlich noch Korrelationen bzgl. \(\varepsilon_i\), bietet es sich an, die Residuen noch entsprechend zu transformieren (siehe Galecki & Burzykowski, 2013, S. 266).

Darüber hinaus unterscheidet man unskalierte (“raw”) “marginal” Residuen (bzw. deren standardisierte Form; diese sind über das Paket HLMdiag verfügbar), bei denen die Beobachtungen relativ zu den vorhergesagten Werten auf Basis der fixed-effects betrachtet werden, also \[\bf e_{(m),i} = \bf y_i - X_i \bf{\hat \beta}\,.\]

Diese Residuen dargestellt gegenüber einzelnen Prädiktoren können dazu genutzt werden, um die Linearität des Modells zu beurteilen.

Zuerst betrachten wir die “raw conditional” (“rohen bedingten”) Residuen als Histogram, separat für die Patientinnen und die Kontrollgruppe (Abb. 2.5) sowie einen entsprechenden QQ-Plot gegenüber einer Normalverteilung (Abb. 2.6):

# raw conditional Residuen
data$res <- resid(lme_out) # direkt anhängen zum leichteren plotten

# Histogramm der conditional Residuen
histogram(~ res | group,
  data = data,
  layout = c(2, 1),
  xlab = "bedingte Residuen"
)

Abbildung 2.5: Histogramm der raw conditional Residuen.

# QQ-Plot gegenüber einer Normalverteilung
qqmath(~ res | group,
  data = data,
  prepanel = prepanel.qqmathline, # ... ab hier nur
  panel = function(x, ...) { # ... für die eingezogene Linie
    panel.qqmathline(x, ...)
    panel.qqmath(x, ...)
  }
)

Abbildung 2.6: QQ-Plot mit Normalverteilung als Referenzverteilung.

Zu erkennen ist, dass die bedingten Residuen nicht wirklich normalverteilt sind, sondern eher rechtssteil mit tendenziell zu großen negativen Residuen. Als nächstes visualisieren wir die bedingten Residuen gegenüber den vorhergesagten Werten, um zu prüfen, ob die Varianzen heterogen sind (siehe Abb. 2.7).³

plot(lme_out,
  resid(.) ~ fitted(.) | group,
  layout = c(2, 1),
  aspect = 2,
  abline = 0
)

Abbildung 2.7: Scatterplot der bedingten Residuen als Funktion der vorhersagten Werte.

Hier sehen wir ein unschönes Muster, welches auch die großen negativen Residuen erklärt. Zwar sind die Residuen teilweise heterogen verteilt, allerdings gibt es in der linken unteren Ecke der beiden Diagramme in Abbildung 2.7 Residuen, die auf einer Art diagonalen “Linie” verlaufen. Der Grund hierfür ist, dass manche Personen keine sportliche Aktivität angegeben haben, sodass sich bzgl. des Kriteriums eine untere Grenze ergibt (sportliche Aktivität kann nun mal nicht negativ sein). Diese “Ballung” von Aktivitätswerten am unteren Ende hatten wir auch bei der explorativen Datenanalyse zu Beginn festgehalten (siehe Abb. 2.1 und hier tritt sie nun als Problem im Zuge des LMMs ans Tageslicht. Klare Ausreißer sind in den Plots allerdings nicht erkennbar.

Schließlich betrachte wir noch die unskalierten (“raw”) marginal Residuen und plotten diese als Funktion der Prädiktoren (siehe Abb. 2.8):

data$raw_res <- data$exercise - predict(lme_out,
  re.form = NA
)
xyplot(raw_res ~ age | group, data,
  type = c("p", "r")
)

Abbildung 2.8: Unskalierte marginal Residuen als Funktion der Prädiktoren.

Hier sind die Daten nahezu gleichmäßig in \(y\)-Richtung verteilt (abgesehen von den wenigen Datenpunkten rechts unten im rechten Plot), sodass es keine großen Zweifel an der Linearität der Regression gibt. Auch hier scheint es keine neuen Hinweise auf weitere Ausreißer zu geben (abgesehen von den bereits angesprochenen, teilweise großen, negativen Residuen durch Teilnehmerinnen mit wenig sportlicher Aktivität)

2.4.2 Einflussreiche Datenpunkte

Neben dem Prüfen nach reinen Ausreißern ist es stets sinnvoll auch nach einflussreichen Datenpunkten zu suchen. Wie auch im Fall der normalen Regression im Zuge des ALMs, bestimmt sich der Einfluss eines Datenpunkts nicht nur durch sein Residuum, sondern durch eine Kombination aus “Leverage” und Residuum (siehe hier). Da die Daten im Rahmen von LMMs allerdings oft auf mehreren “Ebenen” betrachtet werden können (bspw. bezüglich einzelner Beobachtungen oder bzgl. ganzer Gruppen) und zudem einflussreiche Datenpunkte sich sowohl auf die fixed-effects auch auf die Parameter der random-effects bzw. der Residuen auswirken können, gibt es auch hier mehrere Ansätze (siehe West et al., 2022, S. 43ff für eine Übersicht) und nicht alle Ansätze sind “out-of-the-box” in R verfügbar.

Ein durchaus verbreiteter Ansatz ist es, gewisse Gruppen (also bspw. ganze Versuchspersonen) nacheinander wegzulassen und auf jeden so erhaltenen Teildatensatz das Modell neu zu fitten. Analysiert man dann, wie stark das Weglassen einer Gruppe zu Veränderungen in den geschätzten Parametern führt, kann der Einfluss einzelner Gruppen identifiziert werden. Eine Metrik ist hierbei wieder Cook’s Distance, bei der die Differenz der Parameter zwischen dem Modell-Fit mit allen Daten und den Daten ohne eine Gruppe relativ zur geschätzten Variabilität der Parameter zu einer Kenngröße verrechnet wird (für mehr Informationen siehe auch Nieuwenhuis et al., 2012). Eine Funktion zur Berechnung der Cook’s Distance für die fixed-effects liefert das Paket influence.ME (für den Einfluss auf die random effects können sog. DFBETAs für einzelne Parameter berechnet werden; siehe die Funktion dfbetas()).⁴ Da hierbei aber pro Gruppe das Modell neu geschätzt wird, kann es je nach Komplexität des Modells etwas dauern, bis die Funktion fertig gerechnet hat. Als Daumenregel kann ein Datenpunkt als einflussreich gesehen werden, wenn dessen Cook’s Distance größer ist als \(4/K\), wobei \(K\) die Anzahl der Gruppen pro Gruppierungsvariable ist (in unserem Beispiel also die Anzahl der Versuchspersonen). Eine Visualisierung ist in Abbildung 2.9 zu finden.

library(influence.ME)
estex_lme <- influence(lme_out, # lmer-Output
  group = "subject"
) # Welche Gruppen sollen weggelassen werden?
n <- length(unique(data$subject)) # Anzahl der Versuchspersonen
plot(estex_lme,
  which = "cook", # Cook's Distance soll geplottet werden
  cutoff = 4 / n, # cut-off Daumenregel
  sort = TRUE,
  xlab = "Cook´s Distance",
  ylab = "Versuchsperson"
)

Abbildung 2.9: Visualisierung von Cook’s Distance für das LMM der Längsschnittdaten.

In diesem Datensatz ergeben sich Ausreißer (die roten Dreiecke) und die genauen Werte pro Versuchsperson lassen sich auch ausgeben:

all_cooks_dists <- cooks.distance(estex_lme,
  sort = TRUE
)
tail(all_cooks_dists,
  n = 10
) # die höchsten 10 Werte

##            [,1]
## 198  0.01503385
## 311  0.01507431
## 151  0.01511313
## 279b 0.01512782
## 253  0.01602581
## 241  0.01746677
## 251  0.02122207
## 212  0.02511654
## 284  0.03007137
## 303  0.03259715

Auf Basis dieser Ergebnisse lohnt es sich, die Versuchspersonen 241, 251, 212, 284 und 303 genauer anzusehen und evtl. als Ausreißer zu entfernen. In diesem Zuge lohnt es sich auch ggf. die Cook’s Distance erneut zu berechnen, aber diesmal separat für einzelne Parameter um so herauszubekommen, auf welchen Teil des Modells eine Versuchsperson bzw. Gruppe besonderen Einfluss nimmt.

Anmerkung zu den Daten: Im vorletzten Code-Abschnitt werden eigentlich mehrere Warnungen geworfen, die anzeigen, dass das Modell nicht sauber konvergiert, wenn einzelne Gruppen (also hier: Teilnehmerinnen) aus den Daten entfernt werden; wir haben allerdings hier die Ausgabe der Warnungen unterdrückt. Allerdings sind diese Warnungen ein klarer Hinweis darauf, dass die verwendeten Daten nicht ganz ideal für das verlangte LMM sind und ggf. neben der Ausreißerdiagnostik sowie Maßnahmen gegen die teilweise sehr großen negative Residuen (siehe vorherigen Abschnitt) auch eine Reduktion der Komplexität der random-effect Struktur notwendig ist (siehe hierzu auch das nächste Beispiel in Kap. 3). Beispielsweise werden die Warnungen deutlich seltener, wenn man die Korrelation zwischen den random-effects nicht berücksichtigt. Hieran erkennt man auch, dass es meist mehrere “Iterationen” von Modelldiagnostik, Datentransformation bzw. -bereinigung und Modellspezifikation gibt, bis dann ein schließlich finales Modell existiert und berichtet werden kann.

2.4.3 Verteilung der random-effects

Als letztes schauen wir noch die Verteilung der random-effects an. Ein häufiger Plot ist hierbei ein QQ-Plot eines jeden random-effects gegenüber einer Normalverteilung als Referenz. Dabei ist allerdings anzumerken, dass dieser Plot nicht unbedingt dazu genutzt werden sollte, um die Normalverteiltheit der random-effects zu prüfen, denn auf Stichprobenebene müssen diese nicht zwangsläufig normalverteilt sein (vgl. Verbeke & Molenberghs, 2000). Dennoch sind solche Plots sinnvoll, um mögliche Ausreißer oder bimodale, nicht-kontinuierliche Verteilungen zu entdecken. Ein letztendlicher Test der Normalverteiltheit kann durch einen Modellvergleich erzielt werden, bei dem ein Modell mit Normalverteilungsannahme gegenüber einem Modell ohne Normalverteilungsannahme getestet wird (siehe auch hier Verbeke & Molenberghs, 2000). In R kann ein QQ-Plot der random effects wie folgt anfordern (siehe Abb. 2.10):

qqmath(ranef(lme_out)) ## for comparison

## $subject

Abbildung 2.10: QQ-Plot der random-effects mit einer Normalverteilung als Referenz.

2.4.4 Mögliche Gegenmaßnahmen bei Ausreißern oder unerwünschten Trends der Residuen

Werden in einem LMM unerwünschte Eigenschaften entdeckt, so gibt es leider keine allgemeinen Regeln, wie man mit diesen umgeht. Grundsätzlich kann aber empfohlen werden, weitere Analysen anzustreben und darüber nachzudenken, was die Ursache für die möglichen Probleme sind.

• Bei Ausreißern könnte es bspw. sein, dass Werte einer Versuchsperson falsch kodiert wurden. Ggf. lohnt es sich auch die Komplexität des Modells zu reduzieren, damit es weniger anfällig für einzelne einflussreiche Gruppen/Beobachtungen wird.

• Bei nicht-linearen Trends könnte es bspw. sein, dass weitere Prädiktoren aufgenommen oder ggf. Interaktionen zwischen Prädiktoren zugelassen werden müssen.

• Bei Trends in den Varianzen der Residuen müssen ggf. weitere Kovarianzparameter geschätzt werden, um diese Trends akkurat zu berücksichten.

• Auch kann es sein, dass Versuchspersonen schlichtweg ausgeschlossen werden müssen, weil sie Daten produziert haben, die nicht typisch sind (z.B. könnte es sein, dass eine Versuchsperson aufgrund von Übernächtigung besonders lange RTs produziert hat).

Sollte es allerdings am Ende der Fall sein, dass die Daten an sich zu sehr von den Modellannahmen abweichen, muss ggf. vom klassischen LMM abgewichen werden. Beispielsweise könnte es sein, dass man eine andere Residualverteilung modellieren muss, da diese eher Gamma- als normalverteilt ist. Dann müssen sog. generalized linear mixed models (GLMMs) verwendet werden, die mit der Funktion glmer() berechnet werden können. Alternativ kann es auch sein, dass man zu non-parametrischen Schätzverfahren wechseln sollte (siehe hierfür bspw. Waterman et al., 2015).

2.5 Exkurse mit weiteren Informationen und Erläuterungen

2.5.1 Wie und warum sind Beobachtungen nun korreliert?

Nachdem wir die Grundlagen von LMMs und deren Tests und Diagnostik nun behandelt haben, möchten wir nun einmal genauer beleuchten, warum und zu welchem Ausmaß die Berücksichtigung von random-effects eigentlich dazu führt, dass Korrelationen zwischen Beobachtungen modelliert werden. Einer der wesentlichen Gründe für die Verwendung von LMMs war ja die Notwendigkeit, korrelierte Beobachtungen sauber zu berücksichtigen. In RT-Studien unterscheiden sich bspw. die Versuchspersonen in ihrem allgemeinen RT-Level, sodass RTs einer Person sich ähnlicher sind als solche verschiedener Versuchspersonen, wodurch sie dann korrelieren.

Zur Klärung, wie die Korrelationen ins Spiel kommen sobald random-effects berücksichtigt werden, betrachten wir eine Versuchsperson im AN-Beispieldatensatz etwas genauer. Grundsätzlich haben wir festgehalten, dass ein LMM als \[Y_i = {\bf X_i} {\bf \beta} + {\bf Z_i} b_i + \varepsilon_i \quad i\in\{1,..,K\}\] formuliert wird, wobei \(b_i \sim N_q({\bf 0}, {\bf {\mathcal D}})\) und \(\varepsilon_i \sim N_{n_i}({\bf 0}, {\bf {\mathcal R_i}})\).

Wie ist dieses Modell im vorliegenden Fall für eine Versuchsperson genau parametrisiert? Schauen wir uns hierfür zuerst die (Prädiktor-)Daten für die erste Teilnehmerin an (hier mit subject-Nummer 100; wir werden sie aber im folgenden als Versuchsperson 1 bezeichnen):

data[data$subject == 100, ]

Da in unserem vorherigen LMM die Gruppenzugehörigkeit, das Alter und deren Interaktion als fixed-effects betrachtet wurden, ergibt sich die folgende Designmatrix auf Basis der Daten (mit den Spalten für den Achsenabschnitt, den Prädiktor “Gruppe”, das “Alter” und deren Interaktion als multiplikative Verknüpfung): \[\bf X_1 = \begin{pmatrix} 1 & 1 & 0 & 0 \\ 1 & 1 & 2 & 2 \\ 1 & 1 & 4 & 4 \\ 1 & 1 & 6 & 4 \\ 1 & 1 & 7.92 & 7.92 \end{pmatrix}\,.\] Zur Erinnerung: Die erste Einser-Spalte modelliert den Achsenabschnitt und die zweite Einserspalte ergibt sich, da die Gruppenzugehörigkeit dummykodiert ist und die Kontrollgruppe dabei die Referenzgruppe darstellt (d.h. die betrachtete Versuchsperson erhält für den Prädiktor “Gruppe” eine 1). Die Matrix der random-effects \(\bf Z_1\) erhält die gleichen Spalten wie \(\bf X_1\) für den Achsenabschnitt und den Prädiktor “Alter”, da wir einen random intercept und einem random slope für das Alter modelliert haben. \[\bf Z_1 = \begin{pmatrix} 1 & 0 \\ 1 & 2 \\ 1 & 4 \\ 1 & 6 \\ 1 & 7.92 \end{pmatrix}\,.\]

Ausgeschrieben ergibt sich damit folgende Grundgleichung: \[\begin{pmatrix} Y_{1,1} \\ Y_{1,2} \\ Y_{1,3} \\Y_{1,4} \\ Y_{1,5} \end{pmatrix} = \begin{pmatrix} 1 & 1 & 0 & 0 \\ 1 & 1 & 2 & 2 \\ 1 & 1 & 4 & 4 \\ 1 & 1 & 6 & 4 \\ 1 & 1 & 7.92 & 7.92 \end{pmatrix} \cdot \begin{pmatrix} \beta_{1} \\ \beta_{2} \\ \beta_{3} \\ \beta_{4} \\ \beta_{5} \end{pmatrix} + \begin{pmatrix} 1 & 0 \\ 1 & 2 \\ 1 & 4 \\ 1 & 6 \\ 1 & 7.92 \end{pmatrix} \cdot \begin{pmatrix} b_{1,1} \\ b_{1,2} \end{pmatrix} + \begin{pmatrix} \varepsilon_{1,1} \\ \varepsilon_{1,2} \\ \varepsilon_{1,3} \\\varepsilon_{1,4} \\ \varepsilon_{1,5} \end{pmatrix}\,.\]

Bei der Nutzung der Funktion lmer() werden keinerlei Restriktionen für die Varianz-Kovarianzmatrix der random-effects eingeführt. Es gilt also ganz allgemein: \[\bf{\mathcal D} = \begin{pmatrix} \Psi_{11} & \Psi_{12} \\ \Psi_{21} & \Psi_{22} \end{pmatrix} \quad \text{mit} \quad \Psi_{12} = \Psi_{21}\,.\] Außerdem geht lmer() stets davon aus, dass alle Residuen die gleiche Varianz besitzen und paarweise unkorreliert sind, also \[\bf {\mathcal R_1} = \sigma^2 \cdot I_5 = \sigma^2 \begin{pmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1\end{pmatrix}\,.\]

Was gilt nun für die (geschätzte) Varianz-Kovarianzmatrix der Beobachtungen für diese Versuchsperson, wenn man sie nicht isoliert, sondern in der Gesamtheit aller Daten betrachtet? Die allgemeine Formel haben wir bei den Details zur Schätzung des klassischen LMMs ausgelagert (siehe hier) Auf Populationsebene gilt \[V(Y_i) = {\bf {\mathcal V_i}} = {\bf Z_i}{\bf {\mathcal D} } {\bf Z_i}'+ {\bf {\mathcal R_i}}\] und analog dazu gilt auf Stichprobenebene \[{\bf {\hat {\mathcal V_i}}} = {\bf Z_i}{\bf {\hat {\mathcal D}}} {\bf Z_i}'+ {\bf {\mathcal {\hat R}_i}}\,.\]

Die Formel für jedes Element der (geschätzten) Kovarianzmatrix des Kriteriums formulieren wir hier lieber nicht explizit aus, da dies in eine recht mühsame Schreibarbeit ausarten würde (man findet sie aber in Fox, 2015, S. 722). Es sei aber angemerkt, dass das Ergebnis keine Diagonalmatrix sein wird, da auch \({\bf Z_i}{\bf {\hat {\mathcal D}}} {\bf Z_i}'\) keine Diagonalmatrix ergibt. In anderen Worten kommt über die Kovarianzmatrix der random-effects in Kombination mit den Werten in der zugehörigen Matrix \(\bf Z_i\) eine Korrelation zwischen den Beobachtungen einer Person ins Spiel.

Schauen wir uns nun die Werte für \(\bf{\hat {\mathcal V_i}}\) in R für Versuchsperson 1 genauer an. Zunächst kodieren wir \(\bf Z_1\), \(\bf {\hat {\mathcal D}}\), und \(\bf {\hat {\mathcal R_1}}\). Es gibt auch elegantere und exaktere Varianten diese Informationen direkt aus dem lmer()-Objekt zu extrahieren, aber da diese nicht ganz so intuitiv sind, belassen wir es bei der eher “händischen” Umsetzung (mehr hierzu findet sich allerdings bei Galecki & Burzykowski, 2013, S. 370):

# Z_1 bauen
Z_1 <- matrix(c(1, 0, 1, 2, 1, 4, 1, 6, 1, 7.92),
  nrow = 5,
  ncol = 2,
  byrow = TRUE
)
Z_1

##      [,1] [,2]
## [1,]    1 0.00
## [2,]    1 2.00
## [3,]    1 4.00
## [4,]    1 6.00
## [5,]    1 7.92

# D (geschwungen geschrieben) bauen
VarCorr(lme_out)

##  Groups   Name        Std.Dev. Corr  
##  subject  (Intercept) 1.4436         
##           age         0.1648   -0.281
##  Residual             1.2441

# hieraus ergeben sich die Varianzen und auf Basis der Korrelation auch die Kovarianz
D <- matrix(
  c(
    1.4436^2, # Varianz random intercept
    -0.281 * 1.4436 * 0.1648, # Kovarianz
    -0.281 * 1.4436 * 0.1648, # Kovarianz
    0.1648^2
  ), # Varianz random slope
  nrow = 2,
  ncol = 2
)
D

##             [,1]        [,2]
## [1,]  2.08398096 -0.06685138
## [2,] -0.06685138  0.02715904

# R_1 (geschwungen geschrieben) bauen
sigma <- sigma(lme_out)
R_1 <- sigma^2 * diag(1,
  nrow = 5,
  ncol = 5
)
R_1

##          [,1]     [,2]     [,3]     [,4]     [,5]
## [1,] 1.547773 0.000000 0.000000 0.000000 0.000000
## [2,] 0.000000 1.547773 0.000000 0.000000 0.000000
## [3,] 0.000000 0.000000 1.547773 0.000000 0.000000
## [4,] 0.000000 0.000000 0.000000 1.547773 0.000000
## [5,] 0.000000 0.000000 0.000000 0.000000 1.547773

# vom Modell implizierte Kovarianzmatrix (Z_1' * D * Z_1 + R_1)
V_1 <- Z_1 %*% D %*% t(Z_1) + R_1
V_1

##          [,1]     [,2]     [,3]     [,4]     [,5]
## [1,] 3.631754 1.950278 1.816575 1.682873 1.554518
## [2,] 1.950278 3.472985 1.900145 1.875078 1.851014
## [3,] 1.816575 1.900145 3.531488 2.067284 2.147511
## [4,] 1.682873 1.875078 2.067284 3.807263 2.444007
## [5,] 1.554518 1.851014 2.147511 2.444007 4.276417

# Korrelationen
R_y <- cov2cor(V_1)

# Zeilen bzw. Spalten nach den Messzeitpunkten benennen
colnames(R_y) <- rownames(R_y) <- paste0("t_", Z_1[, 2])
R_y

##              t_0       t_2       t_4       t_6    t_7.92
## t_0    1.0000000 0.5491448 0.5072434 0.4525709 0.3944550
## t_2    0.5491448 1.0000000 0.5425708 0.5156577 0.4803064
## t_4    0.5072434 0.5425708 1.0000000 0.5637868 0.5526071
## t_6    0.4525709 0.5156577 0.5637868 1.0000000 0.6056979
## t_7.92 0.3944550 0.4803064 0.5526071 0.6056979 1.0000000

Die geschätzten bzw. implizierten Korrelationen durch die Berücksichtigung der random-effects sind relativ hoch und positiv: Personen die höher in ihrer sportlichen Aktivität zu Kindeszeiten waren, sind es auch später im Jugendalter. Außerdem ergibt sich hier ein interessantes Muster das intuitiv auch Sinn ergibt, denn die Korrelationen zwischen zwei Messungen sind höher, je näher sich die Messzeitpunkte sind. Wichtig ist allerdings zu betonen, dass diese Struktur nicht immer genau so aussieht, sondern (wie in der oberen Gleichung impliziert) von den geschätzten Varianz- und Kovarianzparametern der random-effects sowie den Werten in \(\bf Z_i\) abhängt. Aus diesem Grund ist die Kovarianz zwischen zwei Messzeitpunkten, die gleich weit auseinander liegen, auch nicht gleich. Auch ändert sich das Korrelationsmuster drastisch, wenn random-effects ein- bzw. ausgeschlossen werden. Es ergibt sich z.B. eine sog. compound symmetry-Struktur mit gleichen Kovarianzen und Varianzen, wenn der random slope nicht mehr berücksichtigt wird:

# kein random slope
lme_out_test <- lme4::lmer(exercise ~ group + age + group * age + (1 | subject),
  data = data
)

# Z_1 bauen (nur noch ein Einsvektor, da nur noch ein random intercept)
Z_1 <- matrix(c(1, 1, 1, 1, 1),
  nrow = 5,
  ncol = 1
)

# D (geschwungen geschrieben) "bauen", was diesmal aber nur eine Varianz ist
VarCorr(lme_out_test)

##  Groups   Name        Std.Dev.
##  subject  (Intercept) 1.3695  
##  Residual             1.3443

D <- matrix(1.3695^2,
  nrow = 1,
  ncol = 1
)

# R_1 (geschwungen geschrieben) bauen
sigma <- sigma(lme_out_test)
R_1 <- sigma^2 * diag(1,
  nrow = 5,
  ncol = 5
)

# vom Modell implizierte Kovarianzmatrix (Z_1' * D * Z_1 + R_1)
V_1 <- Z_1 %*% D %*% t(Z_1) + R_1
V_1

##          [,1]     [,2]     [,3]     [,4]     [,5]
## [1,] 3.682632 1.875530 1.875530 1.875530 1.875530
## [2,] 1.875530 3.682632 1.875530 1.875530 1.875530
## [3,] 1.875530 1.875530 3.682632 1.875530 1.875530
## [4,] 1.875530 1.875530 1.875530 3.682632 1.875530
## [5,] 1.875530 1.875530 1.875530 1.875530 3.682632

2.5.2 Warum lmer() nicht ideal für Längsschnittstudien ist

Im vorherigen Abschnitt haben wir gesehen, dass durch die Berücksichtigung von random-effects Korrelationen zwischen Beobachtungen impliziert bzw. modelliert werden können. Wir haben auch festgehalten, dass die genaue Struktur von \(\bf Z_i\), den geschätzten Varianzen und Kovarianzen in \(\bf {\mathcal D}\) sowie den Varianzen und Kovarianzen in \(\bf{\mathcal R_i}\) abhängt: \[V(Y_i) = \bf{\mathcal V_i} = {\bf Z_i}{\bf {\mathcal D}} {\bf Z_i}'+ {\bf {\mathcal R_i}}\,.\] Hier wird nun eine Limitation von lmer() deutlich, denn mit lmer() können wir die Matrix \(\bf {\mathcal R_i}\) nicht beliebig parametrisieren; es gilt stets \(\bf { \mathcal R_i} = \sigma^2 I_{n_i}\). In anderen Worten wird die Korrelations- bzw. Kovarianzstruktur nur vom Term \({\bf Z_i}{\bf{ \mathcal D}} {\bf Z_i}'\) beeinflusst, wobei das Ergebnis dieses Produkts nicht wirklich intuitiv abschätzbar oder kontrollierbar ist, sobald mehrere random-effects und kontinuierliche Prädiktoren in \(\bf Z_i\) vorliegen. Dies wird dann zum Problem, wenn man mehr Kontrolle haben möchte über die implizierten Korrelationen, bspw. bei der statistischen Auswertung von Längsschnittstudien, bei denen die modellierten Korrelationen zwischen zwei Zeitpunkten wesentlich sind. In diesen Fällen ist es ggf. notwendig auf die Funktion lme() aus dem nlme Paket auszuweichen, mit deren Hilfe man explizit die Kovarianzstruktur von \({\bf { \mathcal D}}\) und \({\bf {\mathcal R_i}}\) definieren kann. Beispielsweise könnte man einen random intercept schätzen mit einer gleichzeitigen sog. autoregressiven Korrelations- bzw. Kovarianzstruktur bzgl. \({\bf {\mathcal R_i}}\), bei dem die Korrelationen zwischen Zeitpunkten exponentiell abnehmen. Wie solche Fälle dann parametrisiert und in R umgesetzt werden findet sich in Galecki und Burzykowski (2013, Kap. 10 und 14).

3 Auswertung einer komplexeren Experimentalstudie

Mit einem zweiten Beispiel wollen wir uns einen experimentalpsychologischen Datensatz anschauen und hierbei auf unser Vorwissen auf Basis des ersten Datensatzes aufbauen. Insbesondere wollen wir bei diesem Datensatz beleuchten, wie man idealerweise damit umgeht, wenn Modelle nicht fehlerfrei geschätzt werden (und Warnungen ausgegeben werden).

Wir verwenden dazu wieder die Daten einer Simon-Aufgabe, die wir bereits in Teil 1 zu Mathematischen Modellen und Simulationen und in den Grundlagen zu LMMs genutzt haben. Dieser Datensatz umfasst 32 Versuchspersonen, die jeweils mehr als 200 Trials bearbeitet haben (man bedenke, das man mehr als 6 Gruppen pro Gruppierungsvariable braucht um ein LMM sinnvoll zu schätzen; also hier mehr als 6 Versuchspersonen). In der Regel treten Probleme und Warnungen im Modellfit dann auf, wenn das Modell zu komplex für die gegebenen Daten ist. Etwas salopp gesagt verlangt man den Daten für ihre Menge zu viel ab. Aus diesem Grunde führen wir in unseren Simon-Datensatz zwei neue Variablen ein, nämlich die Kongruenz des vorherigen Trials \(N-1\) und die Kongruenz des vorvorherigen Trials \(N-2\). Auf diese Art und Weise haben wir hinreichend viele Prädiktoren bzw. random-effects, um das geschätzte LMM einmal gezielt an dessen Grenzen zu bringen.

data <- read.table("./Daten/daten_Simon.txt",
  sep = ";",
  header = TRUE
)

# Kongruenz abkürzen und in vorherigen trials kodieren
data$cong_N1 <- dplyr::lag(data$cong) # lag explizit aus dem dplyr Paket
data$cong_N2 <- dplyr::lag(data$cong_N1) # lag explizit aus dem dplyr Paket

Im nächsten Schritt bereiten wir die Daten in einer Art und Weise auf, die typisch für experimentalpsychologische Studien mit einem Fokus auf RTs ist. Hierzu zählt zum einen eine grobe Ausreißerkorrektur, sowie der Ausschluss fehlerhafter Durchgänge. Außerdem entfernen wir die Übungsdurchgänge und die ersten beiden Trials pro Block (da erst ab Trial 3 die Kongruenz in Trial \(N-2\) betrachtet werden kann).⁵

# Filtern der Daten: Auschließen der Übungsblöcke, fehlerhafter Trials,
# und zu kurzer/langer RTs ("Ausreißer")
data <- subset(
  data,
  blocknr > 0 & trialnr > 2 & error == 0 & errortype == 0 &
    rt >= 200 & rt <= 1200
)

Wie wir bereits in Teil 1 gesehen haben, sind RTs (eher) nicht normalverteilt. Dies stellt zunächst ein Problem dar, da klassische LMMs normalverteilte Residuen annehmen. Wie schon hier wenden wir daher die übliche log-Transformation an:

data$logrt <- log(data$rt)

3.1 Eine erste Datenexploration

Bevor wir nun ein LMM schätzen, verschaffen wir uns zunächst einen Überblick über die Daten. Zuerst prüfen wir, ob überhaupt von allen Versuchspersonen mehrere Beobachtungen pro Bedingungskombination vorliegen. Andernfalls können wir ein LMM ggf. gar nicht schätzen.

H <- table( # abs. Häufigkeit der Zellen pro VP
  data$vpid, 
  data$cong,
  data$cong_N1,
  data$cong_N2
)
H <- as.vector(H) # Tabelle in einen Vektor umwandeln, um dann ....
range(H) # ... die enthaltenen Häufigkeiten zu analysieren

## [1] 14 41

Zu erkennen ist, dass pro Versuchsperson mindestens 14 Beobachtungen für die Kombination aus Kongruenz in Trial \(N\), Trial \(N-1\) und Trial \(N-2\) vorliegen. Nun kodieren wir noch die entsprechenden kategorialen Variablen als Faktoren. In einem varianzanalytischen Design mit mehreren Faktoren, ausgewertet als Regression, ist es dabei wichtig, dass die kategorialen Variablen effektkodiert werden. Andernfalls entsprechen die Koeffizienten nicht den Haupt- und Interaktionseffekten. Die Funktion mixed() aus dem Paket afex, die wir gleich verwenden werden, kodiert Faktoren auch per default mit einer Effektkodierung.⁶ Aus diesem Grund führen wir im folgenden Code-Chunk eine Effektkodierung explizit ein und umgehen auch so unnötige Warnungen und Informationen im Folgenden. Auch definieren wir unsere Versuchspersonennummern als Faktor um eine Warnung beim Aufruf von ezANOVA() zu vermeiden:

# contr.sum() setzen
data$cong <- as.factor(data$cong)
data$cong_N1 <- as.factor(data$cong_N1)
data$cong_N2 <- as.factor(data$cong_N2)

contrasts(data$cong) <- contr.sum(2)
contrasts(data$cong_N1) <- contr.sum(2)
contrasts(data$cong_N2) <- contr.sum(2)

# VP als Faktor um Warnungen zu vermeiden
data$vpid <- as.factor(data$vpid)

Als nächstes visualisieren wir die Daten und rechnen eine klassische mehrfaktorielle Varianzanalyse (ANOVA), um uns mit dem Einfluss der Kongruenz in Trial \(N\), Trial \(N-1\) und Trial \(N-2\) auf die RTs vertraut zu machen (siehe Abb. 3.1). Hierfür nutzen wir das Paket ez in Kombination mit dem Paket schoRsch:

library(ez)
# Einfacher Plot mit ez 
ezPlot(
  data = data,
  dv = rt,
  wid = vpid,
  within = .(cong, cong_N1, cong_N2),
  x = cong_N1,
  split = cong,
  col = cong_N2
)

Abbildung 3.1: Einfache Visualisierung der Simon-Daten in Abhängigkeit der Kongruenz in Trial \(N\), Trial \(N-1\) und Trial \(N-2\).

Zu erkennen ist, dass alle drei Faktoren (cong, cong_N1 und cong_N2) einen Einfluss auf die RT haben. Was diese Ergebnisse genau bedeuten, muss uns hier nicht wirklich interessieren, da wir aus didaktischen Gründen einen Fokus auf die anschließenden LMM Analysen haben. Dennoch sei angemerkt, dass wir ein sehr typisches Muster vorfinden: Nach einem inkongruenten Trial \(N-1\) (x-Achse im Plot) ist der Einfluss der Kongruenz im aktuellen Trial \(N\) (die vertikale Differenz zwischen blauen und roten Punkten) kleiner als nach einem kongruenten Trial \(N-1\). Dies wird als sequenzielle Modulation bezeichnet. Inhaltlich wird dieses Muster damit erklärt, dass nach einem inkongruenten Trial mehr kognitive Kontrolle verwendet wird, sodass der Einfluss der irrelevanten Information in der Simon-Aufgabe im nächsten Trial reduziert wird. Dadurch verschwindet der Kongruenzeffekt (oder er dreht sich manchmal sogar um).

Die entsprechende ANOVA gibt uns erste Hinweise über die Signifikanz der Haupt- und Interaktionseffekte:

library(schoRsch)
anova_out(ezANOVA(
  data = data,
  dv = rt,
  wid = vpid,
  within = .(cong, cong_N1, cong_N2),
  type = 3,
  detailed = TRUE
))

## $`--- ANOVA RESULTS     ------------------------------------`
##                 Effect        MSE df1 df2       F     p petasq getasq
## 1          (Intercept) 19488.1553   1  31 3598.48 0.000   0.99   0.99
## 2                 cong  1454.7822   1  31    8.34 0.007   0.21   0.02
## 3              cong_N1   502.8067   1  31    5.20 0.030   0.14   0.00
## 4              cong_N2   299.1362   1  31    5.08 0.031   0.14   0.00
## 5         cong:cong_N1   792.8648   1  31  169.84 0.000   0.85   0.15
## 6         cong:cong_N2   544.2194   1  31   11.83 0.002   0.28   0.01
## 7      cong_N1:cong_N2   397.6867   1  31    6.44 0.016   0.17   0.00
## 8 cong:cong_N1:cong_N2   543.4828   1  31    0.02 0.887   0.00   0.00
## 
## $`--- SPHERICITY TESTS  ------------------------------------`
## [1] "N/A"
## 
## $`--- FORMATTED RESULTS ------------------------------------`
##                 Effect                                    Text
## 1          (Intercept) F(1, 31) = 3598.48, p < .001, np2 = .99
## 2                 cong F(1, 31) =    8.34, p = .007, np2 = .21
## 3              cong_N1 F(1, 31) =    5.20, p = .030, np2 = .14
## 4              cong_N2 F(1, 31) =    5.08, p = .031, np2 = .14
## 5         cong:cong_N1 F(1, 31) =  169.84, p < .001, np2 = .85
## 6         cong:cong_N2 F(1, 31) =   11.83, p = .002, np2 = .28
## 7      cong_N1:cong_N2 F(1, 31) =    6.44, p = .016, np2 = .17
## 8 cong:cong_N1:cong_N2 F(1, 31) =    0.02, p = .887, np2 < .01
## 
## $`NOTE:`
## [1] "Reporting unadjusted p-values."

Tatsächlich sind hier einige Effekte signifikant: Wir finden einen signifikanten Haupteffekt von Trial \(N\) und Trial \(N-1\) (Variablen cong und cong_N1), allerdings (gerade so) keinen signifikanten Haupteffekt von Trial \(N-2\) (Variable cong_N2). Außerdem sind alle Zweifach-Interaktionen zwischen den Variablen signifikant, allerdings ist die Dreifachinteraktion nicht signifikant.

3.2 Auswertung mit einem LMM und Umgang mit fehlender Konvergenz

Zur Auswertung mit einem LMM nutzen wir hier die Funktion mixed() aus dem Paket afex (Singmann et al., 2024; siehe auch Singmann & Kellen, 2019). Diese Funktion stellt einen Wrapper um die Funktion lmer() dar und ist daher von der grundlegenden Funktionalität her identisch. Sie hat aber eine Reihe von Vorteilen, die v.a. bei der Auswertung von Experimentaldaten eines ANOVA-Designs nützlich sind:

1. mixed() prüft defaultmäßig Faktoren und kodiert diese ggf. im Sinne einer Kontrastkodierung (contr.sum). Eine solche Kodierung ist zwingend notwendig, damit die Koeffizienten des LMM konzeptuell zu den Haupt- und Interaktionseffekten im Sinne einer ANOVA passen.

2. mixed() erlaubt auf eine sehr einfache Art und Weise die Unterdrückung von Kovarianzen zwischen den random effects, wenn kategoriale Variablen im Modell verwendet werden. Prinzipiell geht dies auch mit lmer(), allerdings müssten hier händische Umkodierungen vorgenommen werden, um die korrekte Bestimmung von \(\bf Z_i\) zu gewährleisten. Die Funktion mixed() nimmt uns diese Umkodierung ab.

3. mixed() bietet Argumente, um leicht die verschiedenen Varianten der Signifikanztests der fixed-effects auszuwählen. Dies umfasst sowohl die Ansätze zur Berechnung der p-Werte, als auch die Auswahl des “Typs” (1 bis 3), mit denen Haupt- und Interaktioneffekte getestet werden. Per Default werden Prädiktoren bzw. die jeweiligen fixed-effects mit Hilfe der approximierten Freiheitsgrade nach Satterthwaite sowie Typ 3 getestet.

Wie bereits im letzten Abschnitt erwähnt, treten \(\alpha\)-Fehler beim Testen der fixed-effects vermehrt auf, wenn Variabilität in den random-effects nicht berücksichtigt wird, obwohl sie eigentlich vorhanden ist. Aus diesem Grund wird mitunter empfohlen, alle random-effects initial aufzunehmen, was dann als full random-effect Struktur bezeichnet wird (Barr et al., 2013). Gleichzeitig hat aber eine komplexe random-effect Struktur den Nachteil, dass diese nur schwer schätzbar ist. Aus diesem Grund erhalten wir bei der Berechnung eines LMMs mit vielen random-effects auch oft Warnungen. Um Konvergenzprobleme zu minimieren, setzen wir im folgenden Aufruf die Anzahl der Iterationen des Suchalgorithmus mit Hilfe des lmerControl Arguments höher, damit die Suche nicht vorzeitig abbricht. Allerdings können auch andere Warnungen auftreten, wie wir gleich sehen werden.

Wir beginnen nun mit dem Modell mit der maximalen random-effect Struktur und geben uns das Ergebnis kompakt mit der Funktion print() aus:

m1 <- mixed(logrt ~ cong * cong_N1 * cong_N2 + (cong * cong_N1 * cong_N2 | vpid),
  data = data,
  control = # Anzahl der Iterationen zur Optimierung hochstellen
    lmerControl(optCtrl = list(maxfun = 1e6))
)
print(m1)

## Warning: lme4 reported (at least) the following warnings for 'full':
##   * boundary (singular) fit: see help('isSingular')

## Mixed Model Anova Table (Type 3 tests, S-method)
## 
## Model: logrt ~ cong * cong_N1 * cong_N2 + (cong * cong_N1 * cong_N2 | 
## Model:     vpid)
## Data: data
##                 Effect       df          F p.value
## 1                 cong 1, 30.38  13.74 ***   <.001
## 2              cong_N1 1, 40.71     6.74 *    .013
## 3              cong_N2 1, 91.99     4.89 *    .029
## 4         cong:cong_N1 1, 31.36 224.26 ***   <.001
## 5         cong:cong_N2 1, 35.60  16.16 ***   <.001
## 6      cong_N1:cong_N2 1, 50.73     4.88 *    .032
## 7 cong:cong_N1:cong_N2 1, 42.41       0.01    .943
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Die Warnung besagt, dass das Modell einen singulären fit ergeben hat (boundary (singular) fit). Genauer gibt uns dies an, dass die Kovarianzmatrix der random-effects singulär ist und damit keine Inverse berechnet werden kann. Geometrisch kann man sich vorstellen, dass somit das Ellipsoid der random-effects, welches durch die ensprechende Kovarianzmatrix bestimmt wird, in einer Dimension die Länge 0 besitzt (vgl. auch Kapitel 11.2 aus Statistik 2). In einfachen Worten: Es konnte keine Kovarianzmatrix der random-effects bestimmt werden, die sinnvoll ist, da die Parametrisierung der random-effects (und der Residuen; wenn man nicht gerade lmer() nutzt) nicht beliebig ausfallen kann. Am Ende der Schätzung muss \(\bf {\mathcal{\hat D}}\) (oder eben auch \(\bf {\mathcal {\hat R_i}}\)) eine sinnvolle Kovarianzmatrix ergeben; andernfalls leiden sowohl die inhaltliche Interpretation des Modells als auch die hieraus abgeleiteten Inferenzen.

Im Allgemeinen gibt es keinen Konsens darüber, wie man mit singulären und nicht-konvergierenden LMMs umgehen soll (vgl. auch Barr et al., 2013; Matuschek et al., 2017):

• Im Falle einer Konvergenzwarnung (weniger einer Warnung über Singularität) hilft es die Anzahl der Iterationen des Such- bzw. Minimierungsalgorithmus hochzustellen (vgl. auch das definierte Argument maxfun des control Arguments im vorherigen Code-Beispiel) bzw. den Algorithmus selbst zu wechseln (indem man das Argument optimizer im Zuge des control Arugments anpasst; siehe help("lmerControl")).

• Insofern man sich hiermit auskennt, kann es auch helfen ein LMM bayesianisch zu schätzen (siehe bspw. das R-Paket brms; Bürkner, 2017).

• Bleibt man bei einer nicht-bayesianischen Schätzung, muss allerdings in den meisten Fällen die Komplexität der random-effects reduziert werden. Da bei LMMs ein primärer Fokus auf den fixed-effects liegt, folgen wir hier den eher pragmatischen Vorschlägen von Singmann und Kellen (2019) bzw. basieren alle folgenden Erklärungen auf der “Anleitung” zur Analyse von RT-Daten im Paket afex (siehe vignette(package = "afex", topic = "afex_mixed_example")). Die Idee ist dabei, mit dem komplexesten Modell anzufangen und dann sukzessive das Modell zu reduzieren, bis das Modell ohne Warnung konvergiert.

Als ersten Schritt können wir die Korrelationen bzw. Kovarianzen zwischen den random-effects aus dem Modell nehmen, was bereits eine große Anzahl an Parametern reduziert (man bedenke, dass es ggf. viele Elemente außerhalb der Diagonalen von \(\bf {\mathcal {\hat D}}\) gibt und diese i.d.R. ohne Einschränkungen geschätzt werden). Hier hilft nun die Funktion mixed(), denn sie gewährleistet eine korrekte Schätzung unkorrelierter random-effects auch beim Vorhandensein kategorialer Variablen mit Hilfe der ||-Schreibweise. Hierfür muss aber als zusätzliches Argument expand_re = TRUE gesetzt werden:

m2 <- mixed(logrt ~ cong * cong_N1 * cong_N2 + (cong * cong_N1 * cong_N2 || vpid),
            data = data,
            control = # Anzahl der Iterationen zur Optimierung hochstellen
              lmerControl(optCtrl = list(maxfun = 1e6)),
            expand_re = TRUE
)
print(m2)

## Warning: lme4 reported (at least) the following warnings for 'full':
##   * boundary (singular) fit: see help('isSingular')

## Mixed Model Anova Table (Type 3 tests, S-method)
## 
## Model: logrt ~ cong * cong_N1 * cong_N2 + (cong * cong_N1 * cong_N2 || 
## Model:     vpid)
## Data: data
##                 Effect         df          F p.value
## 1                 cong   1, 30.33  14.27 ***   <.001
## 2              cong_N1   1, 29.15    7.68 **    .010
## 3              cong_N2 1, 6783.59     6.04 *    .014
## 4         cong:cong_N1   1, 29.54 231.22 ***   <.001
## 5         cong:cong_N2   1, 31.64  18.86 ***   <.001
## 6      cong_N1:cong_N2 1, 6649.86     6.61 *    .010
## 7 cong:cong_N1:cong_N2   1, 29.68       0.00   >.999
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Das geschätzte Modell besitzt weiterhin eine singuläre Kovarianzmatrix der random-effects. Daher muss in den nächsten Schritten weiter reduziert werden. Eine Hilfestellung bieten hierbei die Varianzen der random-effects. Dazu verwenden wir die Funktion VarCorr(), die (eigentlich) auf das Ergebnis der lmer()-Funktion angewendet wird. Da mixed() mehr zurück gibt als nur das Ergebnis der Funktion lmer(), muss dies hier explizit mit $full_model adressiert werden:

VarCorr(m2$full_model)

##  Groups   Name                              Std.Dev.  
##  vpid     (Intercept)                       8.6561e-02
##  vpid.1   re1.cong1                         2.2438e-02
##  vpid.2   re1.cong_N11                      6.0578e-03
##  vpid.3   re1.cong_N21                      0.0000e+00
##  vpid.4   re1.cong1_by_cong_N11             1.1104e-02
##  vpid.5   re1.cong1_by_cong_N21             7.4437e-03
##  vpid.6   re1.cong_N11_by_cong_N21          1.4884e-05
##  vpid.7   re1.cong1_by_cong_N11_by_cong_N21 3.9816e-03
##  Residual                                   1.7658e-01

Offenbar hat der random-effect für den Prädiktor cong_N2 eine geschätzte Varianz von \(0\).⁷ Da es allerdings meist wenig Sinn ergibt, individuelle Variabilität in höheren Interaktionen in einem Modell zu belassen, wenn die entsprechende Variable als Haupteffekt in den random-effects fehlt, wird empfohlen, zuerst höhere Interaktionen herauszunehmen, in die die jeweilige Variable involviert ist. Daher reduzieren wir die random-effect Struktur zunächst um die Dreifachinteraktion. Einen einfachen Weg dies zu erreichen, ist es R mit Hilfe einer ^2-Schreibweise anzuweisen, nur 2-fach Interaktionen zu berücksichtigen :

m3 <- mixed(logrt ~ cong * cong_N1 * cong_N2 + ( (cong + cong_N1 + cong_N2)^2 || vpid),
            data = data,
            control = # Anzahl der Iterationen zur Optimierung hochstellen
              lmerControl(optCtrl = list(maxfun = 1e6)),
            expand_re = TRUE
)
print(m3)

## Warning: lme4 reported (at least) the following warnings for 'full':
##   * boundary (singular) fit: see help('isSingular')

## Mixed Model Anova Table (Type 3 tests, S-method)
## 
## Model: logrt ~ cong * cong_N1 * cong_N2 + ((cong + cong_N1 + cong_N2)^2 || 
## Model:     vpid)
## Data: data
##                 Effect         df          F p.value
## 1                 cong   1, 30.35  14.32 ***   <.001
## 2              cong_N1   1, 29.14    7.69 **    .010
## 3              cong_N2 1, 6792.19     6.08 *    .014
## 4         cong:cong_N1   1, 29.54 231.52 ***   <.001
## 5         cong:cong_N2   1, 31.64  18.90 ***   <.001
## 6      cong_N1:cong_N2 1, 6785.83     6.61 *    .010
## 7 cong:cong_N1:cong_N2 1, 6765.66       0.00    .999
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Auch dieses Modell liefert weiterhin eine Singularitätswarnung und auch hier besitzt der random-effect für den Haupteffekt von cong_N2 eine Varianz von 0:

VarCorr(m3$full_model)

##  Groups   Name                     Std.Dev. 
##  vpid     (Intercept)              0.0865525
##  vpid.1   re1.cong1                0.0223911
##  vpid.2   re1.cong_N11             0.0060633
##  vpid.3   re1.cong_N21             0.0000000
##  vpid.4   re1.cong1_by_cong_N11    0.0110876
##  vpid.5   re1.cong1_by_cong_N21    0.0074308
##  vpid.6   re1.cong_N11_by_cong_N21 0.0000000
##  Residual                          0.1766267

Wir entfernen daher nun alle Zweifach-Interaktionen aus den random-effects, in denen die entsprechende Variable vorkommt:

m4 <- mixed(logrt ~ cong * cong_N1 * cong_N2 + ( cong + cong_N1 + cong_N2 || vpid),
            data = data,
            control = # Anzahl der Iterationen zur Optimierung hochstellen
              lmerControl(optCtrl = list(maxfun = 1e6)),
            expand_re = TRUE
)
print(m4)
VarCorr(m4$full_model)

## Warning: lme4 reported (at least) the following warnings for 'full':
##   * boundary (singular) fit: see help('isSingular')

## Mixed Model Anova Table (Type 3 tests, S-method)
## 
## Model: logrt ~ cong * cong_N1 * cong_N2 + (cong + cong_N1 + cong_N2 || 
## Model:     vpid)
## Data: data
##                 Effect         df          F p.value
## 1                 cong   1, 30.35  14.38 ***   <.001
## 2              cong_N1   1, 29.08    7.69 **    .010
## 3              cong_N2 1, 6817.74     6.63 *    .010
## 4         cong:cong_N1 1, 6792.95 422.11 ***   <.001
## 5         cong:cong_N2 1, 6807.29  26.36 ***   <.001
## 6      cong_N1:cong_N2 1, 6794.54    6.65 **    .010
## 7 cong:cong_N1:cong_N2 1, 6821.31       0.00    .979
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

##  Groups   Name         Std.Dev. 
##  vpid     (Intercept)  0.0864558
##  vpid.1   re1.cong1    0.0223850
##  vpid.2   re1.cong_N11 0.0063597
##  vpid.3   re1.cong_N21 0.0000000
##  Residual              0.1770990

Die Warnung ist weiterhin vorhanden und die geschätzte Varianz des random-effects für den Haupteffekt cong_N2 ist immer noch \(0\). Wir entfernen daher nun auch den Haupteffekt als random-effect:

m5 <- mixed(logrt ~ cong * cong_N1 * cong_N2 + ( cong + cong_N1 || vpid),
            data = data,
            control = # Anzahl der Iterationen zur Optimierung hochstellen
              lmerControl(optCtrl = list(maxfun = 1e6)),
            expand_re = TRUE
)
VarCorr(m5$full_model)

print(m5)

## Mixed Model Anova Table (Type 3 tests, S-method)
## 
## Model: logrt ~ cong * cong_N1 * cong_N2 + (cong + cong_N1 || vpid)
## Data: data
##                 Effect         df          F p.value
## 1                 cong   1, 30.35  14.38 ***   <.001
## 2              cong_N1   1, 29.08    7.69 **    .010
## 3              cong_N2 1, 6817.73     6.63 *    .010
## 4         cong:cong_N1 1, 6792.94 422.11 ***   <.001
## 5         cong:cong_N2 1, 6807.29  26.36 ***   <.001
## 6      cong_N1:cong_N2 1, 6794.54    6.65 **    .010
## 7 cong:cong_N1:cong_N2 1, 6821.31       0.00    .979
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

VarCorr(m5$full_model)

##  Groups   Name         Std.Dev. 
##  vpid     (Intercept)  0.0864526
##  vpid.1   re1.cong1    0.0223844
##  vpid.2   re1.cong_N11 0.0063586
##  Residual              0.1770991

Dieses Mal erhalten wir keine Warnung mehr und das Modell konvergiert scheinbar ohne Probleme. Hierzu passt natürlich, dass alle verbleibenden random-effects eine Varianz \(>0\) aufweisen. Als letzten Schritt führen wir nun wieder Korrelationen zwischen den random-effects ein, um zu prüfen, ob das Modell noch etwas komplexer gestaltet werden kann:

m6 <- mixed(logrt ~ cong * cong_N1 * cong_N2 + ( cong + cong_N1 | vpid),
            data = data,
            control = # Anzahl der Iterationen zur Optimierung hochstellen
              lmerControl(optCtrl = list(maxfun = 1e6)),
            expand_re = TRUE
)
print(m6)

## Warning: lme4 reported (at least) the following warnings for 'full':
##   * Model failed to converge with max|grad| = 0.0022802 (tol = 0.002, component 1)

## Mixed Model Anova Table (Type 3 tests, S-method)
## 
## Model: logrt ~ cong * cong_N1 * cong_N2 + (cong + cong_N1 | vpid)
## Data: data
##                 Effect         df          F p.value
## 1                 cong   1, 30.35  14.37 ***   <.001
## 2              cong_N1   1, 29.10    7.86 **    .009
## 3              cong_N2 1, 6820.07    6.74 **    .009
## 4         cong:cong_N1 1, 6793.79 421.38 ***   <.001
## 5         cong:cong_N2 1, 6807.66  26.60 ***   <.001
## 6      cong_N1:cong_N2 1, 6796.84    6.75 **    .009
## 7 cong:cong_N1:cong_N2 1, 6824.73       0.00    .968
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Allerdings erhalten wir hier jetzt eine Konvergenzwarnung. Das komplexeste LMM, welches wir problemlos verwenden können, ist also Modell 5 und damit haben wir den Modellselektionsprozess bzgl. der random-effects beendet (zumindest vorläufig, denn natürlich sollten wir im nächsten Schritt die Residuen etc. ansehen und die Daten auf einflussreiche Datenpunkte überprüfen).

Bevor nun im weiteren Schritt die Ergebnisse der Signifikanztests der fixed-effects berichtet werden, erinnern wir uns, dass fixed-effects leichter (fälschlicherweise) signifikant werden, wenn die entsprechenden random-effects in Wahrheit Variabilität aufweisen (aber nicht mehr im Modell sind). Aus diesem Grund lohnt es sich, die Ergebnisse aus dem komplexesten Modell (Modell 1) mit den Ergebnissen des finalen Modells (Modell 5) zu vergleichen. Auch wenn (aus statistischer Sicht) das komplexeste Modell vielleicht nicht adäquat ist, kann es eine sinnvolle Referenz liefern (die \(p\)-Werte wurden auf Basis der Methode von Satterthwaite berechnet; ggf. möchte man hier den Funktionsaufruf von mixed() modifizieren um bspw. die Methode nach Kenward und Roger, 1997, auszuwählen). Natürlich können wir die beiden Modelle nacheinander nocheinmal anschauen:

print(m1)
print(m5)

Schöner ist es allerdings (auch wenn es hier platzmäßig nicht ganz aufgeht), die Ergebnisse beider Modelle direkt nebeneinander zu platzieren. Dies kann z.B. mit der Funktion left_join() aus dem dplyr Paket erfolgen:

dplyr::left_join(nice(m1), nice(m5),
  by = "Effect",
  suffix = c("_full", "_final")
)

Vergleicht man die Inferenzen, erkennt man, dass bzgl. der eigentlichen Signifikanz beide Modelle zur gleichen Konklusion kommen. Im Übrigen ergeben sich auch nahezu die gleichen Ergebnisse wie in der klassischen ANOVA (siehe weiter oben den Aufruf von ezANOVA()).

4 Referenzen

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.

Bürkner, P.-C. (2017). brms: An R package for bayesian multilevel models using Stan. Journal of Statistical Software, 80(1), 1-28.

Claeskens, G. (2013). Lack of fit, graphics, and multilevel model diagnostics. In M. Scott, J. Simonoff, & B. Marx(Eds.). The Sage handbook of multilevel mModeling, (pp. 425-443). Sage Publications.

Davis, C., Blackmore, E., Katzman, D. K., & Fox, J. (2005). Female adolescents with anorexia nervosa and their parents: A case-control study of exercise attitudes and behaviours. Psychological Medicine, 35(3), 377-386.

Fox, J. (2015). Applied regression analysis and generalized linear models. Sage Publications.

Galecki, A. & Burzykowski, T. (2013). Linear mixed-effects model using R. Springer New York.

Kenward, M.G. & Roger, J. (1997). Small sample inference for fixed effects from restricted maximum likelihood. Biometrics, 53, 983-997.

Matuschek, H., Kliegl, R., Vasishth, S., Baayen, H. & Bates, D. (2017). Balancing type I error and power in linear mixed models. Journal of Memory and Language, 94, 305–315.

Nieuwenhuis, R., te Grotenhuis, M. & Pelzer, B. (2012). influence.ME: Tools for detecting influential data in mixed effects models. R Journal, 4(2), 38-47.

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.

Verbeke, G. & Molenberghs, G. (2000). Linear mixed models for longitudinal data. Springer.

Waterman, M. J., Birch, J. B. & Abdel-Salam, A.-S. G. (2015). Several nonparametric and semiparametric approaches to linear mixed model regression. Journal of Statistical Computation and Simulation, 85(5), 956-977.

West, B. T., Welch, K. B. & Galecki, A. T. (2022). Linear mixed models: A practical guide using statistical software. Chapman and Hall.

---

1. In der Syntax wird ein random intercept immer implizit mit berücksichtigt, sobald man einen random slope schätzt. Die explizite Schreibweise mit am Anfang des Skripts hatten wir nur aus didaktischen Gründen verwendet. Möchte man einen random intercept explizit nicht schätzen, wohl aber einen random slope, dann kann man die Schreibweise nutzen.↩︎

2. Die REML-Schätzung ist hierbei der default und kann über beim Aufruf von explizit ein- bzw. umgestellt werden.↩︎

3. In der Modellsyntax tritt die Schreibweise auf. Der “Punkt” ist hierbei ein Platzhalter und R setzt für ihn das Objekt ein. Die Schreibweise ist also “nur” eine Abkürzung, um den Output der Funktion und nicht vorher selbst explizit abzuspeichern.↩︎

4. DFBETAs funktionieren hierbei analog zur Cook’s Distance, beziehen sich aber immer nur auf einen bestimmten Parameter. Sie berechnen sich als standardisierte Differenz eines Parameters zwischen den Modellen mit und ohne einer Beobachtung bzw. Gruppe.↩︎

5. Häufig werden noch weitere Maßnahmen zur Datenbereinigung ergriffen. Beispielweise ist es auch üblich, Trials, die einem Fehler folgen, auszuschließen. Hier haben wir es aber bei einer oberflächlichen Bereinigung belassen.↩︎

6. Dies zu wissen ist besonders wichtig, wenn man ein Design wie in Kapitel 2 mit der Funktion auswerten würde. Durch die automatische, aber dann unnötige, Umkodierung bedeuten die Koeffizienten etwas anderes als man typischerweise anstrebt und auch die Signifikanzwerte ändern sich. Das Paket dient aber auch vornehmlich zur Analyse faktorieller Experimente und der Name steht für “analysis of factorial experiments”.↩︎

7. Kleine Standardabweichungen deuten natürlich auch auf wenig interindividuelle Variabilität hin (dabei ist aber auch auf die Skalierung der Variablen immer zu achten). Ein Wert von exakt 0 ist aber eher als Hinweis auf Schätzprobleme zu verstehen.↩︎