1 Modelle mit Level-1-Prädiktoren

In diesem Kapitel behandeln wir Hierarchische Lineare Modelle mit Level-1-Prädiktoren. Aufbau und Systematik der Darstellung orientieren sich am Lehrbuch “Statistik und Forschungsmethoden” von Eid, Gollwitzer, und Schmitt (2017).

Es handelt sich um Modelle mit einer 2-Ebenen-Struktur (Individuen genestet in Gruppen), und zunächst werden nur Modelle ohne jede Prädiktorvariable (Intercept-Only-Modell) sowie Modelle mit ausschliesslich Prädiktorvariablen auf Level 1 (Random-Intercept-Modell und Random-Coefficients-Modell) behandelt.

Laden der in diesem Kapitel benötigten Packages:

pacman::p_load(lme4, nlme, tidyverse, lmerTest, gridExtra, ggplot2)

1.1 HLM Packages

Zu Beginn ein kurzer Überblick über die beiden Packages, die wir in diesem Kurs für Hierarchische Lineare Modelle benutzen werden.

nlme

Dies ist das ursprüngliche HLM- oder Mixed-Model-Package. Es hat mehr Optionen, insbesondere im Hinblick auf die Struktur der Level-1-Fehler bei Daten mit wiederholten Messungen. Wir verwenden nlme im Kapitel zu HLM-Modellen mit Messwiederholung. Die Syntax-Struktur wird dort behandelt.

lme4

Dieses modernere HLM-Package ist sehr gut geeignet für “normale” Hierarchische Lineare Modelle mit genesteten Gruppen. Wir verwenden es im Folgenden für die Analysen der Modelle mit Level-1- und Level-2-Prädiktoren.

Die “lme4”-Syntax baut auf der model syntax auf, die wir schon von lm() kennen. Also:

lm(formula = AbhängigeVariable ~ UnabhängigeVariable, data = dataframe)

Bei der Funktion lmer() (aus dem Package lme4) kommt die Spezifizierung der Gruppenvariable hinzu, und welche zufälligen Effekte geschätzt werden sollen. Die Syntax in der Klammer unterscheidet die lmer()-Syntax von der simplen Syntax der lm()-Funktion. In der Klammer wird links vom | spezifiziert, welche zufälligen Effekte geschätzt werden sollen. Rechts des | steht die Variable, die die Gruppenstruktur (Nestungsstruktur) der Daten definiert.

Eine 1 im linken Teil der Klammer bedeutet, dass random intercepts geschätzt werden sollen. Wenn die Unabhängige Variable (also die Prädiktorvariable) auch links vor dem | steht, bedeutet dies, dass zusätzlich random slopes geschätzt werden.

Was bedeutet also folgende Syntax?

lmer(data = dataframe, AbhängigeVariable ~ UnabhängigeVariable + (1 | Gruppenvariable))

Antwort: die Abhängige Variable wird von der Unabhängigen Variable vorhergesagt. Gleichzeitig wird der Intercept separat für jede Ausprägung der Gruppenvariable geschätzt bzw. die Varianz der Level-2-Residuen des Intercepts ist ein Parameter des Modells.

Diese Syntax wird später noch erweitert.

1.2 Daten laden und aufbereiten

Wir laden die Daten und speichern sie im Objekt df (für “data frame”). Wir können das Datenfile direkt aus dem Internet (GitHub) herunterladen:

df <- read_csv(
  url("https://raw.githubusercontent.com/methodenlehre/data/master/statistik_IV/salary-data.csv")
)

Beschreibung

Bei diesen Daten handelt es sich um Lohn-Daten von 20 Firmen mit je 30 Angestellten (fiktive, simulierte Daten). Zu jedem Angestellten haben wir Informationen zum Gehalt (salary, in 1000 CHF), bei welcher Firma die Person angestellt ist (company), und wie lange die Person schon dort arbeitet (experience). Ausserdem haben wir noch Informationen darüber, in welchem Sektor (public oder privat) die Firma tätig ist (sector). Mit dieser Variable beschäftigen wir uns aber erst in Übung 2.

Bevor wir uns die Daten anschauen, müssen company und sector noch als Faktoren definiert werden.

df <- df |>
  mutate(
    company = as.factor(company),
    sector = as.factor(sector)
  )

Jetzt können wir uns den Daten-Frame anschauen:

# Die Funktion "tail()" gibt uns die letzten 6 rows aus.
tail(df)

# A tibble: 6 × 4
  company    experience salary sector 
  <fct>           <dbl>  <dbl> <fct>  
1 Company 20       3.58   6.84 Private
2 Company 20       3.18   7.60 Private
3 Company 20       3.39   5.71 Private
4 Company 20       7.12  10.1  Private
5 Company 20       2.98   6.94 Private
6 Company 20       6.45   9.33 Private

# Die "summary()" Funktion gibt uns deskriptive Statistiken
# zum Datensatz aus.
summary(df)

       company      experience         salary           sector   
 Company 01: 30   Min.   : 0.000   Min.   : 4.587   Private:300  
 Company 02: 30   1st Qu.: 4.027   1st Qu.: 7.602   Public :300  
 Company 03: 30   Median : 5.170   Median : 8.564                
 Company 04: 30   Mean   : 5.190   Mean   : 8.738                
 Company 05: 30   3rd Qu.: 6.402   3rd Qu.: 9.840                
 Company 06: 30   Max.   :10.000   Max.   :15.418                
 (Other)   :420

1.3 Intercept-Only-Modell (Modell 1)

In diesem Modell ist keine Prädiktorvariable enthalten, es handelt sich also um eine Art Nullmodell, das nur die Variation des Achsenabschnitts zwischen den Gruppen modelliert. Und in einem Modell ohne Prädiktor beziehen sich die Gruppenabschnitte auf die Gruppenmittelwerte (in diesem Fall der 20 Firmen).

Level-1-Modell: \(\boldsymbol {y_{mi}=\beta_{0i}+\varepsilon_{mi}}\)

Level-2-Modell: \(\boldsymbol {\beta_{0i}=\gamma_{00} + \upsilon_{0i}}\)

Gesamtmodell: \(\boldsymbol {y_{mi} = \gamma_{00} + \upsilon_{0i}+\varepsilon_{mi}}\)

1.3.1 Modell berechnen

# Fitten des Modells und Abspeichern in der Variable intercept.only.model
intercept.only.model <- lmer(salary ~ 1 + (1 | company), data = df, REML = TRUE)

# Abspeichern von Prädiktionen, die dieses Modell macht.
# In diesem Falle ist dies jeweils der durchschnittliche Lohn jeder Firma
# Die Funktion predict() gibt pro Zeile (in diesem Fall pro Angestelltem)
# eine Prädiktion aus, je nachdem, in welcher Firma die Person angestellt ist.
# Wir speichern diese vorhergesagten Werte, um das Modell später mit "ggplot2"
# zu veranschaulichen.
df$intercept.only.preds <- predict(intercept.only.model)

# Model output anschauen
summary(intercept.only.model)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: salary ~ 1 + (1 | company)
   Data: df

REML criterion at convergence: 2158

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.9816 -0.6506 -0.0494  0.5779  4.2131 

Random effects:
 Groups   Name        Variance Std.Dev.
 company  (Intercept) 0.8512   0.9226  
 Residual             1.9547   1.3981  
Number of obs: 600, groups:  company, 20

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)   8.7376     0.2141 19.0000   40.82   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mit der Funktion ranef() kann man sich die einzelnen Random-Effects (Level-2-Residuen des Intercepts, also die \(\upsilon_{0i}\)) ausgeben lassen:

ranef(intercept.only.model)

