## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## Multilevel-Modelle
## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#### ++ Import der Daten und wichtige Pakete laden ####
# setwd("")
daten <- read.table( "SnijdersBosker2012.txt", header = TRUE )
# ggf. Pakete installieren:
#install.packages("lme4")
#install.packages("lmerTest")
#install.packages("sjPlot")
library( lme4 ) # zur Schaetzung von Multilevel-Modellen
## Lade nötiges Paket: Matrix
library( lmerTest ) # fuer die Signifikanztests der festen Effekte
##
## Attache Paket: 'lmerTest'
## Das folgende Objekt ist maskiert 'package:lme4':
##
## lmer
## Das folgende Objekt ist maskiert 'package:stats':
##
## step
library( sjPlot ) # berichten eines LMMs
#### ++ Beispiel ###
# mit N = 3758 Schüler*innen (Level 1) an 211 Schulen (Level 2) in den
# Niederlanden wurde ein Sprachtest durchgeführt, zusätzlich dazu liegen Daten
# zur verbalen Intelligenz und zum sozialen Status der Schule vor.
# school: Nummer der Schule (i = 1, …, 211)
# sprache: Ergebnisse der Schüler*innen im Sprachtest → Level 1-Prädiktor
# iq_verb: Ergebnisse der Schüler*innen im IQ-Test (zentriert) → Level 1-Prädiktor
# school_ses: sozialer Status der Schule (zentriert) → Level 2-Prädiktor
###############################
# Berechnung ICC ----
# Schritt 1: Random Intercept-Only Modell fitten
# Kurzbeschreibung zu Schritt 1
# In unserem Beispiel haben wir Schueler*innen in Klassen erfasst und die
# Klassenvariable heisst "school".
# Regressionskonstanten werden in R mit der "1" angeschrieben, wenn Random
# Intercept. Die lmer Funktion erlaubt es uns das Modell zu berechnen.
m1 <- sprache ~ 1 + (1 | school) # 1 steht für Intercept
fit1 <- lme4::lmer( m1, data = daten, REML = TRUE )
summary( fit1 )
## Linear mixed model fit by REML ['lmerMod']
## Formula: sprache ~ 1 + (1 | school)
## Data: daten
##
## REML criterion at convergence: 26595.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1850 -0.6417 0.0905 0.7226 2.5281
##
## Random effects:
## Groups Name Variance Std.Dev.
## school (Intercept) 18.24 4.271
## Residual 62.85 7.928
## Number of obs: 3758, groups: school, 211
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 41.0038 0.3257 125.9
# Schritt 2: ICC ablesen
tab_model(fit1)
|
|
sprache
|
|
Predictors
|
Estimates
|
CI
|
p
|
|
(Intercept)
|
41.00
|
40.37 – 41.64
|
<0.001
|
|
Random Effects
|
|
σ2
|
62.85
|
|
τ00 school
|
18.24
|
|
ICC
|
0.22
|
|
N school
|
211
|
|
Observations
|
3758
|
|
Marginal R2 / Conditional R2
|
0.000 / 0.225
|
# alternativ per 'Hand' mit ANgaben aus der Summary:
(icc <- 18.24/(18.24+62.85))
## [1] 0.2249353
# Interpretation: → bedeutet das 22,49% der Unterschiede von Y (= hier der
# Sprachtestvariable) auf Unterscheide zwischen den Schulen
# zurückgehen.
###############################
# Random Intercept-Fixed Slopes Modell mit Level 1 & Level 2 -Prädiktoren ----
# In diesem Modell wird iq_verb und school_ses zur Vorhersage der Punktzahl im
# Sprachtest benutzt. Die Regressionskonstante zwischen den Schulen darf zufällig
# variieren.
#Schritt 1: Modell definieren und fitten.
m2 <- sprache ~ 1 + iq_verb + school_ses + (1 |school)
fit2 <- lmer( m2, data = daten,
REML = TRUE,
lmerControl(optimizer = "bobyqa") ) # hilft bei Konvergenzproblemen
summary( fit2 )
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: m2
## Data: daten
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 24905.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2285 -0.6403 0.0661 0.7122 3.2056
##
## Random effects:
## Groups Name Variance Std.Dev.
## school (Intercept) 9.162 3.027
## Residual 40.455 6.360
## Number of obs: 3758, groups: school, 211
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.111e+01 2.369e-01 2.006e+02 173.561 < 2e-16 ***
## iq_verb 2.485e+00 5.466e-02 3.707e+03 45.469 < 2e-16 ***
## school_ses 1.538e-01 3.793e-02 2.120e+02 4.054 7.07e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) iq_vrb
## iq_verb -0.003
## school_ses 0.049 -0.112
tab_model(fit2)
|
|
sprache
|
|
Predictors
|
Estimates
|
CI
|
p
|
|
(Intercept)
|
41.11
|
40.64 – 41.57
|
<0.001
|
|
iq verb
|
2.49
|
2.38 – 2.59
|
<0.001
|
|
school ses
|
0.15
|
0.08 – 0.23
|
<0.001
|
|
Random Effects
|
|
σ2
|
40.45
|
|
τ00 school
|
9.16
|
|
ICC
|
0.18
|
|
N school
|
211
|
|
Observations
|
3758
|
|
Marginal R2 / Conditional R2
|
0.365 / 0.483
|
# Schritt 2: Interpretation
# Fixed Effects
# Intercept c00: Die vorhergesagte Testleistung im Sprachtest bei einem durchschnittlichen
# Ergebnis des IQ-Tests (da zentriert) beträgt über alle Schulen mit einem
# mittleren SES hinweg 41.11 Punkte.
#
# iq verb c10: Die vorhergesagte Testleistung im Sprachtest steigt über alle Schulen
# hinweg um 2.49 Punkte, wenn Schüler:innen einen um einen Punkt höheren
# verbalen IQ haben unter Kontrolle für den SES der Schule.
#
# school ses c01: Die erwartete Testleistung im Sprachtest steigt über alle Schulen
# hinweg um 0.15, wenn der SES einer Schule sich um 1 Einheit erhöht
# bei Konstanthaltung des verbalen IQs.
#
# Random Effects
#
# intercept u0j : Die Schulen unterscheiden sich mit einer Varianz von 9.162 in ihrer
# vorhergesagten Testleistung im Sprachtest bei durchschnittlicher Leistung
# im IQ-Tests und SES.
####### BEARBEITUNGSNOTIZ
###############################
## Random Intercept-Random Slopes mit Level 1 Prädiktor & Level 2 -Prädiktoren ----
# Unser Ziel ist es nun Unterschiede (= Variabilität) zwischen den Regressions-
# koeffizienten der Level-2 Einheiten (= hier Schulen) durch den SES der Schule
# (= school_ses), also dem Level 2 Prädiktor zu erklären.
#Schritt 1: Modell definieren und fitten.
m3 <- sprache ~ 1 + iq_verb + school_ses + iq_verb*school_ses + (1 + iq_verb |school)
fit3 <- lmer( m3, data = daten, REML = TRUE, lmerControl(optimizer = "bobyqa") )
summary( fit3 )
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: m3
## Data: daten
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 24880.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2680 -0.6263 0.0698 0.7081 2.7911
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## school (Intercept) 9.0572 3.0095
## iq_verb 0.1809 0.4253 -0.72
## Residual 39.7556 6.3052
## Number of obs: 3758, groups: school, 211
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 41.208027 0.236815 206.710183 174.009 < 2e-16 ***
## iq_verb 2.492022 0.062884 173.594666 39.629 < 2e-16 ***
## school_ses 0.143460 0.037706 214.471488 3.805 0.000185 ***
## iq_verb:school_ses -0.023333 0.009534 152.970729 -2.447 0.015524 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) iq_vrb schl_s
## iq_verb -0.321
## school_ses 0.051 -0.115
## iq_vrb:sch_ -0.110 0.061 -0.289
tab_model(fit3)
|
|
sprache
|
|
Predictors
|
Estimates
|
CI
|
p
|
|
(Intercept)
|
41.21
|
40.74 – 41.67
|
<0.001
|
|
iq verb
|
2.49
|
2.37 – 2.62
|
<0.001
|
|
school ses
|
0.14
|
0.07 – 0.22
|
<0.001
|
|
iq verb × school ses
|
-0.02
|
-0.04 – -0.00
|
0.014
|
|
Random Effects
|
|
σ2
|
39.76
|
|
τ00 school
|
9.06
|
|
τ11 school.iq_verb
|
0.18
|
|
ρ01 school
|
-0.72
|
|
ICC
|
0.20
|
|
N school
|
211
|
|
Observations
|
3758
|
|
Marginal R2 / Conditional R2
|
0.369 / 0.493
|
#Schritt 2: Interpretation
# Fixed Effects
# intercept c00 : Die vorhergesagte Testleistung im Sprachtest bei einem durchschnittlichen
# Ergebnis des IQ-Tests beträgt über alle Schulen mit einem mittleren SES hinweg
# 41.20 Punkte.
#
# iq verb c10: Die vorhergesagte Testleistung im Sprachtest steigt um 2.49 Punkte,
# wenn man den verbalen IQ um eine Einheit erhöht, über alle Schulen
# mit mittlerem SES hinweg.
#
# school ses c01: Die vorhergesagte Testleistung im Sprachtest steigt für Schüler:innen
# mit durchschnittlichem Verbalen IQ bei Erhöhung des SES um eine Einheit um
# 0.14 Punkte an.
#
# iq x ses c11: signifikant, d.h. der SES scheint Unterschiede in den Regressionsgewichten
# erklären zu koennen!
# Wenn der SES um eine Einheit erhöht, dann verringert sich der erwartetete
# Einfluss des verbalen IQ pro IQ Punkt um 0.02 Punkte über alle Schulen hinweg.
# Random Effects
#
# intercept u0j : Die Schulen mit durchschnittlichem SES unterscheiden sich mit einer
# Varianz von 9.06 in ihrer vorhergesagten Testleistung im Sprachtest bei
# durchschnittlicher Leistung im IQ-Test.
#
# iq verb u1j : Die Schulen unterscheiden sich mit einer Varianz von 0.18 darin, inwiefern
# eine Erhöhung im IQ-Test zu einer Steigung der vorhergesagten Testleistung
# im Sprachtest führt (über die vorhergesagten Unterschiede durch SES hinaus).
# Hilfreiche Visualisierung der Interaktion
install.packages("interactions")
## Installiere Paket nach '/home/slk/R/x86_64-pc-linux-gnu-library/4.3'
## (da 'lib' nicht spezifiziert)
library(interactions)
interact_plot(fit3, pred = iq_verb, modx = school_ses, interval = TRUE )
