Mehrebenenanalyse mit R (HLM mit R)
– Robuste Schätzung

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

Mehrebenenanalysen (=Hierarchische lineare Modelle, HLM) haben Voraussetzungen, und wenn diese verletzt sind, können verzerrte Schätzungen resultieren. Ein Weg, um damit umzugehen, ist die Verwendung von robusten Schätzverfahren. Deren Stärke liegt insbesondere im Umgang mit Ausreißern.

Dieses Tutorial zeigt Ihnen, wie Sie mit dem Modul robustlmm eine robuste Schätzung für ein HLM-Modell durchführen können.

INHALT

1. Robuste Schätzung eines Mehrebenenmodells
2. Vollständiger R-Code
3. Literatur

1. Robuste Schätzung eines Mehrebenenmodells

Das Beispiel beruht auf den Daten aus dem 2. Kapitel von Hox et al. (2017). Wir wollen dabei die Popularität von Schülerinnen und Schülern vorhersagen durch deren Geschlecht und deren Extraversion (Level 1 Prädiktoren, also auf Ebene der einzelnen Personen) sowie durch die Berufserfahrung der Lehrkraft (Level 2 Prädiktor, also auf Ebene der Schulklasse).

Zunächst laden wir die nötigen Packages (vor erstmaliger Benutzung müssen diese natürlich einmalig mit install.packages installiert worden sein):

library(lme4)
library(robustlmm)

Bevor wir die robuste Schätzung durchführen, schätzen wir einmal das Modell mit dem gewöhnlichen HLM-Schätzer aus dem lme4-Modul. Dabei verwenden wir für dieses Beispiel ein Modell mit Random Intercept und einer Random Slope für Extraversion. Ausführlichere Informationen zur gewöhnlichen Schätzung finden Sie in meinem Basistutorial zur Mehrebenenanalyse mit dem lme4-Modul. Dort finden Sie auch einen Link zu den für dieses Tutorial verwendeten Daten von Hox.

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

Als Ergebnis erhalten wir:

Vom Ergebnis sind im Vergleich zu den späteren robusten Schätzungen besonders interessant die Varianzen der Level 2 Effekte:

(A) Varianz des Random Intercepts, 1.28

(B) Varianz der Random Slope für Exraversion, 0.03

Jetzt führen wir die robuste Schätzung durch. Der Aufruf ist fast der gleiche wie bei der gewöhnlichen Schätzung, nur dass wir statt der Funktion lmer aus dem lme4-Modul die Funktion rlmer aus dem robustlmm-Modul verwenden:

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

Als Ergebnis erhalten wir:

Die Schätzung dauert mit diesen Beispieldaten sehr lange (je nach Computer können das 5-15 Minuten sein). Das deutet schon auf Schätzprobleme hin.

(C) Diese Warnmeldung zeigt, dass mit der Schätzung etwas nicht stimmt.

In so einem Fall sollte man sich zunächst die Varianzen der Random Effects ansehen, da häufig eine zu geringe Varianz von Random Effects für Schätzprobleme im Rahmen von HLM verantwortlich ist. So ein Problem hatten wir für dieses Modell (Random Intercept sowie Random Slope für Extraversion) mit der gewöhnlichen HLM-Schätzung nicht. Warum tritt dies jetzt bei der robusten Schätzung auf?

Die robuste Schätzung gewichtet einige der Beobachtungen (Level 1 und Level 2) geringer, um den Effekt von Ausreißern zu reduzieren. Und wenn Ausreißer geringer gewichtet werden, sinkt damit die Varianz der Parameterschätzungen, da es die Ausreißer sind, die relativ am stärksten zur Varianz beitragen.

(D) Die Varianz des Random Intercepts beträgt 0.39, hier liegt das Problem vermutlich nicht.

(E) Die Varianz der Random Slope ist sehr gering, 0.0009, und damit nur bei ca. 3% der vergleichbaren Varianz für die Random Slope aus dem gewöhnlichen HLM-Modell (0.03, siehe oben), die auch schon nicht sehr hoch war. Vermutlich ist das die Ursache für unsere Schätzprobleme.

Daher schätzen wir jetzt das gleiche Modell, jedoch ohne die Random Slope für die Extraversion.

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

Als Ergebnis erhalten wir:

Die Schätzung erfolgt sehr schnell und ohne Fehlermeldungen. Im oberen Teil des Output finden wir die üblichen Informationen wie bei anderen Mehrebenenmodellen auch.

(F) Es sieht auf den ersten Blick so aus, dass die Varianz des Random Intercept durch die robuste Schätzung deutlich niedriger ist im Vergleich zur gewöhnlichen Schätzung. Dabei müssen wir aber berücksichtigen, dass wir jetzt auch ein anderes Modell schätzen (ohne Random Slope).

Im unteren Bereich des Output erhalten wir Detailinformationen zur robusten Schätzung.

(G) Robustheits-Gewichte für Level 1. Wir sehen beispielsweise, dass der Fall mit der geringsten Gewichtung in der Schätzung nur ein Gewicht von 0.410 erhalten hat, und dass 400 von 2000 Fällen weniger als mit dem vollen Gewicht in die Schätzung eingegangen sind.

(H) Robustheits-Gewichte für Level 2 Wir sehen hier beispielsweise, dass die Klasse mit der geringsten Gewichtung in der Schätzung nur ein Gewicht von 0.401 erhalten hat, und dass 17 von 100 Klassen nicht mit dem vollen Gewicht in die Schätzung eingegangen sind.

Diese Werte deuten tatsächlich darauf hin, dass wir so gut wie keine extremen Ausreißer im Datensatz haben. Im Rahmen z.B. der robusten Regression habe ich auch schon Schätzungen gesehen, bei denen einige Fälle (extreme Ausreißer) fast auf 0.00 heruntergewichtet worden sind. Eine so starke Korrektur war im vorliegenden Datensatz offenbar nicht nötig, was auch zu den Ergebnissen der Voraussetzungsprüfung für diesen Datensatz passt (siehe Tutorial Mehrebenenanalyse mit R: Voraussetzungen prüfen).

Wenn Sie eine konkrete Liste erhalten möchten, welche Fälle (Level 1) bzw. Gruppen (Level 2) wie gewichtet worden sind, können Sie je eine Liste der Gewichtungsfaktoren anfordern:

getME(rob_fixed_slopes, "w_e") # Level 1
getME(rob_fixed_slopes, "w_b") # Level 2

Die Listen sind zu groß, um sie hier sinnvoll abzudrucken. Sie finden dort jeweils die Fallnummer bzw. Gruppennummer, sowie das Gewicht (1.00 oder weniger). Es kann gerade bei den Level 2 Einheiten aufschlussreich sein, sich die Gruppen (in unserem Beispiel Klassen) mit deutlich geringeren Gewichten einmal daraufhin anzusehen, was dort ggf. besonders ist.

(I) Hier finden Sie technische Schätzparameter. Um diese zu verstehen (und ggf. auch zu ändern), müssten Sie wesentlich tiefer in die Technik und Mathematik der robusten Schätzung einsteigen.

Mit dem folgenden Aufruf können Sie die technische Dokumentation dazu anfordern:

vignette("rlmer")

Hinsichtlich der Stichprobengröße sollten Sie beachten, dass durch die Heruntergewichtung von Fällen bzw. Gruppen bei der Schätzung die Anforderungen an den Stichprobenumfang in der Regel noch höher sind als in der gewöhnlichen HLM-Schätzung.

2. Vollständiger R-Code

library(lme4)
library(robustlmm)

head(popular2, 5)

# Normale Modellschätzung, Random Slopes nur Extraversion

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

# Robuste Modellschätzung, Random Slopes nur Extraversion

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

# Robuste Modellschätzung, keine Random Slopes

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

# Gewichte

getME(rob_fixed_slopes, "w_e") # Level 1
getME(rob_fixed_slopes, "w_b") # Level 2

# Zusatzinformationen:

vignette("rlmer")

3. Literatur

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

Koller, M. (2016). robustlmm: An R package for robust estimation of linear mixed-effects models. Journal of Statistical Software,, 75(6), 1–24. https://doi.org/10.18637/jss.v075.i06

Weitere Tutorials zur Mehrebenenanalyse mit R

• Mehrebenenanalyse (mit R) 1: Grundlagen mit lme4
• Mehrebenenanalyse (mit R) 2: Voraussetzungen prüfen
• Mehrebenenanalyse (mit R) 3: Zentrierung
• Mehrebenenanalyse (mit R) 4: p-Values Random Effects ermitteln
• Mehrebenenanalyse (mit R) 5: Ergebnisse berichten
• Mehrebenenanalyse (mit R) 7: Grundlagen mit nlme
• Formeln Mehrebenenanalyse verstehen: Videotutorial mit graphischer Erklärung der Formeln einer Mehrebenenanalyse (nicht R-spezifisch)