Mehrebenenanalyse (HLM) mit R:
4. Random Effects auf Signifikanz prüfen

Arndt Regorz, Dipl. Kfm. & MSc. Psychologie, 19.05.2022

Wenn Sie im Rahmen einer Mehrebenenanalyse (oder HLM oder linear mixed effects model – im Folgenden werden diese Begriffe austauschbar verwendet) Random Effects berücksichtigen, möchten Sie in der Regel auch prüfen, ob sich diese signifikant von Null unterscheiden. Dieses Tutorial zeigt am Beispiel des lme4-Moduls von R, wie das funktioniert, wobei das Vorgehen im Prinzip auch auf andere HLM-Packages in R übertragbar ist.

Die Berechnungen beruhen auf einem Beispieldatensatz von Hox für das 2. Kapitel des Buchs von Hox et al. (2017). Der Datensatz kann bei Hox (http://joophox.net/mlbook2/DataExchange.zip) heruntergeladen werden.

INHALT

  1. Grundlagen
  2. Beispiel Signifikanztest Random Slope in R
  3. Vollständiger R-Code
  4. Literatur

1. Grundlagen

Für die Fixed Effects erhält man in lme4 bei der zusätzlichen Verwendung des Moduls lmerTest automatisch p-values. Aber wie bekommt man p-values für Random Effects in R?

Im hier vorgestellten Ansatz wird im Rahmen des Hierarchical Linear Modeling die Deviance zwischen verschiedenen genesteten Modellen miteinander verglichen. Wie im Grundlagentutorial geschildert, ist die Deviance ein relatives Maß für die Modellgüte, die man zum Vergleich von genesteten Modellen verwenden kann. Die Differenz der Deviance folgt einer Chi-Quadrat-Verteilung mit so vielen Freiheitsgraden, wie in dem umfassenderen Modell zusätzliche Parameter geschätzt wurden.

Da die Deviance im lme4-Output standardmäßig mit ausgegeben wird, so dass man einfach die Differenz in der Deviance zwischen zwei Modellen berechnen kann, erscheint dies erst einmal vergleichsweise simpel. Allerdings gibt es in der Praxis zwei Komplikationen:

Zum einen erfolgt in der Regel bei der Hinzunahme einer Random Slope zusätzlich auch die Schätzung (mindestens) einer zusätzlichen Kovarianz, so dass mehr als ein Parameter hinzu kommt. Ein entsprechender Test mit Hilfe der Deviance testet dann nicht mehr die Signifikanz der Random Slope alleine, sondern die Signifikanz der gemeinsamen Hereinnahme von Random Slope und Kovarianz. Diese beiden Effekte sollte man also noch auseinander rechnen.

Zum anderen gibt es bei einem Chi-Quadrat-Test für einen Varianzparameter eine weitere Besonderheit: Eine Varianz kann nicht kleiner werden als 0. Damit sind, anders als sonst bei einem Chi-Quadrat-Test, Abweichungen nur in eine Richtung möglich. Um dies zu kompensieren, sollte der resultierende p-Wert halbiert werden (siehe Kapitel 3 aus Hox et al., 2017). Dies gilt aber nur für die Varianz und nicht für eine Kovarianz, weil bei der Kovarianz Abweichungen von 0 in beide Richtungen möglich sind, so dass man dort den p-Wert nicht halbieren darf.

Das ist alles noch sehr theoretisch und abstrakt - im folgenden Beispiel sehen Sie ganz konkret, wie man mit diesen beiden Problemen umgehen kann.

2. Beispiel Signifikanztest Random Slope in R

Dieses Beispiel beruht auf den Modellschätzungen aus dem Tutorial Mehrebenenanalyse mit R (HLM mit R): lme4 Package, dort wird der jeweilige Code ausführlich erklärt. Als Packages wurden im Folgenden verwendet:

library(lme4)
library(lmerTest)

Nachdem dort die Schätzung einer Random Slope für das Geschlecht zu Schätzproblemen geführt hatte, wird hier nur eine Random Slope für die Extraversion der Schülerinnen und Schüler geschätzt. Für die Auswirkungen dieser Random Slope vergleichen wir zunächst ein Modell mit Fixed Slopes (auf Level 1 für Geschlecht und Extraversion der Schülerinnen und Schüler sowie als Level 2 Prädiktor Berufserfahrung der Lehrkraft) mit einem Modell mit zusätzlich einer Random Slope für die Extraversion der Schülerinnen und Schüler.

Die Aufrufe der beiden Modelle waren:

# Fixed Slopes (Random intercept)

fixed_slopes <- lmer(popular ~ 1 + sex + extrav + texp +
(1|class), data=popular2, REML=F)
summary(fixed_slopes)

# Random Slopes nur Extraversion

random_slopes2 <- lmer(popular ~ 1 + sex + extrav + texp +
(1+extrav|class), data=popular2, REML=F)
summary(random_slopes2)

Output Modell Fixed Slopes:

Grafik Output Fixed Slope


Output Modell Random Slope für Extraversion:

Grafik Output Random Slope


Für unsere Fragestellung sind die wesentlichen Informationen:

(A) Deviance des Modells mit Fixed Slopes 4862.3

(B) Deviance des Modells mit Random Slope Extraversion 4812.8

(C) Es wurde zusätzlich zur Random Slope eine Random Covariance geschätzt.

Die Differenz in der Deviance, 4862.3 – 4812.8 = 49.5 lässt sich jetzt schwer interpretieren, weil neben der Random Slope noch eine Kovarianz geschätzt wurde (zwischen Random Slope und Random Intercept) und die Differenz in der Deviance beides betrifft, Random Slope Extraversion und Kovarianz zwischen Random Slope und Random Intercept.

Die Lösung für dieses Problem ist, einen Zwischenschritt einzubauen, bei dem zunächst nur die Random Slope geschätzt wird, die zusätzliche Random Covariance zwischen Random Intercept und Random Slope jedoch auf Null gesetzt wird. Das ist mit dem folgenden R-Code möglich:

random_slopes3 <- lmer(popular ~ 1 + sex + extrav + texp + (1|class)
+ (0+extrav|class), data=popular2, REML=F)
summary(random_slopes3)

Geändert gegenüber dem vorherigen Modell hat sich nur die Schätzung der Random Effects:

(1|class) + (0+extrav|class)

Die erste Klammer ist die Schätzung des Random Intercept, und die zweite Klammer die Schätzung der Random Slope ohne Intercept (durch die 0 +). Und wenn man diese beiden Random Effects so aufruft, wird keine Kovarianz der Random Effects berechnet, weil in jeder der beiden Klammern nur ein Random Effect steht (Random Intercept in der ersten Klammer, Random Slope in der zweiten Klammer).

Als Output für dieses Modell erhalten wir.

Grafik Output nur Random Slope ohne Random Covariance


(D) Deviance des Modells mit Random Slope, ohne Kovarianz 4852.0

(E) Hier ist jetzt keine Schätzung einer Kovarianz der Random Effects mehr enthalten.

Jetzt haben wir eine Abfolge von drei Modellen, bei denen jeweils nur ein Zufallsparameter hinzugekommen ist:

Schritt 1: Nur Fixed Effects und Random Intercept, Deviance = 4862.3

Schritt 2: Random Slope Extraversion, Deviance = 4852.0
Differenz zum vorherigen Schritt: 4862.3 – 4852 = 10.3 für einen zusätzlichen Parameter (Random Slope)

Schritt 3: Random Slope und Random Covariance, Deviance = 4812.8
Differenz zum vorherigen Schritt: 4852.0 – 4812.8 = 39.2 für einen zusätzlichen Parameter (Covariance)

Diese beiden Differenzen sind jeweils Chi-Quadrat-verteilt mit jeweils einem Freiheitsgrad, wobei wir bei der ersten Differenz (für die Varianzschätzung der Random Slope) noch den p-Wert halbieren müssen. Die p-values können wir einfach mit R berechnen:

# Varianz Random Slope
# p-Wert, halbiert wegen einseitiger Test bei Varianz
0.5 * pchisq(10.3, df=1, lower.tail = F)

Hier ergibt sich ein p-Wert = .000665151, wir haben also eine signifikante Varianz in der Slope für Extraversion - und das war es, was wir primär wissen wollten.

# Kovarianz Random Intercept und Random Slope
pchisq(39.2, df=1, lower.tail = F)

Auch hier ergibt sich ein signifikanter p-Wert = 3.825402e-10, wir haben also auch eine signikante Kovarianz der beiden Zufallsparameter (was allerdings nicht ganz so interessant ist).

Alternativ kann man auch das Modul varTestnlme nutzen, um mit Hilfe der Fit-Objekte die Varianzkomponenten auf Signifikanz zu testen. Auch hier empfiehlt sich der Zwischenschritt eines Modells mit nur Random Slope, aber ohne Kovarianz zwischen den Random Effects.

# Alternativ mit Modul varTestnlme (ohne manuelle Differenz der Deviance)
library(varTestnlme)

# Test Random Slope (einseitiger Test für die Varianz automatisch berücksichtigt)
summary(varCompTest(random_slopes3, fixed_slopes))

# Test Random Covariance
summary(varCompTest(random_slopes2, random_slopes3))

Für die Random Slope ergibt sich hier p = .0006499029.

Für die Covariance ergibt sich p = 3.933221e-10.

Sie sehen, dass (bis auf Rundungsdifferenzen) die gleichen Ergebnisse wie oben herauskommen. Dabei ist durch das varTestnlme-Modul für die Random Slope bereits automatisch die Halbierung des p-values vorgenommen worden.

3. Vollständiger R-Code

library(lme4)
library(lmerTest)
head(popular2, 5)

# 2. Fixed Slopes (Random intercept)

fixed_slopes <- lmer(popular ~ 1 + sex + extrav + texp + (1|class), data=popular2, REML=F)
summary(fixed_slopes)

# 3.b Random Slopes nur Extraversion

random_slopes2 <- lmer(popular ~ 1 + sex + extrav + texp + (1+extrav|class),
data=popular2, REML=F)
summary(random_slopes2)

# 3.c Random Slopes nur Extraversion, keine Covarianz Intercept-Slope

random_slopes3 <- lmer(popular ~ 1 + sex + extrav + texp + (1|class) +
(0+extrav|class), data=popular2, REML=F)
summary(random_slopes3)

# Berechnung der p-Werte

# Varianz Random Slope
# p-Wert halbiert wegen einseitiger Test bei Varianz
0.5 * pchisq(10.3, df=1, lower.tail = F)

# Kovarianz Random Intercept und Random Slope
pchisq(39.2, df=1, lower.tail = F)

# Alternativ mit Modul varTestnlme (ohne manuelle Differenz der Deviance)
library(varTestnlme)

# Test Random Slope (einseitiger Test automatisch berücksichtigt)
summary(varCompTest(random_slopes3, fixed_slopes))

# Test Random Covariance
summary(varCompTest(random_slopes2, random_slopes3))

4. Literatur

Hox, J. J., Moerbeek, M., & Van de Schoot, R. (2017). Multilevel analysis: Techniques and applications (3rd edition). Routledge.


Weitere Tutorials zur Mehrebenenanalyse/HLM mit R: