Hierarchische Regression

Daten laden

Dazu laden wir mit load die Daten (bspw. auf dem Desktop von Frau “Musterfrau”) Tipp: Verwenden Sie unbedingt die automatische Vervollständigung von R-Studio, wie in den letzten Sitzung beschrieben.

load("C:/Users/Musterfrau/Desktop/StudentsInClasses.rda")

Genauso können Sie die Daten direkt von der Website laden:

load(url("https://pandar.netlify.app/post/StudentsInClasses.rda"))

Nun sollte in R-Studio oben rechts in dem Fenster unter der Rubrik “Data” unser Datensatz mit dem Namen “StudentsInClasses” erscheinen.

Übersicht über die Daten

In Multilevel-Daten gibt es, wie der Name schon andeutet, mehrere Level. In unserem Datensatz gibt es zwei: Level 1 enthält die Individualdaten (within) und Level 2 enthält die clusterspezifischen Daten (between). Wir verwenden wieder head, um zu schauen, welche Variablen wohl auf welchem Level liegen:

head(StudentsInClasses)

##   MatheL Motivation KFT KlassenG schulklasse
## 1  48.76          4  98       26           1
## 2  46.01          3  96       26           1
## 3  65.96          5 112       26           1
## 4  42.08          4  94       26           1
## 5   0.00          2  78       26           1
## 6  56.52          5 104       26           1

Aus dieser Übersicht ergibt sich folgende Aufteilung:

Ebene 1 (within-level)

• Mathematikleistung als AV: MatheL
• Motivation der Schüler/innen als Prädiktor: Motivation
• Intelligenz der Schüler/innen als Prädiktor: KFT (Kognitiver Fähigkeitstest)

Ebene 2 (between-level)

• Klassengröße als Prädiktor: KlassenG
• Klassenzugehörigkeit als Gruppierungsvariable: schulklasse

Die Daten liegen hier im “Long”-Format vor, was bedeutet, dass Schüler aus einer Klasse untereinander stehen und es eine Clustervariable gibt, die die Schüler jeweils einer Klasse zuordnet. In diesem Fall ist dies die Level 2 Variable schulklasse (hierbei ist es wichtig, dass jedes Schulkind nur einer Schulklasse zugeordnet werden kann und dass unterschiedliche Schulklassen auch unterschiedliche Bezeichnungen haben müssen - die Level 2 Variable KlassenG kann für unterschiedliche Schulklassen die selbe Ausprägung aufweisen, sie ist nämlich nur eine Variable, die das Cluster näher beschreibt). Wir schauen uns noch schnell die Mittelwerte der Variablen an:

colMeans(StudentsInClasses)

##      MatheL  Motivation         KFT    KlassenG schulklasse 
##   53.616047    4.285882  100.001176   27.090588   20.280000

Bei Multi-Level-Analysen ist darauf zu achten, dass für bessere Interpretierbarkeit unabhängige Variablen zentriert werden sollten. Wir können zwischen Grand-Mean-Centering und Group-Mean-Centering unterscheiden. Beim Grand-Mean-Centering wird am Stichprobenmittelwert der jeweiligen Variable zentriert. Somit spricht ein Wert von Null auf der zentrierten Variable für den Mittelwert über alle Erhebungen (hier: Schulkinder). Beim Group-Mean-Centering wird am Mittelwert des jeweiligen Clusters zentriert. Somit entspricht ein Mittelwert von 0 auf der zentrierten Variable für einen durchschnittlichen Wert innerhalb dieses Clusters (hier: für einen durchschnittlichen Wert innerhalb dieser Schulklasse). Beim Group-Mean-Centering können die Ergebnisse gut in Hinsicht auf die Schulklasse interpretiert werden, allerdings geht die Unterschiedlichkeit zwischen Klassen verloren, weswegen häufig zusätzlich zur zentrierten Variable auch der Mittelwert pro Cluster als L2 Variable mit in das Modell aufgenommen wird. Wir schauen uns dies später im Abschnitt Datenzentrierung genauer an.

Hypothesen

Wir wollen folgende Hypothesen untersuchen:

1. Mathematikleistung hängt von der Lernmotivation ab
2. Der Effekt der Motivation unterscheidet sich zwischen Klassen
3. Mathematikleistung hängt von der Klassengröße ab
4. Es gibt eine Wechselwirkung (Cross-Level-Interaktion) zwischen Motivation und Klassengröße

Pakete laden

Wir werden im Folgenden wieder neue R-Pakete benötigen. Diese müssen zunächst installiert werden (install.packages) und anschließend geladen werden: Das wichtigste Paket ist lme4 (linear mixed effects 4, wobei wir annehmen können, dass die 4 hier für four steht, was sich wie “for r” spricht), welches für die Multi-Level Analysen von Nöten ist. lmerTest verwenden wir, um Modellvergleiche sinnvoll auch für Varianzen durchführen zu können. lme4 und lmerTest sollten immer gemeinsam geladen werden (auch wenn in den meisten Fällen lmerTest erzwingt, dass lme4 auch geladen wird), denn lmerTest erweitert lediglich das Paket lme4 um einige wichtige Funktionen. robutmeta verwenden wir, um Daten gruppenspezifisch zu zentrieren.

library(lme4)       # für das Durchführen von Multi-Level Regressionen
library(lmerTest) # lmerTest markiert in der Ausgabe signifikante Koeffizienten und korrigiert Modellvergleiche für Varianzen
library(robumeta)   # Datensatzmanipulation

## Warning: Paket 'robumeta' wurde unter R Version 4.1.1 erstellt

library(ggplot2) # ggplot2 und dplyr werden nur für Grafiken benötigt

## Warning: Paket 'ggplot2' wurde unter R Version 4.1.2 erstellt

library(dplyr)   # für Datensatzmanipulationen

Preliminary Analyses: Nullmodell und Intraklassenkorrelation (ICC)

Wir nennen das Nullmodell-Objekt einfach mal m0 (für Modell 0).

m0 <- lmer(MatheL ~ 1 +  (1 | schulklasse), data = StudentsInClasses)
summary(m0)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: MatheL ~ 1 + (1 | schulklasse)
##    Data: StudentsInClasses
## 
## REML criterion at convergence: 6565.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5448 -0.5642 -0.0015  0.5695  3.5476 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  schulklasse (Intercept)  23.87    4.885  
##  Residual                123.12   11.096  
## Number of obs: 850, groups:  schulklasse, 40
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  54.0235     0.8658 37.0961    62.4   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Die Summary sieht der der normalen linearen Regression sehr ähnlich. Wir schauen uns diese im Detail an.

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: MatheL ~ 1 + (1 | schulklasse)
##   Data: StudentsInClasses

zeigt uns, dass es sich um ein Linear mixed model handelt (also ein Multi-Level-Modell/hierarchische Regression), welches mit der REstricted Maximum Likelihood Methode geschätzt wurde. Bei REML werden zunächst die Random Effects geschätzt und anschließend werden erst die fixed Effects bestimmt. Damit ist das Verfahren stabiler und oft auch statistisch effizienter (kleinere SE). Außerdem wird unsere Modellgleichung nochmals wiederholt.

## 
## REML criterion at convergence: 6565.3
## 
## Scaled residuals:
##     Min      1Q  Median      3Q     Max
## -4.5448 -0.5642 -0.0015  0.5695  3.5476

gibt uns Informationen über die Konvergenz (sozusagen den Log-Likelihood-Wert), sowie Auskunft über die Verteilung de Residuen (ähnlich wie bei der Regression). Diese scheinen hier recht symmetrisch verteilt zu sein.

## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  schulklasse (Intercept)  23.87    4.885
##  Residual                123.12   11.096
## Number of obs: 850, groups:  schulklasse, 40

zeigt uns die Schätzungen der Random effects, also der Varianzen (und Kovarianzen) von \(u\) (dem zufälligen Anteil des Interzept und der Slope, falls vorhanden) sowie von \(\varepsilon\) (dem Residuum). Eine Signifikanzentscheidung wird nicht mit ausgegeben. Wir müssen diese via einem Modellvergleich, dem Likelihood-Differenzen-Test (LRT) untersuchen - dazu später mehr! Wir erkennen, dass uns sowohl die Varianzen des Interzepts und des Residuums, sowie die Standardabweichung (Std.Dev.) ausgegeben werden, wobei dies einfach nur die Wurzel aus Variance ist. Die Varianz zwischen Klassen beträgt 23.868 und die (Residual-) Varianz innerhalb der Klassen 123.117. Mit Hilfe dieser Informationen können wir den Intraklassenkorrelationskoeffizienten bestimmen. Außerdem entnehmen wir diesem Output, dass insgesamt 850 Schulkinder untersucht wurden in 40 Klassen.

# Per Hand
23.87 / (23.87 + 123.12)

## [1] 0.162392

# Mit Zugriff auf das Nullmodell-Objekt
VarCorr(m0)$schulklasse[1] / (VarCorr(m0)$schulklasse[1] + summary(m0)$sigma^2)

## [1] 0.1623841

VarCorr(m0) # nur die Standardabweichungen (also Wurzel aus den Varianzen) werden angezeigt

##  Groups      Name        Std.Dev.
##  schulklasse (Intercept)  4.8855 
##  Residual                11.0958

VarCorr(m0)$schulklasse # auch Varianzen werden angezeigt

##             (Intercept)
## (Intercept)    23.86806
## attr(,"stddev")
## (Intercept) 
##    4.885495 
## attr(,"correlation")
##             (Intercept)
## (Intercept)           1

VarCorr(m0)$schulklasse[1] # das erste Element ist die Varianz des Interzepts

## [1] 23.86806

summary(m0)$sigma # Residualstandardabweichung

## [1] 11.09581

summary(m0)$sigma^2 # Residualvarianz

## [1] 123.1171

names(summary(m0)) # alle Informationen, die wir der Summary entlocken können

##  [1] "methTitle"    "objClass"     "devcomp"      "isLmer"       "useScale"    
##  [6] "logLik"       "family"       "link"         "ngrps"        "coefficients"
## [11] "sigma"        "vcov"         "varcor"       "AICtab"       "call"        
## [16] "residuals"    "fitMsgs"      "optinfo"

## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)  54.0235     0.8658 37.0961    62.4   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pseudo \(R^2\): within

Die Varianz unseres Kriteriums \(Y\) haben wir im Null-Modell in Varianz des Cluster (Interzeptvarianz) und (within) Residualvarianz aufgeteilt: \(\mathbb{V}ar[Y] = \mathbb{V}ar[u_0]+\mathbb{V}ar[\varepsilon]\). Die Variation von \(u_0\) beschreibt somit die between-Variation, die auf Unterschiede zwischen Klassen (und somit Unterschiede, die auf die Clusterung zurückzuführen sind), während die Variation von \(\varepsilon\) gerade die within-Variation ist. Unter ideellen Bedingungen sollte nun ein L1-Prädiktor nur die within-Residualvariation verringern (dies ist dann gegeben, wenn ein L1 Prädiktor sich nicht systematisch zwischen Klassen unterscheidet, es also keine Mittelwertsunterschiede über Klassen hinweg auf dieser L1 Variable gibt!): für den Prädiktor Motivation (\(X_1\)) ergibt sich die Variation von \(Y\) als \[\mathbb{V}ar[Y] = \mathbb{V}ar[u_0]+\underbrace{\gamma_{10}^2\mathbb{V}ar[X_1]+\mathbb{V}ar[\varepsilon^*]}_{\mathbb{V}ar[\varepsilon]}.\] Für mehr Informationen zu Varianzrechenregeln siehe im Appendix B nach. Wir schreiben hier \(\varepsilon^*\), da es sich um ein anderes Residuum handelt, als im leeren Modell (\(X_1\) ist ja mit von der Partie). Dieser ideellen Gleichung ist zu entnehmen, dass die Varianz von \(\varepsilon^*\) kleiner ausfällt, also die von \(\varepsilon\), wenn \(\gamma_{10}\neq 0\) und \(X_1\) nicht konstant ist, also wenn \(X_1\) zur Vorhersage vom Kriterium beiträgt (also folgt \(\mathbb{V}ar[\varepsilon^*]<\mathbb{V}ar[\varepsilon]\)). Wir machen uns diese Gegebenheit zu nutze und quantifizieren die relative Veränderung in der Varianz des Residuums und nennen diese Pseudo-\(R^2\): \[R^2_\text{pseudo:within}=\frac{\mathbb{V}ar[\varepsilon]-\mathbb{V}ar[\varepsilon^*]}{\mathbb{V}ar[\varepsilon]}=1-\frac{\mathbb{V}ar[\varepsilon^*]}{\mathbb{V}ar[\varepsilon]}.\]

In der Realität ist es so, dass Prädiktoren je nach Zentrierung auch within und between Variation von \(Y\) erklären, da der Mittelwert eines Clusters between-Information enthält und die Abweichungen von diesem Mittelwert within-Informationen. Somit ist es so, dass ein nicht-zentrierter Prädiktor sowohl within als auch between Variation am Kriterium erklärt. Hinzu kommt dann noch die Schätzungenauigkeit.

Um Pseudo-\(R^2\) in R umzusetzen, müssen wir lediglich die Residualvarianzen der beiden Modelle abgreifen und verrechnen:

VarE0 <- summary(m0)$sigma^2 # Varianz des Residuums im Nullmodell
VarE1 <- summary(m1)$sigma^2 # Varianz des Residuums im Modell mit Motivation als Prädiktor

1- VarE1/VarE0

## [1] 0.3383368

# oder kurz:
1 - summary(m1)$sigma^2/summary(m0)$sigma^2

## [1] 0.3383368

Insgesamt können also 33.83% der Variation der Matheleistung innerhalb einer Klasse erklärt werden. Die Motivation scheint substantielle Anteile zu erklären. Die Residualvarianz im Modell mit der Motivation als Prädiktor ist um 33% geringer als im Nullmodell.

Pseudo \(R^2\): between

Genauso wie wir soeben die relative Veränderung der (within) Residualvarianz quantifiziert haben, können wir auch die relative Veränderung der Interzeptvariation (between Variation, Variation, die auf Unterschiede des Clusters/zwischen Klassen zurückzuführen sind) bestimmen. Da wir die Motivation am Gesamtmittelwert zentriert haben, sind einige Unterschiede zwischen den Klassen immer noch in der Variable Motivation_c enthalten (dazu später mehr im Abschnitt zur Datenzentrierung). Diesen Varianzanteil wollen wir nun quantifizieren und vergleichen dazu die beiden Modelle m0 und m1: \[\begin{align} \texttt{m0}:&\quad \mathbb{V}ar[Y] = \mathbb{V}ar[u_0]+\mathbb{V}ar[\varepsilon]\\ \texttt{m1}:&\quad \mathbb{V}ar[Y] = \mathbb{V}ar[u_0^*]+\gamma_{10}^2\mathbb{V}ar[X_1]+\mathbb{V}ar[\varepsilon^*] \end{align}\] Falls der Prädiktor der Motivation (Grand-Mean zentriert) auch noch Unterschiede zwischen den Klassen enthält, dann sollten gelten: \(\mathbb{V}ar[u_0^*]<\mathbb{V}ar[u_0]\); also sollte die Interzeptvariation von m0 zu m1 kleiner werden. Analog zum Pseudo-\(R^2\) zuvor wird Pseudo-\(R^2\) für die between Ebene (L2) definiert: \[R^2_\text{pseudo:between}=\frac{\mathbb{V}ar[u_0]-\mathbb{V}ar[u_0^*]}{\mathbb{V}ar[u_0]}=1-\frac{\mathbb{V}ar[u_0^*]}{\mathbb{V}ar[u_0]}.\] Wie wir die Interzeptvarianz erhalten, hatten wir uns schon angesehen, als es um den ICC ging. Folglich berechnen wir \(R^2_\text{pseudo:between}\) so:

VarU0 <- VarCorr(m0)$schulklasse[1]  # Varianz des Interzepts im Nullmodell
VarU1 <- VarCorr(m1)$schulklasse[1]  # Varianz des Interzepts im Modell mit Motivation als Prädiktor

1- VarU1/VarU0

## [1] -0.1130776

# oder kurz:
1 - VarCorr(m1)$schulklasse[1]/VarCorr(m0)$schulklasse[1]

## [1] -0.1130776

Am Ergebnis erkennen wir auch, wieso es sich hierbei nur um ein Pseudo-\(R^2\) handelt, denn die Varianz des Interzepts steigt tatsächlich von m0 zu m1. Dies kann auf die simultane Schätzung der beiden Varianzen mit Hilfe von (restricted) Maximum Likelihood zurückgeführt werden, welche nicht nur von den Varianzen der Mittelwerte abhängt, sondern auch vom Stichprobenfehler (vgl. Eid, et al., 2017, S. 752). In so einem Fall wird Pseudo-\(R^2\) auf 0 gesetzt. Folglich gehen wir davon aus, dass der erklärte between Varianzanteil der Matheleistung durch die Motivation 0 ist. Diesem Problem kann durch die zusätzliche Gewichtung der Varianzen durch die Anzahl der Level 1 und Level 2 Einheiten entgegengewirkt werden, was Sie vielleicht beim Schreiben Ihrer Masterarbeit dann durchführen könnten (vgl. Eid, et al., 2017, S. 752 und folgend).

Pseudo \(R^2\): between and within

Zum Schluss bestimmen wir nun noch den Anteil der Variation des Kriteriums der insgesamt (between und within) durch den Prädiktor erklärt werden kann. Dazu schauen wir uns die gemeinsame Veränderung der Varianzen auf Level 1 (within) und Level 2 (between) an. Erklärt der Prädiktor substantielle Anteile, so sollte gelten: \(\mathbb{V}ar[u_0^*]+\mathbb{V}ar[\varepsilon^*]<\mathbb{V}ar[u_0]+\mathbb{V}ar[\varepsilon]\); also sollte die gesamte nicht erklärte Varianz vom Kriterium von m0 zu m1 kleiner werden! Die relative Veränderung quantifizieren wir mit Hilfe von Pseudo-\(R^2\) für between und within: \[\begin{align*} R^2_\text{pseudo:bw}&=\frac{(\mathbb{V}ar[u_0]+\mathbb{V}ar[\varepsilon])-(\mathbb{V}ar[u_0^*]+\mathbb{V}ar[\varepsilon^*])}{\mathbb{V}ar[u_0]+\mathbb{V}ar[\varepsilon]}\\ &=1-\frac{\mathbb{V}ar[u_0^*]+\mathbb{V}ar[\varepsilon^*]}{\mathbb{V}ar[u_0]+\mathbb{V}ar[\varepsilon]}. \end{align*}\] Folglich berechnen wir \(R^2_\text{pseudo:bw}\) so:

VarE0 <- summary(m0)$sigma^2 # Varianz des Residuums im Nullmodell
VarE1 <- summary(m1)$sigma^2 # Varianz des Residuums im Modell mit Motivation als Prädiktor

VarU0 <- VarCorr(m0)$schulklasse[1]  # Varianz des Interzepts im Nullmodell
VarU1 <- VarCorr(m1)$schulklasse[1]  # Varianz des Interzepts im Modell mit Motivation als Prädiktor

1- (VarU1 + VarE1)/(VarU0 + VarE0)

## [1] 0.2650343

# oder kurz:
1 - (VarCorr(m1)$schulklasse[1] + summary(m1)$sigma^2)/(VarCorr(m0)$schulklasse[1] + summary(m0)$sigma^2)

## [1] 0.2650343

Insgesamt lassen sich also 26.5% der Variation der Matheleistung durch die Motivation erklären. Hierbei ist zu beachten, dass \(R^2_\text{pseudo:bw}\) nicht einfach der Mittelwert von \(R^2_\text{pseudo:between}\) und \(R^2_\text{pseudo:within}\) ist, sondern tatsächlich das gewichtete Mittel der beiden, wobei die Varianzen die Gewichte darstellen, denn ist die Residualvarianz sehr groß, dann fällt diese auch insgesamt bei der Gesamtvarianz vom Kriterium stärker ins Gewicht! Lassen wir zusätzlich unterschiedliche Steigungen zu (Random Slope), so müssen wir auch diese Varianzquelle in die Berechnung mit einbeziehen.

Test des Zufallseffekts

Um zu testen, ob die Varianz des Motivationseffekts signifikant wird, wird das Modell m1 ohne random slope mit dem Modell m2 mit random slope verglichen. Dies machen wir mit Hilfe der anova Funktion, welche in R dazu genutzt wird, um Modelle zu vergleichen. Hierbei ist es essentiell, dass wir zuvor das Paket lmerTest geladen haben. Der anova-Funktion müssen wir zunächst das restriktivere Modell (mit weniger \(df\)) übergeben und als zweites Argument das weniger restriktive Modell. Außerdem wählen wir explizit die Methode, mit welcher die Modelle verglichen werden sollen via test = "LRT". Somit verwenden wir den Likelihood-Ratio-Test. Allerdings sei an dieser Stelle erwähnt, dass dieser LRT nur eine Näherung darstellt. Um mögliche Fehlschlüsse zu reduzieren, sollten große Stichproben verwendet werden. Hierbei ist die Clusteranzahl die entscheidende Größe. Es sollten so viele Clustereinheiten wie möglich verwendet werden. Das wir Varianzen gegen Null testen stellt auch generelles Problem dar. Dies liegt daran, dass Varianzen nicht negativ sein können, wir daher eine gerichtete Hypothese formulieren müssten. Dies gilt allerdings ausschließlich für Varianzen (vgl Eid, et al., 2017, Seite 749 und folgend).

Da wir hier zunächst nur Random Effects auf Signifikanz prüfen wollen, sollten wir mit REML Schätzungen rechnen. Dazu müssen wir refit = F wählen und so der anova-Funktion zu sagen, dass REML-Schätzungen verwendet werden sollen. Würden wir fixed Effects auf Signifikanz via LRTs prüfen wollen, so müssten wir das Modell “refitten” – und zwar mit der vollen Likelihood. Wir würden dann ML verwenden.

anova(m1, m2, test = "LRT", refit = F) 

## Data: StudentsInClasses
## Models:
## m1: MatheL ~ 1 + Motivation_c + (1 | schulklasse)
## m2: MatheL ~ 1 + Motivation_c + (1 + Motivation_c | schulklasse)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## m1    4 6239.1 6258.1 -3115.6   6231.1                         
## m2    6 5781.6 5810.1 -2884.8   5769.6 461.51  2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Die LRT-Teststatistik (empirischer \(\chi^2\)-Wert) liegt bei \(\chi^2(df=2)=461.51\). Der Modellvergleichsoutput enthält folgende Informationen (dieser ist annähernd identisch aufgebaut und enthält jeweils unterschiedliche Teststatistiken, je nach dem welchen Test wir anfordern - bspw. bei der Regression würde Omnibustest/F-Test verwendet werden, entsprechend stünde dort die \(F\)-Statistik).

## Data: StudentsInClasses
## Models:
## m1: MatheL ~ 1 + Motivation_c + (1 | schulklasse)
## m2: MatheL ~ 1 + Motivation_c + (1 + Motivation_c | schulklasse)

gibt uns an, welche Daten verwendet werden und welche Modelle verglichen werden. Hier sollte das komplexere Modell an 2. Stelle stehen, da sonst die Likelihood-Differenz negativ sein könnte.

##     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
##  m1    4 6239.1 6258.1 -3115.6   6231.1                         
##  m2    6 5781.6 5810.1 -2884.8   5769.6 461.51  2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

gibt uns Kennwerte zum Modellvergleich an. Ganz vorne stehen die Namen der Modelle (der Objekte), die verglichen werden. npar gibt die zu schätzenden Parameter an. Für das Modell m1 ist dies \(\gamma_{00},\gamma_{10}\) (fixed) und \(\mathbb{V}ar[u_0],\mathbb{V}ar[\varepsilon]\) (random), während für m2 \(\gamma_{00},\gamma_{10}\) (fixed) und \(\mathbb{V}ar[u_0], \mathbb{V}ar[u_1], \mathbb{C}ov[u_0,u_1],\mathbb{V}ar[\varepsilon]\) (random) geschätzt werden müssen. Es kommen also die Varianz von \(u_1\) sowie die Kovarianz zwischen \(u_0\) und \(u_1\) hinzu. Der AIC (Akaike’s Information Criterion) sowie der BIC (Bayes’ Information Criterion) sind Informationskriterien, die die Komplexität des Modells sowie die Stichprobengröße mit einbeziehen. Kleinere Werte stehen für besseren Modell-Fit. Allerdings kann keine Aussage über die Signifikanz einer möglichen Differenzen in diesen getroffen werden. logLik gibt die Loglikelihood (hier restricted Loglikelihood) der Modelle an. Größere Werte sind wünschenswert. deviance ist hier einfach \(-2*\)logLik und entspricht quasi der Abweichung der Daten vom Modell (kleine Werte sind zu bevorzugen). Chisq ist die Likelihood-Differenz \(\chi^2_{emp}=-2(LL_{m_1}-LL_{m_2})=D_{m_1}-D_{m_2}\), die sich entweder über -2 mal der Loglikelihood (\(LL\)) Differenz bestimmen lässt oder via der Differenz der Devianzen (\(D\)). Df gibt die Freiheitsgrade der Likelihood-Differenz an. Dies ist hier gerade die Differenz der Anzahl der Parameter der beiden Modelle (siehe npar). Pr(>Chisq) ist der zugehörige \(p\)-Wert.

Interpretation: Modell m1 passt signifikant schlechter zu den Daten, \(\chi^2(df=2)=461.51; p<.001\), also ist die Varianz des Effekts signifikant von Null verschieden. Wir entscheiden uns somit für ein Modell inklusive Random Slope für die Motivation! Hypothese 2 wird somit durch die Daten gestützt. Dies bedeutet, dass sich die Beziehung zwischen Motivation und Matheleistung über die Klassen hinweg unterscheidet (in manchen Klassen hängt die Matheleistung somit stärker von der Lernmotivation eines Kindes ab als in anderen Schulklassen).

anova(m0, m1) # KEIN refit

## refitting model(s) with ML (instead of REML)

## Data: StudentsInClasses
## Models:
## m0: MatheL ~ 1 + (1 | schulklasse)
## m1: MatheL ~ 1 + Motivation_c + (1 | schulklasse)
##    npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)    
## m0    3 6572.8 6587.1 -3283.4   6566.8                        
## m1    4 6240.1 6259.1 -3116.1   6232.1 334.7  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Grafische Veranschaulichung