$company
           (Intercept)
Company 01  0.78965416
Company 02  1.10501807
Company 03  1.92302618
Company 04 -1.13650080
Company 05  0.95354595
Company 06 -0.95893185
Company 07 -0.01014627
Company 08 -0.25484048
Company 09 -0.65131802
Company 10  0.76848590
Company 11 -0.50626620
Company 12  0.94040709
Company 13 -0.74284383
Company 14 -0.97545482
Company 15 -1.16136931
Company 16 -0.09768008
Company 17  0.66196052
Company 18 -0.16811195
Company 19  0.35123926
Company 20 -0.82987351

with conditional variances for "company"

1.3.2 Intraklassenkorrelation

Die Intraklassenkorrelation quantifiziert das Ausmass der „Nicht-Unabhängigkeit” zwischen den Messwerten aufgrund systematischer Level-2-Unterschiede im gemessenen Merkmal. Je grösser der Anteil der Level-2-Varianz (Varianz der Gruppenmittelwerte) an der Gesamtvarianz (Summe der Level-2-Varianz und der Level-1-Varianz), desto grösser sind die Ähnlichkeiten innerhalb der Level-2-Einheiten im Vergleich zu zwischen den Level-2-Einheiten.

Die Intraklassenkorrelation ist definiert als: \(\boldsymbol {\rho = \frac{\sigma^2_{Level-2}}{\sigma^2_{Level-2}+ \sigma^2_{Level-1}}}\)

Die Intraklassenkorrelation erhält man bei Schätzung eines Nullmodells (Intercept-Only-Modell, s.o.), bei dem sowohl die (zufällige) Varianz des Intercepts (in einem Modell ohne Prädiktoren also die Varianz der Mittelwerte der Level-2-Einheiten) als auch die Level-1-Residualvarianz ausgegeben wird.

Wenn die Level-2-Varianz sich nicht überzufällig von 0 unterscheidet, spielen die Ähnlichkeiten/Abhängigkeiten innerhalb der Level-2-Einheiten keine Rolle und ein Mehrebenenmodell ist nicht unbedingt notwendig.

Intraklassenkorrelation: \(\hat{\rho} = \frac{\hat{\sigma}^2_{\upsilon_0}}{\hat{\sigma}^2_{\upsilon_0}+ \hat{\sigma}^2_\varepsilon} = \frac {0.8512}{0.8512+1.9547} = 0.303\)

-> ca. 30 % der Gesamtvarianz sind auf Level-2-(Firmen-)Unterschiede zurückzuführen.

1.3.3 Signifikanztest

Signifikanztest für \(\sigma^2_{\upsilon_0}\): Modellvergleich mit absolutem Nullmodell

Wir kennen nun also die Level-2-Varianz, aber wir haben noch keinen Signifikanztest für diesen Parameter. Diesen erhalten wir über einen Modellvergleich (Likelihood-Ratio-Test) des Intercept-Only-Modells mit einem Modell, dass keinen Random-Intercept enthält, d.h. mit einem “normalen” linearen Modell ohne Prädiktor (mit dem Gesamt-Mittelwert \(\gamma_{00}\) als einzigem Modellparameter neben der Level-1-Residualvarianz \(\sigma^2_\varepsilon\)). Dieses Modell kann man auch als “absolutes Nullmodell” bezeichnen. Wir müssen ein solches Modell nicht extra spezifizieren, sondern können die Funktion ranova() auf das Outputobjekt intercept.only.model anwenden. ranova() (aus dem lme4-Hilfspackage lmerTest) führt automatisch Modellvergleiche (nur) für Random-Effects durch, indem die vorhandenen Random-Effects Schritt für Schritt entfernt werden und das Ausgangsmodell dann mit dem so reduzierten Modell verglichen wird. In diesem Fall (Intercept-Only-Modell als Ausgangsmodell) kann nur ein Random-Effect entfernt werden, nämlich der des Intercepts.

# Vergleich des Intercept-Only-Modells mit dem "absoluten Nullmodell"
ranova(intercept.only.model)

ANOVA-like table for random-effects: Single term deletions

Model:
salary ~ (1 | company)
              npar  logLik    AIC    LRT Df Pr(>Chisq)    
<none>           3 -1079.0 2164.0                         
(1 | company)    2 -1157.7 2319.4 157.44  1  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Der Vergleich ist signifikant \((p<0.001)\). Wir können den p-Wert sogar noch halbieren, da es sich um einen einseitigen Test handelt. Der zu testende Parameter \(\sigma^2_{\upsilon_0}\) weist nämlich einen bounded parameter space auf (kann nicht kleiner als 0 werden). Damit haben wir eine signifikante Level-2-Varianz des Intercepts und die Verwendung eines hierarchischen linearen Modells für die folgenden Analysen ist damit angezeigt.

Weiter unten werden wir zeigen, wie Modellvergleiche mittels Likelihood-Ratio-Test auch von Hand durchgeführt werden können.

Visuell kann man sich den Vergleich zwischen dem absoluten Nullmodell und dem Intercept-Only-Modell folgendermassen vorstellen:

Code

xlabel <- "Experience in years"
ylabel <- "Salary per month"


# Visual comparison between the two models
null.model <- lm(salary ~ 1, data = df)
df$null.model.preds <- predict(null.model)

ggplot(data = df, aes(x = experience, y = null.model.preds)) + # , group=groupVar, colour=groupVar))+
  geom_smooth(method = "lm", fullrange = TRUE, se = F, size = .3) +
  geom_jitter(aes(group = company, colour = company), show.legend = FALSE, alpha = .2) +
  labs(x = xlabel, y = ylabel) +
  ggtitle("Absolutes Nullmodell") +
  scale_colour_discrete() +
  ggthemes::theme_tufte()

Code

df$intercept.only.preds <- predict(intercept.only.model)

ggplot(data = df, aes(
  x = experience,
  y = intercept.only.preds,
  colour = company
)) +
  geom_smooth(method = "lm", fullrange = TRUE, se = F, size = .3) +
  geom_jitter(aes(y = salary), alpha = .2) +
  labs(x = xlabel, y = ylabel) +
  ggtitle("Intercept-Only-Modell") +
  scale_colour_discrete() +
  geom_abline(intercept = fixef(intercept.only.model), slope = 0) +
  ggthemes::theme_tufte()

Vertiefung: Shrinkage

Wir haben oben bereits angemerkt, dass in einem Intercept-Only-Modell die Gruppen-Intercepts den salary Firmenmittelwerten “entsprechen”. Die geschätzten Gruppenmittelwerte \((\hat\gamma_{00} + \upsilon_{0i})\) erhält man mit coef(). Vergleichen wir diese mit den empirischen Durchschnittsgehältern pro Firma:

est_means <- coef(intercept.only.model)$company[, 1]
emp_means <- df |>
  group_by(company) |>
  summarise(emp_means = mean(salary)) |>
  ungroup() |>
  dplyr::select(emp_means)

salary <- tibble(est_means, emp_means) |>
  mutate(diff_means = est_means - emp_means)

salary

# A tibble: 20 × 3
   est_means emp_means diff_means
       <dbl>     <dbl>      <dbl>
 1      9.53      9.59  -0.0604  
 2      9.84      9.93  -0.0846  
 3     10.7      10.8   -0.147   
 4      7.60      7.51   0.0870  
 5      9.69      9.76  -0.0730  
 6      7.78      7.71   0.0734  
 7      8.73      8.73   0.000777
 8      8.48      8.46   0.0195  
 9      8.09      8.04   0.0499  
10      9.51      9.56  -0.0588  
11      8.23      8.19   0.0388  
12      9.68      9.75  -0.0720  
13      7.99      7.94   0.0569  
14      7.76      7.69   0.0747  
15      7.58      7.49   0.0889  
16      8.64      8.63   0.00748 
17      9.40      9.45  -0.0507  
18      8.57      8.56   0.0129  
19      9.09      9.12  -0.0269  
20      7.91      7.84   0.0635

Man kann erkennen, dass die geschätzten Intercepts pro Firma nicht mit den tatsächlichen empirischen Gruppenmitteln übereinstimmen. In einem Mehrebenenmodell werden diese Schätzungen mit “Shrinkage” (Schrumpfung) berechnet, d.h. sie werden in Richtung des Gesamtmittelwertes “gezogen”. Dies liegt daran, dass die Level-2-Residuen in diesem Modell als normalverteilt angenommen werden und die Schätzunsicherheit in dem Sinne berücksichtigt wird, dass die Firmenmittelwerte der AV salary durch Gewichtung (“Pooling”) mit dem Gesamtmittelwert korrigiert werden, um zu den Estimates \(\beta_{0i}=\hat\gamma_{00} + \upsilon_{0i}\) zu gelangen.

Dies ist auch der Grund, warum Mehrebenenmodelle manchmal als “Partial-Pooling-Modelle” bezeichnet werden (z.B. bei Gelman und Hill 2007), im Gegensatz zu “No-Pooling-Modellen” (lineare Modelle mit Gruppen als fixed effects ohne jegliche Annahmen bezüglich der Varianz auf Ebene 2).

Schauen wir uns zum Vergleich ein solches “No-Pooling-Modell” an:

no.pooling <- lm(salary ~ company, data = df)
summary(no.pooling)