Die Wechselwirkung kann veranschaulicht werden, indem die Regressionsgeraden nach Klassengröße unterschieden werden (hier: über die Farbe). Die in der folgenden Grafik goldene/gelbe Geraden repräsentieren (überdurchschnittlich) große Klassen, die blauen kleine (unterdurchschnittlich große/ überdurchschnittlich kleine) Klassen. Die goldenen/gelben Linien sind steiler als die blauen Linien.

Der Grafik ist dieser Effekt deutlich zu entnehmen. Final können wir sagen, dass die individuelle Matheleistung in Schulklassen mit mehr Schülerinnen und Schülern stärker von der Lernmotivation der/des Einzelnen abhängt als in kleinen Schulklassen.

Eigentlich würde eine Cross-Level-Interaktion vor allem dann in ein Modell aufgenommen werden, wenn wir die Random Slope Variation durch eine L2 Variable erklären wollen. Wenn wir ein Modell mit ebenfalls Random Intercept Random Slope untersuchen, erkennen wir allerdings, dass die Wechselwirkung und der Haupteffekt der Klassengröße nun keinen Effekt mehr haben:

m4c <- lmer(MatheL ~ 1 + KlassenG_c + Motivation_c  + KlassenG_c:Motivation_c + 
              (1 + Motivation_c | schulklasse), 
            data=StudentsInClasses)
summary(m4c)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: MatheL ~ 1 + KlassenG_c + Motivation_c + KlassenG_c:Motivation_c +  
##     (1 + Motivation_c | schulklasse)
##    Data: StudentsInClasses
## 
## REML criterion at convergence: 5767.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2062 -0.6242  0.0259  0.6379  2.8811 
## 
## Random effects:
##  Groups      Name         Variance Std.Dev. Corr
##  schulklasse (Intercept)  28.26    5.316        
##              Motivation_c 36.74    6.062    0.39
##  Residual                 39.85    6.313        
## Number of obs: 850, groups:  schulklasse, 40
## 
## Fixed effects:
##                         Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)              53.7744     0.8889 36.9895  60.495  < 2e-16 ***
## KlassenG_c               -0.2158     0.1366 37.8116  -1.580    0.122    
## Motivation_c              6.1883     1.0062 37.6948   6.150 3.65e-07 ***
## KlassenG_c:Motivation_c   0.2148     0.1542 38.1319   1.394    0.171    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) KlssG_ Mtvtn_
## KlassenG_c  0.206               
## Motivatin_c 0.367  0.080        
## KlssnG_c:M_ 0.080  0.365  0.211

anova(m4, m4c, test = "LRT", refit = F)

## Data: StudentsInClasses
## Models:
## m4: MatheL ~ 1 + KlassenG_c + Motivation_c + KlassenG_c:Motivation_c + (1 | schulklasse)
## m4c: MatheL ~ 1 + KlassenG_c + Motivation_c + KlassenG_c:Motivation_c + (1 + Motivation_c | schulklasse)
##     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## m4     6 6212.9 6241.4 -3100.5   6200.9                         
## m4c    8 5783.3 5821.2 -2883.6   5767.3 433.66  2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Der Modellvergleich zeigt, dass ein Modell nur mit Klassengröße, Motivation und Interaktion zwischen Klassengröße und Motivation signifikant schlechter zu den Daten passt als ein Modell mit zusätzlicher Random Slope. Final würden wir uns also für ein Modell mit Random Intercept und Random Slope entscheiden, allerdings die Klassengröße als Prädiktor (und Cross-Level Interaktion) herauslassen.

Grand-Mean-Centering

Beim Grand-Mean-Centering, also der Zentrierung am Stichprobenmittelwert, ziehen wir einfach den globalen Mittelwert der Variable über alle Erhebungen und Cluster ab (das haben wir oben ja schon gemacht!): \[X^*_{ij}=X_{ij}-\bar{X}_{\cdot\cdot}\] Ein Wert von 0 auf \(X_{ij}^*\) bedeutet, dass diese Messung gerade genau dem Mittelwert der gesamten Stichprobe (über alle Personen \(i\) und alle Clustereinheiten \(j\)) entspricht.

Group-Mean-Centering

Bei der gruppenspezifischen (auch clusterspezifischen) Zentrierung (Group-Mean-Centering), ziehen wir von jeder Messung \(X_{ij}\) den Mittelwert des \(j\) Clusters auf dieser Variable ab (also bspw. den Mittelwert der Matheleistung einer Klasse ziehen wir von den individuellen Matheleistungen aller Schüler dieser Klasse ab):

\[X_{ij}^{**}=X_{ij}-\bar{X}_{\cdot j}\] Ein Wert von 0 auf \(X_{ij}^{**}\) entspricht nun dem Durchschnitt in Cluster/Gruppe \(j\). Die Variable \(X_{ij}^{**}\) sollte somit within Variation erklären (im Beispiel also die Variation innerhalb einer Klasse), was dafür spricht, dass es sinnvoll ist, diese Varianzaufklärung mit \(R^2_\text{pseudo:within}\) zu quantifizieren. Außerdem ist es sinnvoll anschließend die Variable \(\bar{X}_{\cdot j}\) mit in die Analysen aufzunehmen, da sie between-Variation enthält (also Variation, die Unterschiede zwischen Clustern darstellt). Die Varianzaufklärung von \(\bar{X}_{\cdot j}\) kann somit mit \(R^2_\text{pseudo:between}\) quantifiziert werden. Die gemeinsame Varianzaufklärung (quasi den Gesamtanteil between und within) können wir mit \(R^2_\text{pseudo:bw}\) quantifizieren.