Call:
lm(formula = salary ~ company, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.2291 -0.8918 -0.0401  0.7795  5.8300 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        9.58770    0.25526  37.560  < 2e-16 ***
companyCompany 02  0.33950    0.36099   0.940 0.347369    
companyCompany 03  1.22013    0.36099   3.380 0.000774 ***
companyCompany 04 -2.07359    0.36099  -5.744 1.49e-08 ***
companyCompany 05  0.17644    0.36099   0.489 0.625201    
companyCompany 06 -1.88243    0.36099  -5.215 2.57e-07 ***
companyCompany 07 -0.86102    0.36099  -2.385 0.017393 *  
companyCompany 08 -1.12444    0.36099  -3.115 0.001931 ** 
companyCompany 09 -1.55127    0.36099  -4.297 2.03e-05 ***
companyCompany 10 -0.02279    0.36099  -0.063 0.949687    
companyCompany 11 -1.39512    0.36099  -3.865 0.000124 ***
companyCompany 12  0.16229    0.36099   0.450 0.653188    
companyCompany 13 -1.64980    0.36099  -4.570 5.96e-06 ***
companyCompany 14 -1.90022    0.36099  -5.264 1.99e-07 ***
companyCompany 15 -2.10036    0.36099  -5.818 9.84e-09 ***
companyCompany 16 -0.95525    0.36099  -2.646 0.008362 ** 
companyCompany 17 -0.13747    0.36099  -0.381 0.703488    
companyCompany 18 -1.03108    0.36099  -2.856 0.004441 ** 
companyCompany 19 -0.47197    0.36099  -1.307 0.191585    
companyCompany 20 -1.74349    0.36099  -4.830 1.75e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.398 on 580 degrees of freedom
Multiple R-squared:  0.3154,    Adjusted R-squared:  0.293 
F-statistic: 14.06 on 19 and 580 DF,  p-value: < 2.2e-16

Das Modell ohne Pooling entspricht einem ANOVA-Modell, und die Effekte vergleichen die Gehaltsmittelwerte der Firmen 02-20 mit dem von Firma 01 (d.h. Effekte von dummy-codierten Variablen mit Firma 01 als Referenzkategorie). Das \(R^2\) dieses Modells mit festen Effekten (fixed effects ANOVA) ist etwas höher als die Intraklassen-Korrelation des Intercept-Only-Modells (0.3154 vs. 0.3034), was eben die im Mehrebenenmodell vorhandene Shrinkage widerspiegelt.

Um ein Modell mit Schätzungen für alle Unternehmensmittelwerte (statt der Dummy-Effekte) zu erhalten, können wir den Intercept des lm()-Modells weglassen:

no.pooling2 <- lm(salary ~ company - 1, data = df)
summary(no.pooling2)


Call:
lm(formula = salary ~ company - 1, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.2291 -0.8918 -0.0401  0.7795  5.8300 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
companyCompany 01   9.5877     0.2553   37.56   <2e-16 ***
companyCompany 02   9.9272     0.2553   38.89   <2e-16 ***
companyCompany 03  10.8078     0.2553   42.34   <2e-16 ***
companyCompany 04   7.5141     0.2553   29.44   <2e-16 ***
companyCompany 05   9.7641     0.2553   38.25   <2e-16 ***
companyCompany 06   7.7053     0.2553   30.19   <2e-16 ***
companyCompany 07   8.7267     0.2553   34.19   <2e-16 ***
companyCompany 08   8.4633     0.2553   33.16   <2e-16 ***
companyCompany 09   8.0364     0.2553   31.48   <2e-16 ***
companyCompany 10   9.5649     0.2553   37.47   <2e-16 ***
companyCompany 11   8.1926     0.2553   32.09   <2e-16 ***
companyCompany 12   9.7500     0.2553   38.20   <2e-16 ***
companyCompany 13   7.9379     0.2553   31.10   <2e-16 ***
companyCompany 14   7.6875     0.2553   30.12   <2e-16 ***
companyCompany 15   7.4873     0.2553   29.33   <2e-16 ***
companyCompany 16   8.6324     0.2553   33.82   <2e-16 ***
companyCompany 17   9.4502     0.2553   37.02   <2e-16 ***
companyCompany 18   8.5566     0.2553   33.52   <2e-16 ***
companyCompany 19   9.1157     0.2553   35.71   <2e-16 ***
companyCompany 20   7.8442     0.2553   30.73   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.398 on 580 degrees of freedom
Multiple R-squared:  0.9761,    Adjusted R-squared:  0.9753 
F-statistic:  1185 on 20 and 580 DF,  p-value: < 2.2e-16

Nebenbemerkung: Da das Modell keinen Intercept enthält, beinhaltet die F-Statistik des Modells nun auch den Test der Gruppenmittelwerte gegen 0, so dass wir sie nicht für den Vergleich der Gehaltsmittelwerte zwischen den Firmen verwenden können.

Überprüfen wir nun, ob die Schätzungen dieses Modells wirklich mit den Gruppenmittelwerten übereinstimmen:

# "no-pooling" Mittelwerte aus no.pooling2 extrahieren
np_means <- as.numeric(coef(no.pooling2))
# nicht benötigte Attribute entfernen
emp_means <- as.numeric(unlist(emp_means))
# Äquivalenz überprüfen
all.equal(np_means, emp_means)

[1] TRUE

Zum Schluss plotten wir noch die mit dem “partial pooling” Intercept-Only-Modell geschätzten Mittelwerte und die “no pooling” Mittelwerte des linearen Modells, also die empirischen Firmenmittelwerte:

Code

salary |>
  pivot_longer(-diff_means, names_to = "means", values_to = "salary") |>
  mutate(means = as.factor(means)) |>
  ggplot(aes(x = means, y = salary)) +
  geom_boxplot(aes(color = means)) +
  geom_point(alpha = 0.6) +
  guides(color = "none") +
  scale_color_brewer(palette = "Set1") +
  theme_classic()

Wir sehen, dass die Verteilung der empirischen Mittelwerte etwas breiter ist als die der geschätzten (“geschrumpften”) Mittelwerte. Im Falle ungleich grosser Unternehmensstichproben wäre die Shrinkage/Schrumpfung noch ausgeprägter, wobei die Mittelwerte von Firmen mit kleineren Stichproben stärker zum Gesamtmittelwert “gezogen” werden würden als die Mittelwerte von Firmen mit grossen Stichproben (was die geringere Genauigkeit der Schätzungen aus kleinen Stichproben widerspiegelt).

1.4 Random-Intercept-Modell (Modell 2)

Dieses Modell enthält einen Level-1-Prädiktor (experience), der jedoch als fester Effekt konzeptualisiert wird (kein random slope, nur random intercept). Dieses Modell wird insbesondere benötigt, um den Anteil der durch den Level-1-Prädiktor erklärten Varianz zu ermitteln (vgl. Eid, Gollwitzer, und Schmitt 2017).

Level-1-Modell:
\(\boldsymbol {y_{mi}=\beta_{0i} + \beta_{1i} \cdot x_{mi} + \varepsilon_{mi}}\)

Level-2-Modell:
\(\boldsymbol {\beta_{0i}=\gamma_{00} + \upsilon_{0i}}\)
\(\boldsymbol {\beta_{1i}=\gamma_{10}}\)

Gesamtmodell:
\(\boldsymbol {y_{mi}=\gamma_{00} + \gamma_{10} \cdot x_{mi} + \upsilon_{0i}+\varepsilon_{mi}}\)

1.4.1 Modell berechnen

# Fitten und Abspeichern des Modells.
# Wenn eine Prädiktorvariable im Modell ist,
# muss der Intercept (also 1) nicht zusätzlich angegeben werden.
random.intercept.model <- lmer(salary ~ experience + (1 | company),
  data = df, REML = TRUE
)

# Abspeichern von Prädiktionen, die dieses Modell macht.
# Wir speichern diese Prädiktionen (vorhergesagten Werte),
# um das Modell später mit "ggplot2" zu veranschaulichen.
df$random.intercept.preds <- predict(random.intercept.model)


# Modell output anschauen
summary(random.intercept.model)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: salary ~ experience + (1 | company)
   Data: df

REML criterion at convergence: 1865.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.8109 -0.6884  0.0005  0.5980  3.8833 

Random effects:
 Groups   Name        Variance Std.Dev.
 company  (Intercept) 0.6144   0.7838  
 Residual             1.1845   1.0883  
Number of obs: 600, groups:  company, 20

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)   5.96418    0.22941  48.00060   26.00   <2e-16 ***
experience    0.53434    0.02721 589.48148   19.64   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
           (Intr)
experience -0.615

Die Ergebnisse zeigen, dass der Fixed-Effect von experience positiv signifikant ist \((\hat\gamma_{10} = 0.534, p < 0.001)\). Im Durchschnitt aller Firmen steigt das vorhergesagte Einkommen (salary) also mit jedem Jahr Arbeitserfahrung um ca. 0.534 Einheiten (534 CHF) an.

Mit der Funktion ranef() kann man sich wieder die einzelnen Random-Effects \(\upsilon_{0i}\) ausgeben lassen:

ranef(random.intercept.model)

$company
           (Intercept)
Company 01  0.20430372
Company 02  0.64614732
Company 03  1.49200151
Company 04 -0.91078990
Company 05  0.38916512
Company 06 -0.92463977
Company 07  0.57766959
Company 08 -0.51651767
Company 09 -0.63824646
Company 10  0.76848113
Company 11 -0.61955111
Company 12  1.09133530
Company 13 -0.77367207
Company 14 -0.73817926
Company 15 -0.65294087
Company 16  0.05733923
Company 17  0.45805487
Company 18 -0.08938416
Company 19  0.94422822
Company 20 -0.76480474

with conditional variances for "company"

Die Werte scheinen im Vergleich zum Intercept-Only-Modell (s.o.) deutlich kleiner geworden zu sein. Was das bedeutet, sehen wir gleich bei der Berechnung der erklärten Gesamtvarianz durch den Prädiktor experience.

1.4.2 Varianzerklärung durch Level-1-Prädiktor

Wieviel Level-1-Varianz wurde durch experience erklärt?

\(R^2_{Level-1} = \frac {\sigma^2_{\varepsilon1} - \sigma^2_{\varepsilon2}}{\sigma^2_{\varepsilon1}} = \frac {{1.9547}-{1.1845}}{1.9547} = {0.394}\)

Wieviel Gesamtvarianz wurde durch experience erklärt?

\(R^2_{Gesamt} = \frac {(\sigma^2_{\upsilon_01} + \sigma^2_{\varepsilon1}) - (\sigma^2_{\upsilon_02} + \sigma^2_{\varepsilon2})}{\sigma^2_{\upsilon_01} + \sigma^2_{\varepsilon1}} = \frac{({0.8512 + 1.9547}) - ({0.6144 + 1.1845}) } { 0.8512 + 1.9547} = 0.359\)

Wir sehen, dass durch Aufnahme des Level-1 Prädiktors experience neben der Level-1-Varianz auch die Level-2-Varianz kleiner wurde. Das liegt daran, dass sich die Firmen nicht nur systematisch bezüglich der AV salary unterscheiden, sondern auch bezüglich des Prädiktors experience. Wie in einer Kovarianzanalyse verändert sich der Effekt (bzw. die Level-2-Varianz) von company durch Einbezug einer mit dem Faktor korrelierten Kovariate. Mehr dazu werden wir in der nächsten Übung erfahren.

\(~\)

Code

df$random.intercept.preds <- predict(random.intercept.model)

ggplot(data = df, aes(
  x = experience,
  y = random.intercept.preds,
  colour = company
)) +
  geom_smooth(method = "lm", fullrange = TRUE, se = F, size = .3) +
  geom_jitter(aes(y = salary), alpha = .2) +
  labs(x = xlabel, y = ylabel) +
  ggtitle("Random-Intercept-Modell") +
  geom_abline(
    intercept = fixef(random.intercept.model)["(Intercept)"],
    slope = fixef(random.intercept.model)["experience"]
  ) +
  scale_colour_discrete() +
  ggthemes::theme_tufte()

1.5 Random-Coefficients-Modell (Modell 3)

Dieses Modell enthält sowohl einen Random-Intercept, als auch einen Random-Slope. Mit Hilfe dieses Modells kann die Varianz des Level-1-Slopes geschätzt werden. Damit erhalten wir ein Mass für die Unterschiedlichkeit des Effekts (Slope) des Level-1-Prädiktors (X, experience) auf die AV (Y, salary) zwischen den Level-2-Einheiten (company).

Level-1-Modell:
\(\boldsymbol {y_{mi}=\beta_{0i} + \beta_{1i} \cdot x_{mi} + \varepsilon_{mi}}\)

Level-2-Modell:
\(\boldsymbol {\beta_{0i}=\gamma_{00} + \upsilon_{0i}}\)
\(\boldsymbol {\beta_{1i}=\gamma_{10} + \upsilon_{1i}}\)

Gesamtmodell:
\(\boldsymbol {y_{mi}=\gamma_{00} + \gamma_{10} \cdot x_{mi} + \upsilon_{0i} + \upsilon_{1i} \cdot x_{mi} + \varepsilon_{mi}}\)

1.5.1 Modell berechnen

# Fitten und Abspeichern des Modells.
# Wenn für eine Prädiktorvariable ein Random-Slope im Modell ist, muss der
# Intercept (also 1) nicht mehr zusätzlich im Klammerausdruck für die zufälligen
# Effekte angegeben werden.
random.coefficients.model <- lmer(salary ~ experience + (experience | company),
  data = df, REML = TRUE
)

# Abspeichern von Prädiktionen, die dieses Modell macht.
# Wir speichern diese Prädiktionen, um das Modell später mit ggplot2
# zu veranschaulichen.
df$random.coefficients.preds <- predict(random.coefficients.model)


# Modell output anschauen
summary(random.coefficients.model)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: salary ~ experience + (experience | company)
   Data: df

REML criterion at convergence: 1855.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.8307 -0.6804  0.0037  0.5999  3.3608 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 company  (Intercept) 0.72257  0.8500        
          experience  0.01839  0.1356   -0.51
 Residual             1.13629  1.0660        
Number of obs: 600, groups:  company, 20

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  5.93372    0.24038 18.88905   24.68 7.77e-16 ***
experience   0.53085    0.04059 18.95008   13.08 6.20e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
           (Intr)
experience -0.690

Die Ergebnisse zeigen, dass der Fixed-Effect von experience wiederum positiv signifikant ist \((\hat\gamma_{10} = 0.531, p < 0.001)\) (leicht unterschiedlich im Vergleich zum Random-Intercept-Modell). Das Hauptinteresse im Random-Coefficients-Modell gilt aber der Varianz bzw. Standardabweichung des Random-Slope: \(\hat\sigma^2_{\upsilon_1} = 0.0184\) und damit \(\hat\sigma_{\upsilon_1} = \sqrt{0.0184} = 0.1356\). Die durchschnittliche Abweichung (Level-2-Residuum \(\upsilon_{1i}\)) des experience-Effekts einer Firma vom mittleren Effekt \(\hat\gamma_{10}\) beträgt also 0.1356 Einheiten (135.6 CHF).

Mit der Funktion ranef() kann man sich wieder die einzelnen Random-Effects (jetzt Level-2-Residuen des Intercepts \(\upsilon_{0i}\) und des Slopes \(\upsilon_{1i}\)) ausgeben lassen:

ranef(random.coefficients.model)

$company
            (Intercept)   experience
Company 01 -1.003794102  0.203747776
Company 02 -0.160820541  0.143400958
Company 03  0.705132570  0.139886714
Company 04 -0.611829652 -0.054018876
Company 05 -0.033975063  0.076641659
Company 06 -0.123906149 -0.150704419
Company 07  0.669506844 -0.013616984
Company 08 -0.830781477  0.065813067
Company 09 -1.020757779  0.086008002
Company 10  0.392944052  0.081964065
Company 11 -0.769063062  0.037944670
Company 12  1.248886366 -0.025028992
Company 13 -0.593370830 -0.025453610
Company 14 -0.186191134 -0.109527150
Company 15  0.007222191 -0.150163600
Company 16  0.772829070 -0.140812590
Company 17  0.562032356 -0.011037598
Company 18  0.508160168 -0.112469594
Company 19  1.006168875 -0.006676886
Company 20 -0.538392700 -0.035896609

with conditional variances for "company"

Im Output befindet sich auch eine Angabe zum dritten Random-Effect-Parameter, der in diesem Modell geschätzt wurde, nämlich zur Kovarianz \(\hat\sigma_{\upsilon_0\upsilon_1}\) zwischen den Level-2-Residen des Intercepts und den Level-2-Residuen des Slopes. Allerdings wird nicht die Kovarianz selbst ausgegeben, sondern ihre standardisierte Variante, die Korrelation \(r_{\upsilon_0\upsilon_1} = -0.51\). Diese besagt, dass das durchschnittliche Gehalt eines Berufsanfängers in einer Firma (vorhergesagter Wert von salary an der Stelle experience = 0) negativ mit der Steigung von experience korreliert ist. Zur Interpretation: Wenn schon Berufsanfänger in einer Firma gut verdienen, gibt es offenbar weniger Spielraum für Gehaltssteigerungen als in Firmen mit weniger hohem Einstiegsgehalt.

1.5.2 Signifikanztest für \(\sigma^2_{\upsilon_1}\): Modellvergleich mit Random-Intercept-Modell

Wenn wir dieses Modell mittels LR-Test mit dem Random-Intercept-Modell (s.o.) vergleichen, werden zwei Parameter gleichzeitig getestet: die Level-2-Varianz des Slopes (\(\sigma^2_{\upsilon_1}\)) und die Level-2-Kovarianz zwischen Intercept und Slope (\(\sigma_{\upsilon_0\upsilon_1}\)). Laut gegenwärtigem Forschungsstand ist dies die beste Methode, um zu überprüfen, ob es eine signifikante Slope-Varianz gibt (vgl. Eid, Gollwitzer, und Schmitt 2017).

1.5.2.1 Modellvergleich von Hand

Zunächst benötigen wir die -2-Log-Likelihoods (Devianzen) der beiden zu vergleichenden Modelle. Diese wird im Modelloutput des jeweiligen Modells unter der Bezeichnung “REML criterion at convergence” ausgegeben.

Wir folgen hier ausserdem den Bezeichnungen für uneingeschränkte Modelle \((M_u)\) und für eingeschränkte (restringierte) Modelle \((M_e)\), wie wir sie auch für die Modellvergleiche bei der logistischen Regression in Statistik III verwendet haben. \(M_u\) steht hier für das Random-Coefficients-Modell, das die zu testenden Parameter (Slope-Varianz \(\sigma^2_{\upsilon_1}\) und Intercept-Slope-Kovarianz \(\sigma_{\upsilon_0\upsilon_1}\)) enthält, \(M_e\) für das Random-Intercept-Model, das genau diese Parameter nicht enthält.

-2-Log-Likelihood random.intercept.model (\(M_e\)): \(-2\cdot ln(L_e) = 1865.4\)

-2-Log-Likelihood random.coefficients.model (\(M_u\)): \(-2\cdot ln(L_u) = 1855.4\)

\(~\)

Berechnen des Modellvergleichstests \(LR(M_u - M_e)\):

\(\chi^2_{emp} = [-2 \cdot ln(L_e)] - [-2 \cdot ln(L_u)] = 1865.4 - 1855.4 = 10.0\)

Anzahl Parameter \((nPar) \space M_u: \space \textbf{6} \space \space (\gamma_{00},\space\gamma_{10},\space\sigma^2_{\varepsilon}, \space\sigma^2_{\upsilon_0}, \space\sigma^2_{\upsilon_1}, \space\sigma_{\upsilon_0\upsilon_1})\)

Anzahl Parameter \((nPar) \space M_e:\space \textbf{4} \space \space (\gamma_{00}, \space\gamma_{10}, \space\sigma^2_{\varepsilon}, \space\sigma^2_{\upsilon_0})\)

\(LR(M_u - M_e)\) ist \(\chi^2\)-verteilt mit \(df = nPar(M_u) - nPar(M_e) = 6 - 4 = 2\)

Kritischer Wert: \(\chi^2_{(0.95, df = 2)} = 5.9915\)

Aber: Es handelt sich genau genommen um eine Mischverteilung aus \(\chi^2_{df=1}\) und \(\chi^2_{df=2}\), da einer der beiden getesteten Parameter \((\sigma^2_{\upsilon_1})\) einen bounded parameter space aufweist (kann nicht kleiner als 0 werden). Daher ist der korrekte kritische Chi-Quadrat-Wert hier das 95 %-Quantil dieser Mischverteilung und liegt bei (gerundet) 5.14 (vgl. Eid, Gollwitzer, und Schmitt 2017).

Der Test ist signifikant, da \(\chi^2_{emp} \space(=10.0) > \chi^2_{krit} \space(=5.14)\).

Beide Parameter zusammen unterscheiden sich also signifikant von 0. Der Fokus der Interpretation liegt hier aber auf dem Parameter der Slope-Varianz \(\sigma^2_{\upsilon_1}\) (nur wenn dieser > 0 ist, kann auch die Intercept-Slope-Kovarianz \(\sigma_{\upsilon_0\upsilon_1}\) überhaupt ungleich 0 sein).

Inhaltlich bedeutet dies, dass sich die Stärke des Effekts von experience auf salary signifikant zwischen den Firmen (company) unterscheidet. Die oben im Random-Coefficients-Modell geschätzte Slope-Varianz von experience \((\hat \sigma^2_{\upsilon_1} = 0.0184)\) ist damit als signifikant zu bewerten.

1.5.2.2 Modellvergleich in R

Wie beim Signifikanztest der Random-Intercept-Varianz \(\sigma^2_{\upsilon_0}\) im Intercept-Only-Modell können wir den Modellvergleich hier wieder mit ranova() vornehmen. Wie oben wird ranova() auf das uneingeschränkte Modell angewandt, die Angabe eines (eingeschränkten) Vergleichsmodells ist nicht notwendig.

ranova(random.coefficients.model)

ANOVA-like table for random-effects: Single term deletions

Model:
salary ~ experience + (experience | company)
                                     npar  logLik    AIC    LRT Df Pr(>Chisq)
<none>                                  6 -927.70 1867.4                     
experience in (experience | company)    4 -932.71 1873.4 10.011  2   0.006699
                                       
<none>                                 
experience in (experience | company) **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

\(~\)

Alternativ können wir den Modellvergleich auch mit der uns bereits bekannten Funktion anova() durchführen, indem wir die beiden zu vergleichenden Modelle explizit angeben:

anova(random.coefficients.model, random.intercept.model, refit = FALSE)

Data: df
Models:
random.intercept.model: salary ~ experience + (1 | company)
random.coefficients.model: salary ~ experience + (experience | company)
                          npar    AIC    BIC  logLik deviance  Chisq Df
random.intercept.model       4 1873.4 1891.0 -932.71   1865.4          
random.coefficients.model    6 1867.4 1893.8 -927.70   1855.4 10.011  2
                          Pr(>Chisq)   
random.intercept.model                 
random.coefficients.model   0.006699 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Vertiefung

Bei der anova() der beiden mit lmer() geschätzten Modelle werden in der Default-Einstellung (refit = TRUE) beide Modelle vor Durchführung des LR-Tests automatisch nochmals mittels ML (statt REML) geschätzt (falls vorher REML verwendet wurde) und der LR-Test dann auf die ML-gefitteten Modellobjekte angewendet. Dies ist eine Vorsichtsmassnahme der Entwickler des Packages, um den häufigen Fehler zu vermeiden, REML-geschätzte Modelle, die sich (auch) in ihrer Fixed-Effects-Struktur voneinander unterscheiden, per LR-Test zu vergleichen. Da wir hier zwei Modelle vergleichen, die sich nur in ihrer Random-Effects-Struktur unterscheiden, wählen wir refit = FALSE, um den Vergleich der mit REML gefitteten Modelle zu erhalten. Beim Modellvergleich mit ranova() spielt das keine Rolle, da dort ein mit REML bzw. mit ML gefittetes Modell automatisch mit dem ensprechend gefitteten reduzierten Modell verglichen wird. Der Grund dafür ist, dass ranova() nur Modellvergleiche zum Test von Random-Effects durchführt und die Refit-Problematik daher dort nicht existiert.

Das Ergebnis zeigt einen signifikanten LR-Test mit \(\chi^2_{emp} = 10.011\). Der (Rundungs-)Unterschied zum Modellvergleich von Hand ergibt sich wegen der ungenaueren Angabe der Devianzen in den Modelloutputs.

Hinweis

Man kann den empirischen Chi-Quadrat-Wert des Modellvergleichs auch mit den Angaben zu den Log-Likelihoods (logLik) der Modelle aus dem Output von ranova() bzw. anova() nachrechnen:

\(\chi^2_{emp} = [-2 \cdot (-932.71)] - [-2 \cdot (-927.70)] = 1865.42 - 1855.40 = 10.02\)

Hier ergibt sich dann eine Ungenauigkeit in die andere Richtung…

Der angegebene \(p\)-Wert ist allerdings noch nicht ganz korrekt, da einer der beiden Parameter - \(\sigma^2_{\upsilon_1}\) - wie oben bereits festgestellt einen bounded parameter space aufweist. Der korrekte \(p\)-Wert auf Basis der Mischverteilung ist daher noch etwas kleiner als der ausgegebene.

Vertiefung: \(p\)-Wert der Mischverteilung

Die exakte Bestimmung des p-Werts aufgrund der Mischverteilung erfolgt über die Mittelung der p-Werte der beiden Verteilungen \(\chi^2_{df=1}\) und \(\chi^2_{df=2}\):

0.5 * pchisq(10.011, 1, lower.tail = FALSE) +
  0.5 * pchisq(10.011, 2, lower.tail = FALSE)

[1] 0.004128535

\(~\)

Code

# sal.visualisation$random.intercept.random.slope.plot
df$random.coefficients.preds <- predict(random.coefficients.model)

ggplot(data = df, aes(x = experience, y = random.coefficients.preds, colour = company)) +
  geom_smooth(method = "lm", fullrange = TRUE, se = F, size = .3) +
  labs(x = xlabel, y = ylabel) +
  geom_jitter(aes(y = salary), alpha = .2) +
  geom_abline(
    intercept = fixef(random.coefficients.model)["(Intercept)"],
    slope = fixef(random.coefficients.model)["experience"]
  ) +
  ggtitle("Random Coefficients Model") +
  scale_colour_discrete() +
  ggthemes::theme_tufte()

1.6 Übung

Wir möchten (fiktive) Daten zur Pisa-Studie analysieren. Dafür haben wir einen Datensatz mit den Variablen:

id : Durchnummerierung aller Schüler (1:1000)

name : Name der Schülerin / des Schülers

grade : Die durchschnittliche Vornote jedes Schülers (gerundet auf 0.25)

motivation : Wie motiviert die Schüler für den Pisa-Test sind

class : Variable der Schulklasse (durchnummeriert, 50 Klassen wurden zufällig ausgewählt)

teacher_salary : Monatslohn des Lehrers (jede Klasse hat einen anderen Lehrer)

school : Kategoriale Variable der Schule (10 Schulen wurden zufällig ausgewählt)

pisa : Erreichte Punktzahl beim Pisa-Test.

Aus jeder der 10 Schulen wurden 5 Klassen ausgesucht. Insgesamt sind es also 50 Klassen. Für die folgenden Aufgaben beachten wir die Nestung der Klassen in den verschiedenen Schulen nicht, da wir dafür ein 3-Ebenen-Modell bräuchten. Unser Sample auf Level-2 sind also die 50 Schulklassen. Die Daten kann man sich etwa so vorstellen (subsample, nur 15 Klassen für die Visualisierung):

Code

school_df <- read_csv(url("https://raw.githubusercontent.com/methodenlehre/data/master/statistik_IV/simulated_school_data.csv"))

subsample <- school_df[1:300, ]
subsample <- subsample |> mutate(class = as.factor(class))

ggplot(data = subsample, aes(x = grade, y = pisa, group = class, color = class)) +
  geom_smooth(method = "lm", alpha = .5, se = FALSE, show.legend = F, fullrange = TRUE, size = .5) +
  geom_point(alpha = .3) +
  ggthemes::theme_tufte()

1. Aufbereitung der Daten

Laden Sie die Daten in R und schauen Sie sich die ersten paar Zeilen des Datensatzes an. Die Daten können entweder zuerst heruntergeladen und im Ordner versorgt werden, oder man kann sie auch direkt von R aus dem Internet laden (weil sie öffentlich im Internet zugänglich sind). Die Daten sind hier:

https://raw.githubusercontent.com/methodenlehre/data/master/statistik_IV/simulated_school_data.csv

Tipp: Verwenden Sie dieselbe Funktion wie weiter oben im Skript um Daten direkt aus dem Internet zu laden.

Lösung

# 1a
pacman::p_load(lme4, nlme, tidyverse, lmerTest, gridExtra, ggplot2)

school_df <- read_csv(url("https://raw.githubusercontent.com/methodenlehre/data/master/statistik_IV/simulated_school_data.csv"))

Rows: 1150 Columns: 8
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): name, school
dbl (6): id, motivation, grade, class, teacher_salary, pisa

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(school_df)

# A tibble: 6 × 8
     id name   motivation grade class teacher_salary school   pisa
  <dbl> <chr>       <dbl> <dbl> <dbl>          <dbl> <chr>   <dbl>
1     1 Denis        4        2     1           4.74 Bitzius   513
2     2 Talia        5.25     4     1           4.74 Bitzius   563
3     3 Ahmed        6        4     1           4.74 Bitzius   610
4     4 Lamija       4.75     4     1           4.74 Bitzius   337
5     5 Titus        4.5      4     1           4.74 Bitzius   477
6     6 Justus       2        1     1           4.74 Bitzius   456

Stellen Sie sicher, dass die Variablen im richtigen Format sind. Als Tipp: Überlegen Sie sich, welche Variablen kategorial und welche numerisch sind.

Lösung

# 1b
school_df <- school_df |> mutate(
  id = as.factor(id),
  class = as.factor(class)
)

(optional) Visualisieren Sie die Daten mit der Funktion plot(). Als Tipp: Wenn nur numerische Variablen im Datenframe sind, können mit plot(datenframe) direkt die Zusammenhänge der Variablen dargestellt werden.

Lösung

# 1c
plot(school_df[, -c(1:2, 5, 7)])

# entfernt die kategorialen Columns (1, 2, 5 und 7) und plottet den Rest.