Um diese Variablen anschließend auch noch sinnvoll hinsichtlich (z.B. Cross-Level-) Interaktionen interpretieren zu können, sollte die Gruppenmittelwertsvariable idealerweise auch noch am zentralen Mittelwert zentriert sein (der zentrale Mittelwert ist gerade auch der Mittelwert dieser Gruppenmittelwertsvariable): \[\bar{X}_{\cdot j}^*=\bar{X}_{\cdot j} - \bar{X}_{\cdot \cdot}\] Somit entspricht ein Wert von 0 auf dieser Variable gerade einer Erhebung, die dem durchschnittlichen Mittelwert über alle Cluster entspricht. Demnach können wir auch sinnvoll beide zentrierten Variablen auf einmal interpretieren. Wenn eine Person auf der Variable \(X_{ij}^{**}\) und auf der Variable \(\bar{X}_{\cdot j}^*\) einen Wert von Null aufweist, so bedeutet dies, dass die Ausprägung auf dieser Variable durchschnittlich im Bezug auf das Cluster ist und auch auf die Gesamtstichprobe (wir landen also wieder beim gesamten Durchschnitt über die Stichprobe, können aber Effekte nun im Bezug auf Unterschiede zwischen den einzelnen Gruppe und innerhalb der Gruppen auseinanderdröseln). Folgende Aussagen sind möglich: (Kinder aus) Schulklassen mit durchschnittlich höherer Motivation zeigen insgesamt im Durchschnitt eine höhere Matheleistung (quasi: wenn der Durchschnitt der Motivation höher liegt in einer Klasse, so liegt auch der Mittelwert der Matheleistung in dieser Klasse im Mittel höher) und Kinder, die im Vergleich zu ihrer Schulklasse motivierter sind, zeigen ebenfalls eine im Schnitt eine höhere Matheleistung als andere Kinder ihrer Schulklasse. Auch Interaktionen können weitere Einblicke in die Daten bringen. Natürlich könnten diese Effekte auch gegenläufig sein.

Umsetzung in R: Group-Mean-Centering

In R sieht das Ganze so aus:

# group-mean-centering:
StudentsInClasses$Motivation_groupc <- group.center(var = StudentsInClasses$Motivation, 
                                                    grp = StudentsInClasses$schulklasse)

# bestimmen der gruppenspezifischen Mittelwerte (durch Umstellen der Gleichung)
StudentsInClasses$Mot_groupmeans <- group.mean(var = StudentsInClasses$Motivation,
                                               grp = StudentsInClasses$schulklasse)

# Zentrieren der Gruppenmittelwertsvariable
StudentsInClasses$Mot_groupmeans_c <- StudentsInClasses$Mot_groupmeans -
  mean(StudentsInClasses$Motivation)

head(StudentsInClasses)

##   MatheL Motivation KFT KlassenG schulklasse Motivation_c KlassenG_c
## 1  48.76          4  98       26           1   -0.2858824  -1.090588
## 2  46.01          3  96       26           1   -1.2858824  -1.090588
## 3  65.96          5 112       26           1    0.7141176  -1.090588
## 4  42.08          4  94       26           1   -0.2858824  -1.090588
## 5   0.00          2  78       26           1   -2.2858824  -1.090588
## 6  56.52          5 104       26           1    0.7141176  -1.090588
##   Motivation_groupc Mot_groupmeans Mot_groupmeans_c
## 1             -0.36           4.36       0.07411765
## 2             -1.36           4.36       0.07411765
## 3              0.64           4.36       0.07411765
## 4             -0.36           4.36       0.07411765
## 5             -2.36           4.36       0.07411765
## 6              0.64           4.36       0.07411765

# (Spalten-)Mittelwerte (gerundet auf 10 Nachkommastellen)
round(colMeans(StudentsInClasses), 10)

##            MatheL        Motivation               KFT          KlassenG 
##         53.616047          4.285882        100.001176         27.090588 
##       schulklasse      Motivation_c        KlassenG_c Motivation_groupc 
##         20.280000          0.000000          0.000000          0.000000 
##    Mot_groupmeans  Mot_groupmeans_c 
##          4.285882          0.000000

Die Funktion group.center aus dem robumeta-Paket nimmt uns die Arbeit ab, die Daten händisch an den Gruppenmittelwerten zu zentrieren (\(X_{ij}^{**}\)). Sie nimmt 2 Argumente entgegen: var die Variable, die zentriert werden soll (hier die Motivation), und grp die Gruppierungsvariable (hier die Schulklasse). Die Funktion group.mean funktioniert analog zu group.center und gibt uns gruppenspezifische Mittelwerte aus. Zum Schluss wird noch die Gruppierungsvariable zentriert am Gesamtmittelwert: \(\bar{X}_{\cdot j}^*=\bar{X}_{\cdot j} - \bar{X}_{\cdot \cdot}\). Den Spaltenmittelwerten entnehmen wir, dass die gruppenzentrierte Variable einen Mittelwert von 0 hat und die zentrierte Gruppenmittelwertsvariable ebenfalls. Außerdem sehen wir hier nochmals, dass der Mittelwert der Motivation gleich dem Mittelwert der Gruppenmittelwertsmotivation ist (Motivation vs. Mot_groupmeans).

Diese Art der Zentrierung spielt insbesondere in Multi-Level-Analyse mit Messwiederholung von Personen (die Person ist die Clusterung und die L1 Variable sind die Messungen einer Person) eine wichtige Rolle (siehe auch Eid, et al., 2017, Kapitel 20.4).

In Appendix C wird eine Funktion präsentiert, mit welcher wir uns leicht Effekte in Multi-Level-Modellen ansehen können und so bspw. erkennen, wie unterschiedliche within und between Effekte ausfallen können - z.B. dass sie komplett gegenläufig sein können, was gleichzeitig einen wichtigen Vorteil dieser Analysemethode aufzeigt! Es wird auch das Big-Fish Little-Pond Prinzip daran erläutert.

Appendix A

# Histogramm der Klassenspezifischen Koeffizienten
model <- m2
beta <- coef(model)$schulklasse   # Übergeben der Koeffizienten pro Klasse (Interzept und Steigungskoeffizent)
hist(beta[,2], breaks = 10, freq = F, col = "skyblue", border = "blue",
     main = "Verteilung der Steigungskoeffizienten der Motivation",
     xlab = expression("Steigungskoeffizient Motivation:"~beta[1~"j"]~"="~gamma[10]+u[1~"j"])) # Histogramm
gamma10 <- mean(beta[,2]) # Mittlere Steigung
VarU1 <- var(beta[,2]) # Varianz der Steigung
abline(v=gamma10, col = "red", lwd = 5)   # Mittlere Steigung in der Grafik
text(x = gamma10+2, y = 0.05, labels = expression(gamma[10]), cex = 3, col = "red") # Bennenung gamma01
lines(x=seq(-10,20,0.01), dnorm(x = seq(-10,20,0.01), mean = gamma10,
                                sd = sqrt(VarU1)), col = "darkblue", lwd = 3) # Normalverteilung als Vergleich
arrows(y0 = dnorm(x = gamma10+sqrt(VarU1), mean = gamma10,
                  sd = sqrt(VarU1)), y1 = dnorm(x = gamma10+sqrt(VarU1), mean = gamma10,sd = sqrt(VarU1)),
       x0 = gamma10-sqrt(VarU1), x1 = gamma10+sqrt(VarU1),
       code = 3, col = "blue", lwd = 3, angle = 90) # +/- 1 SD in der Grafik
text(x = 17, y= 0.04, labels = "+/- 1SD", col = "blue", cex = 2) # Text in Grafik einfügen

 #### Das Modell zeichnen!
model <- m2
beta <- coef(model)$schulklasse   # Übergeben der Koeffizienten pro Klasse (Interzept und Steigungskoeffizent)
summary_model <- summary(model)
library(ggplot2)
ggplot(data = StudentsInClasses, aes(x=Motivation_c, y = MatheL))+geom_point(col = "grey50")+
  theme_minimal()  + geom_abline(intercept =beta$`(Intercept)`,          # Interzepts
                                 slope = beta$Motivation_c, lwd = .2) +   # und Slopes für Grafiken verwenden
  geom_abline(intercept = summary_model$coefficients[1,1],          # Interzepts und Slope von Durchschnitt nehmen
              slope = summary_model$coefficients[2,1], lwd = 2, col = "blue")

library(ggplot2) # Lade ggplot2
library(dplyr) # Lade dplyr, um Datensätze zu manipulieren (via "mutate" und "%>%")

StudentsInClasses %>%             # Datensatz wird manipuliert
  mutate(Leistung = fitted(m4)) %>% # Vorhergesagte Werte in den Datensatz via Funktion fitted(), nenne diese Variable Leistung
  ggplot(aes(x = Motivation_c, Leistung, group = schulklasse, color = KlassenG_c)) +
  scale_color_gradient(low="blue", high = "gold3") +     # von blau bis gold skalieren
  theme_classic() +                                      # klassisches Thema wählen, sodass der Hintergrund nicht grau ist
  geom_point(aes(y = MatheL), alpha = 0.1, color="grey")+ # beobachtete Werte als Punkte
  geom_line(size=0.5)  # Vorhergesagte Werte als Linien

Appendix B

Appendix C

plot_within_between_effects <- function(nb = 50, nw = 50, between_effect = 1, within_effect = 1)
{
       # Wiederholungsfunktion
        reps <- function(X, r)
        {
                out <- c()
                for(x in X)
                {
                        out <- c(out, rep(x, r))
                }
                out
        }

        # Daten generieren (Normalverteilt)
        Xb <- rnorm(nb, mean = 0)
        Yb <- between_effect*Xb + rnorm(nb)

        # between Effekte
        between <- cbind(Xb, Yb)
        between <- apply(between, 2, FUN = function(x) reps(X = x, r = nw))

        # Within Effekte
        Xw <- .1*rnorm(dim(between)[1])
        Yw <- within_effect*Xw + .1*rnorm(dim(between)[1])
        within <- cbind(Xw, Yw)

        # Gesamteffekte
        total <- between + within
        Cluster <- rep(1:nb, nw); Cluster <- sort(Cluster)

        total <- cbind(total, Cluster)
        total <- data.frame(total)
        names(total) <- c("X", "Y", "Cluster")

        total$Cluster <- as.factor(total$Cluster)

        # Grafik erstellen
        library(ggplot2)
        ggplot(data = total, mapping = aes(x = X, y = Y, col = Cluster))+geom_point()+
          theme_minimal()
}

plot_within_between_effects(nb = 50, nw = 50, between_effect = -1, within_effect = 1)