Man sieht hier keine offensichtlichen Korrelationen.

2. Fragestellung und Intercept-Only-Modell

Wir gehen davon aus, dass die Schulklassen zufällig ausgesucht wurden. Wir sind nicht daran interessiert, wie sich diese spezifischen Klassen voneinander unterscheiden. Uns interessiert vielmehr, ob die Vornote (grade) sowie die Motivation (motivation) einen Einfluss auf die Pisa-Werte haben. Die Schulklassen möchten wir als Random-Effect modellieren.

Berechnen Sie ein Intercept-Only-Modell.

Lösung

# 2a
intOnly_fit <- lmer(
  formula = pisa ~ 1 + (1 | class),
  data = school_df
)
summary(intOnly_fit)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: pisa ~ 1 + (1 | class)
   Data: school_df

REML criterion at convergence: 13768.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.1991 -0.5967  0.0214  0.6522  2.9179 

Random effects:
 Groups   Name        Variance Std.Dev.
 class    (Intercept) 12559    112.07  
 Residual              7983     89.35  
Number of obs: 1150, groups:  class, 50

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)    
(Intercept)   507.52      16.07  49.00   31.59   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Testen Sie, ob der Random-Effect (class) in diesem Modell tatsächlich Varianz aufweist. Ist die Modellstruktur gerechtfertigt, d.h. benötigen wir überhaupt ein HLM-Modell?

Lösung

# 2b
ranova(intOnly_fit)

ANOVA-like table for random-effects: Single term deletions

Model:
pisa ~ (1 | class)
            npar  logLik   AIC    LRT Df Pr(>Chisq)    
<none>         3 -6884.4 13775                         
(1 | class)    2 -7332.0 14668 895.23  1  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Die Klassen unterscheiden sich signifikant im Pisa-Durchschnittswert. Man erkennt das daran, dass der \(p\)-Wert der Varianz für den random-effect kleiner als 0.05 ist. Wir können diesen \(p\)-Wert sogar noch halbieren (s.o.). Die Modellstruktur ist somit also gerechtfertigt.

Wie gross ist die Intraklassenkorrelation?

Lösung

Um die Intraklassenkorrelation zu berechnen, brauchen wir die Varianzen des intOnly_fit Objekts. Statt die Random-Intercept-Varianz und die Level-1-Residualvarianz aus dem Output abzulesen, kann man Sie auch aus dem Modell-Objekt extrahieren und abspeichern.

# Varianz von intOnly_fit extrahieren
intOnly_class_var <- as.data.frame(VarCorr(intOnly_fit))$vcov

classVar <- intOnly_class_var[1] # das erste Element ist die Random-Intercept-Varianz
resVar <- intOnly_class_var[2] # das zweite Element ist die Level-1-Residualvarianz
totalVar <- classVar + resVar # die Summe aus beiden ist die Gesamtvarianz

intraclass_correlation <- (classVar / totalVar) |> round(3)

intraclass_correlation

[1] 0.611

Die Intraclass Correlation ist \(\hat{\rho} =\) 0.611. Klasse erklärt also ca 61.1% der gesamten Varianz der Pisa-Werte in diesem Modell.

Statt die Varianzen so wie oben aus dem Output-Objekt intOnly_fit zu indizieren, können sie auch aus dem Output (s.o. 2a) herausgelesen werden: \[\boldsymbol {\hat \rho = \frac{\hat \sigma^2_{Level-2}}{\hat \sigma^2_{Level-2} + \hat \sigma^2_{Level-1}}} = \frac{1.25587\times 10^{4}}{1.25587\times 10^{4} + 7983.4} = 0.611\].

3. Random-Intercept-Modell

Berechnen Sie ein Random-Intercept-Modell mit den beiden Prädiktoren grade und motivation. Sind die fixed effects (\(\hat\gamma_{10}\) und \((\hat\gamma_{20}\)) beider Prädiktoren signifikant?

Lösung

# 3a
randInt_fit <- lmer(
  formula = pisa ~ grade + motivation + (1 | class),
  data = school_df
)
summary(randInt_fit)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: pisa ~ grade + motivation + (1 | class)
   Data: school_df

REML criterion at convergence: 13699.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.9786 -0.6213  0.0302  0.6310  2.8518 

Random effects:
 Groups   Name        Variance Std.Dev.
 class    (Intercept) 12536    111.96  
 Residual              7555     86.92  
Number of obs: 1150, groups:  class, 50

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  440.336     20.073  118.333  21.937  < 2e-16 ***
grade         14.730      1.882 1100.532   7.826 1.18e-14 ***
motivation     4.890      2.355 1100.484   2.077   0.0381 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
           (Intr) grade 
grade      -0.309       
motivation -0.526  0.036

Beide Prädiktoren sind signifikant. An den Estimates kann man jedoch nicht direkt ablesen, wessen Einfluss auf den Pisa-Wert stärker ist, weil diese Werte von der Skalierung der Prädiktoren abhängen.

Testen Sie, ob ein zusätzlicher Interaktionseffekt zwischen grade und motivation ins Modell aufgenommen werden sollte. Obwohl wir keine spezifische Moderationshypothese haben, können wir im Sinne eines Vermeidens von Underfitting überprüfen, ob eine bedeutsame Interaktion vorhanden ist.

Lösung

# Wir können entweder das gesamte Modell neu definieren, neu auch mit
# grade:motivation, oder wir nutzen die Funktion update() und fügen
# einfach den Prädiktor grade:motivation hinzu.
# Beides führt zum selben Resultat.

randInt_interaction_fit <- lmer(
  formula = pisa ~ grade * motivation + (1 | class),
  data = school_df
)
summary(randInt_interaction_fit)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: pisa ~ grade * motivation + (1 | class)
   Data: school_df

REML criterion at convergence: 13694.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8938 -0.6279  0.0308  0.6441  2.8305 

Random effects:
 Groups   Name        Variance Std.Dev.
 class    (Intercept) 12507    111.8   
 Residual              7552     86.9   
Number of obs: 1150, groups:  class, 50

Fixed effects:
                 Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)       408.368     31.913  552.133  12.796  < 2e-16 ***
grade              24.842      8.074 1098.751   3.077  0.00215 ** 
motivation         12.009      6.009 1098.485   1.999  0.04589 *  
grade:motivation   -2.261      1.755 1098.712  -1.288  0.19808    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) grade  motvtn
grade       -0.802              
motivation  -0.845  0.898       
grade:mtvtn  0.778 -0.972 -0.920

Ein Interaktionseffekt scheint in diesem Modell keinen Sinn zu machen, der \(p\)-Wert für den Interaktionseffekt ist nicht signifikant \((p_{grade:motivation} = 0.1981)\), potentielle moderierende Einflüsse zwischen den beiden Prädiktoren sind also nicht auszumachen.

4. Random-Coefficients-Modell

Zusätzlich zum Random-Intercept können wir auch die Varianzen der Slopes von grade und motivation mit ins Modell aufnehmen. Doch “lohnt” sich das überhaupt, d.h. unterscheiden sich die Effekte der beiden Prädiktoren jeweils signifikant zwischen den Klassen?

Schätzen Sie ein Random-Coefficients Modell mit Random Slopes für grade und motivation.

Lösung

randCoef_fit <- lmer(formula = pisa ~ grade + motivation +
  (grade + motivation | class), data = school_df)

summary(randCoef_fit)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: pisa ~ grade + motivation + (grade + motivation | class)
   Data: school_df

REML criterion at convergence: 13599.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2014 -0.6149  0.0275  0.6347  2.7361 

Random effects:
 Groups   Name        Variance Std.Dev. Corr       
 class    (Intercept) 9519.1   97.57               
          grade        275.3   16.59    -0.10      
          motivation   513.2   22.65    -0.56  0.23
 Residual             6382.5   79.89               
Number of obs: 1150, groups:  class, 50

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)  430.401     18.178  46.240  23.677  < 2e-16 ***
grade         14.977      2.935  47.140   5.103 5.89e-06 ***
motivation     6.970      3.921  48.741   1.777   0.0817 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
           (Intr) grade 
grade      -0.255       
motivation -0.670  0.167

Im Random-Coefficients-Modell ist der Fixed Effect von motivation \((\hat\gamma_{20})\) nun nicht mehr signifikant \((p_{motivation} = 0.0817)\). Das ist nicht ungewöhnlich, die Schätzung von Random Effects Parametern beeinflusst auch diejenigen der Fixed Effects, und derjenige von motivation war im Random-Intercept-Modell nur knapp signifikant.

Testen Sie die Signifikanz beider Random Slopes mit ranova().

Lösung

Die ranova()-Funktion testet die Signifikanz von Level-2-Varianzparametern mithilfe von Likelihood-Ratio-Tests. In unserem Beispiel wird das randCoef_fit Modell mit Modellen verglichen, die jeweils eine Slope-Varianz (und deren Kovarianzen mit dem Intercept und dem jeweils anderen Slope) weniger im Modell haben. Ein signifikanter Test bedeutet hier, dass die Daten schlechter auf die eingeschränkten Vergleichsmodelle passen als auf das uneingeschränkte Modell randCoef_fit.

Die \(p\)-Werte des Likelihood-Ratio-Tests müssten noch angepasst werden, weil hier zwei “unbounded” (die Kovarianz des Slopes des jeweils getesteten Prädiktors mit dem random Intercept und die Kovarianz dieses Slopes mit dem Slope des jeweils anderen Prädiktors) und ein “bounded” Parameter (die jeweilige Slope-Varianz) gleichzeitig getestet werden. Darauf gehen wir jedoch hier nicht weiter ein.

ranova(randCoef_fit)

ANOVA-like table for random-effects: Single term deletions

Model:
pisa ~ grade + motivation + (grade + motivation | class)
                                           npar  logLik   AIC    LRT Df
<none>                                       10 -6799.8 13620          
grade in (grade + motivation | class)         7 -6822.0 13658 44.492  3
motivation in (grade + motivation | class)    7 -6828.7 13671 57.895  3
                                           Pr(>Chisq)    
<none>                                                   
grade in (grade + motivation | class)       1.186e-09 ***
motivation in (grade + motivation | class)  1.655e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Beide Slope-Varianzen sind signifikant: Sowohl das Weglassen von grade in (grade + motivation | class) als auch das Weglassen von motivation in (grade + motivation | class) würde zu einem signifikant schlechteren Fit führen (signifikanter Likelihood-Ratio-Test LRT). Die Klassen unterscheiden sich also signifikant bezüglich der Effekte von motivation und grade. Ein Random-Coefficients-Modell ist daher notwendig und erklärt die Daten besser als das Random-Intercept-Modell.

Der durchschnittliche (fixed) Effekt von motivation über alle Klassen hinweg ist zwar nicht signifikant von 0 unterschieden (s.o.), trotzdem ist der Einfluss der motivation auf die pisa-Werte aber zwischen den Klassen verschieden. Möglicherweise gibt es Merkmale auf Klassenebene (wie z.B. Klassengrösse, um ein Beispiel zu nennen), die erklären könnten, in welchen Klassen motivation einen Effekt auf die pisa-Werte hat und in welchen nicht. Darauf gehen wir aber nicht weiter ein, auch weil wir (ausser teacher_salary) keine Merkmale auf Klassenebene im Datensatz haben.

Vertiefung

Basierend auf dem besten Modell: Was ist der vorhergesagte Wert eines Schülers mit…
- grade = 4.5, motivation = 1 (aus einer unbekannten Klasse)
- grade = 4, motivation = 6, class = 18

Lösung

# 1
(schüler1 <- fixef(randCoef_fit)[1] + 4.5 * fixef(randCoef_fit)[2] +
  1 * fixef(randCoef_fit)[3])

(Intercept) 
   504.7657

# Überprüfen mit predict()
predict(randCoef_fit, newdata = tibble(grade = 4.5, motivation = 1), re.form = NA)

       1 
504.7657

# 2
(schüler2 <- fixef(randCoef_fit)[1] + 4 * fixef(randCoef_fit)[2] +
  6 * fixef(randCoef_fit)[3] + ranef(randCoef_fit)$class$`(Intercept)`[18] +
  4 * ranef(randCoef_fit)$class$`grade`[18] + 6 *
    ranef(randCoef_fit)$class$`motivation`[18])

(Intercept) 
   821.5456

# Überprüfen mit predict()
predict(randCoef_fit, newdata = tibble(grade = 4, motivation = 6, class = 18))

       1 
821.5456

Eid, Michael, Mario Gollwitzer, und Manfred Schmitt. 2017. Statistik und Forschungsmethoden: mit Online-Materialien. 5., korrigierte Auflage. Weinheim Basel: Beltz.

Gelman, Andrew, und Jennifer Hill. 2007. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.