2  Modelle mit Level-2-Prädiktoren

2.1 Setup

Wir laden alle Packages mit pacman.

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

Wir laden unsere Daten und definieren company und die neue Variable sector als Faktoren. sector besitzt die beiden Faktorstufen Private und Public, von denen Private als alphabetisch frühere Kategorie automatisch als erste und damit als Referenzkategorie gesetzt wird.

df <- read_csv(
  url("https://raw.githubusercontent.com/methodenlehre/data/master/statistik_IV/salary-data.csv")
) |>
  mutate(
    company = as.factor(company),
    sector = as.factor(sector)
  )

Jetzt können wir uns die Daten anschauen.

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

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                                                  

2.2 Intercept-as-Outcome-Modell (Modell 4)

2.2.1 Kontexteffekt-Modell

Das Kontexteffekt-Modell ist eine Variante des Intercept-as-Outcome-Modells. Der Level-2-Prädiktor ist hier ein auf Gruppenebene aggregierter Level-1-Prädiktor (hier: exp_groupmean; diese Variable beinhaltet für jede Person den zugehörigen Firmenmittelwert von experience. Um die Einflüsse von Level-1- und Level-2-Prädiktor zu trennen, muss der Level-1-Prädiktor in diesem Modell als gruppenzentrierte Variable aufgenommen werden. Auf Level-1 wird also nur noch die firmenspezifische Abweichung von experience vom jeweiligen exp_groupmean modelliert, die gruppenzentierte Variable nennen wir daher exp_centered.

In unserem Beispiel ist der Level-2-Prädiktor die durchschnittliche Erfahrung pro Firma. Der Effekt dieser Variablen wird gemeinsam mit der gruppenzentrierten Erfahrung untersucht. Zusätzlich zu der Annahme, dass die individuelle Erfahrung mit einem höheren Lohn einhergeht, nehmen wir nun an, dass auch die durchschnittliche Erfahrung der Mitarbeiter einer Firma sich auf den durchschnittlichen Lohn einer Firma auswirkt. In diesem Beispiel ist die Annahme naheliegend, dass beide Effekte vorhanden sind und in die gleiche Richtung gehen. Das ist aber nicht immer der Fall (vgl. “Big-Fish-Little-Pond-Effekt”).

Im Kontexteffekt-Modell wird manchmal kein random slope zugelassen, weil man sich auf die fixed effects auf beiden Untersuchungsebenen konzentriert. Man könnte aber im Prinzip auch hier einen random slope von exp_centered zulassen. Wir verzichten hier darauf und verweisen auf das Allgemeine Intercept-as-Outcome-Modell weiter unten.

\(~\)

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

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

\(\boldsymbol {\beta_{1i}=\gamma_{10}}\)

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

\(~\)

Der Level-2-Prädiktor exp_groupmean muss jedoch zuerst berechnet und jeder Person zugeteilt werden. Wir berechnen diese Variable mit der Pipe |> und benutzen group_by() und mutate() aus dem tidyverse-Package. Gleichzeitig berechnen wir eine gruppenzentrierte Variable exp_centered, indem wir für jede Person ihren jeweiligen Firmenmittelwert von der individuellen Erfahrung abziehen.

# Berechnen der Mittelwerte pro Gruppe (exp_groupmean),
# und der Gruppenmittelwert-zentrierten experience (exp_centered)
df <- df |>
  group_by(company) |>
  mutate(
    exp_groupmean = mean(experience),
    exp_centered = experience - exp_groupmean
  ) |>
  ungroup()
# select() zum Ändern der Variablenreihenfolge

# Jetzt können wir uns die Daten nochmals anschauen - hat alles funktioniert?
head(df)

# A tibble: 6 × 6
  company    experience salary sector  exp_groupmean exp_centered
  <fct>           <dbl>  <dbl> <fct>           <dbl>        <dbl>
1 Company 01       6.87   8.96 Private          6.37        0.496
2 Company 01       5.66   9.54 Private          6.37       -0.714
3 Company 01       4.58   7.30 Private          6.37       -1.79 
4 Company 01       6.56   8.09 Private          6.37        0.186
5 Company 01       8.83  14.3  Private          6.37        2.46 
6 Company 01       7.73  10.3  Private          6.37        1.36 

tail(df)

# A tibble: 6 × 6
  company    experience salary sector  exp_groupmean exp_centered
  <fct>           <dbl>  <dbl> <fct>           <dbl>        <dbl>
1 Company 20       3.58   6.84 Private          5.04        -1.46
2 Company 20       3.18   7.60 Private          5.04        -1.86
3 Company 20       3.39   5.71 Private          5.04        -1.65
4 Company 20       7.12  10.1  Private          5.04         2.08
5 Company 20       2.98   6.94 Private          5.04        -2.06
6 Company 20       6.45   9.33 Private          5.04         1.41

# exp_groupmean anschauen
df |>
  group_by(company) |>
  summarise(exp_groupmean = mean(exp_groupmean)) |>
  ungroup()

# A tibble: 20 × 2
   company    exp_groupmean
   <fct>              <dbl>
 1 Company 01          6.37
 2 Company 02          6.13
 3 Company 03          6.09
 4 Company 04          4.71
 5 Company 05          6.34
 6 Company 06          5.1 
 7 Company 07          4.02
 8 Company 08          5.71
 9 Company 09          5.15
10 Company 10          5.21
11 Company 11          5.40
12 Company 12          4.91
13 Company 13          5.23
14 Company 14          4.70
15 Company 15          4.15
16 Company 16          4.88
17 Company 17          5.61
18 Company 18          5.03
19 Company 19          4.02
20 Company 20          5.04

Modell berechnen

# Fitten und abspeichern des Modells
context.model <- lmer(salary ~ exp_centered + exp_groupmean + (1 | company),
  data = df, REML = TRUE
)

# Durch das Modell vorhergesagte Werte abspeichern
df$context.model.preds <- predict(context.model)

# Modell anschauen mit summary()
summary(context.model)

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

REML criterion at convergence: 1865.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7991 -0.6796  0.0043  0.6008  3.8752 

Random effects:
 Groups   Name        Variance Std.Dev.
 company  (Intercept) 0.623    0.7893  
 Residual             1.185    1.0884  
Number of obs: 600, groups:  company, 20

Fixed effects:
               Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)     4.77925    1.38736  18.00000   3.445  0.00289 ** 
exp_centered    0.53188    0.02735 579.00000  19.446  < 2e-16 ***
exp_groupmean   0.76264    0.26499  18.00000   2.878  0.01001 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) exp_cn
exp_centerd  0.000       
exp_groupmn -0.991  0.000

Der Effekt von exp_centered ist mit \(\hat\gamma_{10} = 0.532\) sehr ähnlich wie im Random-Intercept-Modell. Der Effekt von exp_groupmean ist noch grösser - jedes Jahr durchschnittliche Erfahrung, die die Angestellten einer Firma mehr haben als die Angestellten einer anderen Firma, lässt den vorhergesagten Salary-Wert der Angestellten dieser Firma (und damit den Intercept \(\hat\beta_{0i}\)) um \(\hat\gamma_{01} = 0.763\) Einheiten (ca. 763 CHF) steigen.

Visualisierung des Modells

Die folgenden Plots spielen sich nur auf der Firmenebene (Level-2) ab. Der erste Plot zeigt die Abweichungen (Level-2-Residuen \(\upsilon_{0i}\)) der einzelnen Firmen-Intercepts vom Gesamtmittelwert der Intercepts (\(\hat\gamma_{00}\)) im Intercept-Only-Modell. Die Durchschnittsgehälter der Firmen unterscheiden sich stark vom Mittelwert der Gehälter aller Firmen (hellblaue Punkte). Zur Veranschaulichung sind ausserdem die (ohne Shrinkage geschätzten) empirischen salary-Mittelwerte pro Firma im Plot dargestellt (kleine schwarze Punkte).

Beim zweiten Plot kommt der Level-2 Prädiktor exp_groupmean (also die durchschnittliche Erfahrung, die die Angestellten jeder Firma haben) ins Spiel. Es zeigt sich, dass sich die exp_groupmean zwischen den Firmen unterscheidet, und dass exp_groupmean positiv mit den Intercepts des Intercept-Only-Modells korreliert (Steigung der Geraden = \(\hat\gamma_{01} = 0.763\)). Durch den Prädiktor nimmt die Abweichung der Gruppenmittelwerte von dem vorhergesagten Wert ab: Im Vergleich zum ersten Plot sind die Firmen-Punkte näher an der Linie, Teile der Level-2-Residuen konnten also durch den Level-2 Prädiktor erklärt werden.

Daraus hat das Modell 4 auch seinen Namen: Der Prädiktor auf Level-2 (exp_groupmean) sagt den (gruppenspezifischen) Intercept vorher (“Intercept-as-Outcome”).

Code

# Modell fitten
intercept.only.model <- lmer(salary ~ 1 + (1 | company), data = df, REML = TRUE)
# Daten vorbereiten
plot_df_intonly <- df |>
  mutate(
    intercept_only_predictions = predict(intercept.only.model), # level 1 predictions
    intercept_only_intercept = fixef(intercept.only.model)
  ) |> # overall
  group_by(company) |>
  summarise(
    overall = unique(intercept_only_intercept), # fixef intercept
    prediction = unique(intercept_only_predictions), # fixef intercept + ranef intercept
    observed = mean(salary), # actual mean of salary per company
    overall_to_prediction = prediction - overall, # ranef intercept
    company = unique(company), # name of company
    exp_groupmean = mean(experience)
  )
# Plotten
intercept_only_plot <- ggplot(aes(x = exp_groupmean, y = prediction),
  data = plot_df_intonly
) +
  geom_segment(aes(xend = exp_groupmean, yend = overall), alpha = .3) +
  # Color adjustments made here...
  geom_point(
    aes(
      # color = abs(overall_to_prediction),
      size = abs(overall_to_prediction)
    ),
    color = "#33BBEE"
  ) + # Size mapped to abs(residuals)
  geom_point(aes(y = observed)) +
  # scale_color_continuous(low = "black", high = "red") + # Colors to use here
  guides(color = "none", size = "none") + # Color legend removed
  geom_point(aes(y = overall), shape = 1, alpha = .3) +
  ggthemes::theme_tufte() +
  xlab("Average company experience") +
  ylab("Average salary per month") +
  geom_line(aes(x = exp_groupmean, y = overall, group = 1), alpha = 1) +
  ggtitle("Level-2-Residuen des Intercepts (Intercept-only Modell)")

# use ggplotly only if html:
if (knitr::is_html_output()) {
  intercept_only_plot |> plotly::ggplotly(tooltip = "company")
} else {
  intercept_only_plot +
    geom_text(aes(label = company), # add text to plot
      nudge_x = .3
    ) +
    lims(x = c(4, 7)) # zoom out a little, so everything can be seen
}

\(~\)

In dieser Grafik ist der Firmenmittelwert der Erfahrung auf der X-Achse abgebildet, obwohl diese Variable nicht im Modell ist. Dadurch lassen sich die Firmen und deren Residuen bezüglich der Löhne zwischen diesem und der nächsten Grafik direkt vergleichen, weil sich die Datenpunkte in der Darstellung nicht verschieben. Wir erkennen an der horizontalen Gerade der Modellprädiktion entlang der X-Achse, dass der Einfluss der durchschnittlichen Erfahrung nicht mit-modelliert wurde.

Code

context.model <- lmer(salary ~ exp_centered + exp_groupmean + (1 | company), data = df, REML = TRUE)

companies <- df |>
  select(company) |>
  pull() |>
  levels()


plot_df_context <- tibble(
  company = companies,
  exp_groupmean = df |> group_by(company) |> pull(exp_groupmean) |> unique(),
  overall_nogroupmean = fixef(context.model)["(Intercept)"] + mean(exp_groupmean) * fixef(context.model)["exp_groupmean"],
  overall = fixef(context.model)["(Intercept)"] + exp_groupmean * fixef(context.model)["exp_groupmean"],
  prediction = overall + ranef(context.model)$company |> select("(Intercept)") |> pull(),
  overall_to_prediction = ranef(context.model)$company |> select("(Intercept)") |> pull(),
  observed = df |> group_by(company) |> summarise(observed = mean(salary)) |> pull(observed)
)

context_plot <- ggplot(aes(x = exp_groupmean, y = prediction, label = company),
  data = plot_df_context
) +
  geom_segment(aes(xend = exp_groupmean, yend = overall), alpha = .3) +
  # Color adjustments made here...
  geom_point(
    aes(
      # color = abs(overall_to_prediction),
      size = abs(overall_to_prediction)
    ),
    color = "#33BBEE"
  ) + # Size mapped to abs(residuals)
  geom_point(aes(y = observed)) +
  # scale_color_continuous(low = "black", high = "red") + # Colors to use here
  guides(color = FALSE, size = FALSE) + # Color legend removed
  geom_line(aes(y = overall), alpha = .3) +
  ggthemes::theme_tufte() +
  xlab("Average company experience") +
  ylab("Average salary per month") +
  ggtitle("Level-2-Residuen des Intercepts (Kontexteffekt-Modell)")

# use ggplotly only if html:
if (knitr::is_html_output()) {
  context_plot |> plotly::ggplotly(tooltip = "company")
} else {
  context_plot +
    geom_text(aes(label = company), # add text to plot
      nudge_x = .3
    ) +
    lims(x = c(4, 7)) # zoom out a little, so everything can be seen
}

2.2.2 Allgemeines Intercept-as-Outcome-Modell

Dies ist ein Modell mit Random Intercept, mit oder ohne Random Slope und mit Level-2-Haupteffekt (Prädiktion des Intercepts), aber ohne Prädiktor für den Slope (d.h. ohne Cross-Level-Interaktion).

Im Beispiel überprüfen wir, ob der Sektor (sector), in dem die Firmen tätig sind, einen Effekt auf salary hat, d.h. inwiefern die Höhe des Lohns davon abhängt, ob eine Firma im privaten oder öffentlichen Sektor operiert.

Der Effekt eines Level-2-Prädiktors könnte auch ohne weitere (Level-1-)Prädiktoren im Modell ermittelt werden. Wie in unserem Fall ist man häufig aber gleichzeitig am Effekt eines Level-1-Prädiktors (experience) interessiert. Der Level-1-Prädiktor kann in diesem Modell entweder mit (\(\beta_{1i}=\gamma_{10} + \upsilon_{1i}\)) oder ohne Random-Slope (\(\beta_{1i}=\gamma_{10}\)) modelliert werden (der Unterschied bezüglich dieses Effekts ähnelt dem Unterschied zwischen dem Random-Coefficients-Modell und dem Random-Intercept-Modell, vgl. Kapitel 1 Modelle mit Level-1-Prädiktoren). Wir beschränken uns hier auf ein Intercept-as-Outcome-Modell mit Random-Slope.

\(~\)

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

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

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

\(~\)

Berechnen des Modells

Wie bei den anderen Modellen benutzen wir die lmer()-Funktion aus dem Package lme4:

Im Gegensatz zur Variable exp_groupmean (Level-2-Prädiktor im Kontexteffekt-Modell) ist sector eine kategoriale Variable (Faktor). Praktischerweise übernimmt lme4 die Dummy-Kodierung (mit der ersten Faktorstufe Private als Referenzkategorie) für uns.

intercept.as.outcome.model <- lmer(
  salary ~ experience + sector +
    (experience | company),
  data = df, REML = TRUE
)


# Prädiktionen abspeichern als Variable im Datenframe
df$intercept.as.outcome.preds <- predict(intercept.as.outcome.model)

# Modelloutput anschauen
summary(intercept.as.outcome.model)

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

REML criterion at convergence: 1854.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.8397 -0.6713  0.0018  0.6041  3.3846 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 company  (Intercept) 0.48711  0.6979        
          experience  0.01861  0.1364   -0.31
 Residual             1.13631  1.0660        
Number of obs: 600, groups:  company, 20

Fixed effects:
             Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)   5.65787    0.28568 24.12240  19.805  < 2e-16 ***
experience    0.53242    0.04077 19.06039  13.060 5.86e-11 ***
sectorPublic  0.52951    0.34342 18.15358   1.542     0.14    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) exprnc
experience  -0.491       
sectorPublc -0.662  0.062

In der Stichprobe verdienen Angestellte im öffentlichen Sektor im Schnitt ca. 530 CHF mehr als Angestellte im privaten Sektor (\(\hat\gamma_{01} = 0.530\)), dieser Effekt ist aber nicht signifikant (\(p = 0.1404\)) und darf daher nicht weiter interpretiert werden. Auch die Interpretation des Intercepts \(\hat\gamma_{00} = 5.658\) hat sich verändert: Dieser bedeutet jetzt, dass der durchschnittliche Monatslohn für einen Angestellten einer (unbekannten) Firma im privaten Sektor (Referenzkategorie von sector) ohne Erfahrung (an der Stelle experience = 0) ca. 5658 CHF beträgt.

Vertiefung

Wir können das mit der Funktion predict() überprüfen. Dafür müssen wir einen neuen Datenframe erstellen, der alle Modell-relevanten Variablen enthält (experience, sector & company).

# Neues Daten-Frame für die Prädiktion einer Person aus einer
# unbekannten Firma mit 0 Jahre experience und aus dem
# öffentlichen Sektor (die Referenzkategorien).
df.new <- data.frame(
  company = NA,
  sector = "Private",
  experience = 0
)

prediction <- predict(intercept.as.outcome.model, # Modell für die Prädiktion
  newdata = df.new, #  neue Daten, für die salary vorhergesagt werden soll
  re.form = NA
) |> # damit firma = NA sein kann
  print()

       1 
5.657873 

# Test, dass diese Vorhersage == dem Intercept aus dem
# intercept.as.outcome.model ist. Wir können uns hier dafür die Differenz anschauen.
prediction[[1]] - fixef(intercept.as.outcome.model)[["(Intercept)"]]

[1] 0

Die Differenz ist 0, die zwei Werte sind also identisch.

2.3 Slope-as-Outcome-Modell (Modell 5)

In diesem Modell wollen wir untersuchen, ob die in Modell 3 (Random-Coefficents-Modell mit nur Level-1-Prädiktor experience) gefundene (signifikante) Level-2-Varianz des Slopes von experience (\(\hat\sigma_{\upsilon_1}^2\)) durch den Level-2-Prädiktor sector vorhergesagt werden kann, d.h. ob der Effekt von experience vom Sektor, in dem eine Firma operiert, beeinflusst wird.

Auch in Modell 4 haben wir eine ähnlich hohe Slope-Varianz erhalten, diese aber nicht auf Signifikanz getestet. Streng genommen ist Modell 4 jetzt das Vergleichsmodell (und nicht Modell 3), da dieses (wie Modell 5) den Level-2-Effekt von sector auf den Intercept enthält. Dieser kleine Unterschied spielt aber hauptsächlich für die Berechnung der erklärten Varianz des Slopes eine Rolle (s.u.).

\(~\)

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

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

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

\(~\)

2.3.1 Berechnen des Modells

Mit lmer()-Funktion (Package lme4 bzw. lmerTest):

# Berechnen und Abspeichern des Modells
slope.as.outcome.model <- lmer(
  salary ~ experience * sector +
    (experience | company),
  data = df, REML = TRUE
)

# Modell-output anschauen
summary(slope.as.outcome.model)

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

REML criterion at convergence: 1849.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.8271 -0.6923 -0.0319  0.6088  3.4093 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 company  (Intercept) 0.341721 0.58457      
          experience  0.007713 0.08782  0.13
 Residual             1.137612 1.06659      
Number of obs: 600, groups:  company, 20

Fixed effects:
                        Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)              5.28141    0.29513 21.66283  17.895 1.84e-14 ***
experience               0.64430    0.04789 19.20313  13.455 3.15e-11 ***
sectorPublic             1.21269    0.39405 17.30455   3.078  0.00672 ** 
experience:sectorPublic -0.21695    0.06666 18.05226  -3.255  0.00439 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) exprnc sctrPb
experience  -0.566              
sectorPublc -0.749  0.424       
exprnc:sctP  0.406 -0.718 -0.525

Das Ergebnis zeigt eine signifikante Cross-Level-Interaktion (\(\hat\gamma_{11} = -0.217, p = 0.0044\)) in die erwartete Richtung (stärkere Gehaltssteigerungen im privaten Sektor). Der simple slope (bedingte Effekt) von experience zeigt ausserdem die durchschnittliche Gehaltssteigerung im privaten Sektor an (\(\hat\gamma_{10} = 0.644, p = 0)\). Die durchschnittliche Gehaltssteigerung im öffentlichen Sektor können wir als \(\hat\gamma_{10} + \hat\gamma_{11} = 0.644 + (-0.217) = 0.427\) berechnen.

Im Output ausserdem ablesbar (gewissermassen als “Nebeneffekt”) ist das Ergebnis, dass unerfahrene Angestellte (experience = 0) im öffentlichen Sektor signifikant mehr (\(\hat\gamma_{01} = 1.213, p = 0.0067\)) verdienen als im privaten Sektor (bedingter Effekt/simple slope von sector im Term sectorPublic).

2.3.2 Modellvergleich (Signifikanztest) für die verbliebene Slope-Varianz

Jetzt möchten wir noch testen, ob die in Modell 5 verbliebene Slope-Varianz noch signifikant ist. Das machen wir über einen Modellvergleich mit einem Modell, das sich nur darin von unserem Slope-as-Outcome-Modell unterscheidet, dass es keine Slope-Varianz schätzt:

ranova(slope.as.outcome.model)

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

Model:
salary ~ experience + sector + (experience | company) + experience:sector
                                     npar  logLik    AIC    LRT Df Pr(>Chisq)  
<none>                                  8 -924.53 1865.1                       
experience in (experience | company)    6 -927.04 1866.1 5.0152  2    0.08146 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Nicht mehr signifikant! Das gilt (knapp) auch nach Korrektur des p-Wertes:

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

[1] 0.05329462

Nach Berücksichtigung von sector als Prädiktor des Level-1-Effekts von experience (und damit der Cross-Level-Interaktion) gibt es also keine signifikante Slope-Varianz von experience mehr.

2.3.3 Varianzerklärung durch Cross-Level-Interaktionseffekt

Wir wissen nun bereits, dass die in den Modellen 3 und 4 vorhandenen experience-Steigungsunterschiede zwischen den Firmen auf den Sektor, in dem die Firmen operieren, zurückzuführen sind. D.h. dass nach Berücksichtigung der Tatsache, dass sich Private und Public Firmen im Schnitt bezüglich des Einflusses von Erfahrung auf den Lohn unterscheiden, sich innerhalb der beiden Sektoren keine signifikante Varianz des experience-Effekts mehr zwischen den Firmen zeigt.

Rein deskriptiv bleibt aber noch Slope-Varianz übrig. Welcher Anteil der Level-2-Varianz des experience-Slopes wurde also durch sector erklärt (siehe Eid, Gollwitzer, und Schmitt 2017, pp. 753)? Dazu vergleichen wir mit dem Intercept-as-Outcome-Modell von oben, welches sich nur durch den Interaktionseffekt von Modell 5 unterscheidet.

\(R^2_{Cross-level} = \frac{\sigma^2_{\upsilon_14} - \sigma^2_{\upsilon_15}}{\sigma^2_{\upsilon_14}} = \frac{{0.0186} - {0.0077}}{0.0186} = {0.5856}\)

Im Sinne einer Effektgrösse für die Cross-Level-Interaktion (\(\hat\gamma_{11}\)) wurden 58.56 % der Slope-Varianz von experience durch Sektorenunterschiede im Effekt von experience auf salary erklärt.

In den folgenden beiden Plots wird die Auswirkung der Cross-Level-Interaktion auf die Slope-Residuen von experience veranschaulicht. Der erste Plot zeigt die Slope-Residuen relativ zum durchschnittlichen Effekt von experience im Intercept-as-Outcome Modell (\(\hat\gamma_{10} = 0.532\)). Der zweite Plot zeigt die Slope-Residuen im Modell mit Cross-Level-Interaktion (Slope-as-Outcome-Modell) relativ zu den jeweiligen Simple Slopes von experience in den Sektoren Private und Public. Die kleinen schwarzen Punkte stehen wieder für die pro Gruppe in getrennten Regressionen (d.h. ohne Shrinkage) geschätzten Slopes. Klar erkennbar ist, dass die Residuen über alle Firmen hinweg im Slope-as-Outcome Modell deutlich kleiner ausfallen.

Code

intercept.as.outcome <- lmer(salary ~ exp_centered + sector + (-1 + exp_centered | company) + (1 | company), data = df, REML = TRUE)
slope.as.outcome.lm <- lme4::lmList(formula = salary ~ exp_centered | company, data = df)


plot_df_iao_s <- tibble(
  company = companies,
  sector = df |> group_by(company) |> summarise(sector = unique(sector)) |> pull(),
  overall_nosector = fixef(intercept.as.outcome)["exp_centered"],
  prediction = overall_nosector + {
    ranef(intercept.as.outcome)$company |>
      select("exp_centered") |>
      pull()
  },
  overall_to_prediction = {
    ranef(intercept.as.outcome)$company |>
      select("exp_centered") |>
      pull()
  },
  observed = coef(slope.as.outcome.lm) |> select("exp_centered") |> pull()
) |>
  arrange(sector)
plot_df_iao_s <- plot_df_iao_s |>
  mutate(company = factor(company, c(unique(company))))


ggplot(aes(x = company, y = prediction),
  data = plot_df_iao_s
) +
  geom_segment(aes(xend = company, yend = overall_nosector), alpha = .3) +
  # Color adjustments made here...
  geom_point(aes(size = abs(overall_to_prediction), color = sector)) + # size mapped to abs(residuals)
  geom_point(aes(y = observed)) +
  # scale_color_continuous(low = "black", high = "red") + # Colors to use here
  guides(size = "none") + # Color legend removed
  labs(color = "Sector") + # Change title of color legend
  geom_hline(aes(yintercept = plot_df_iao_s$overall_nosector[[1]])) +
  ggthemes::theme_tufte() +
  scale_color_manual(values = c("#0077BB", "#33BBEE")) +
  xlab("") +
  ylab("Effect (slope) of experience on salary") +
  ggtitle("Level-2 Slope-Residuen (Intercept-as-Outcome Modell)") +
  #  geom_text(aes(label = company), hjust = 1, vjust = 0) +# still have to fix it
  theme(
    axis.text.x = element_text(angle = 35, hjust = 1),
    legend.position = c(.9, .9)
  )

Code

slope.as.outcome <- lmer(salary ~ exp_centered * sector + (exp_centered - 1 | company) + (1 | company), data = df, REML = TRUE)


# to get the "observed" slopes without shrinkage we have to fit a linear model for each company:
slope.as.outcome.lm <- lme4::lmList(formula = salary ~ exp_centered | company, data = df)

plot_df_sao <- tibble(
  company = companies,
  sector = df |> group_by(company) |> summarise(sector = unique(sector)) |> pull(),
  overall_nosector = fixef(slope.as.outcome)["exp_centered"] +
    mean(dummy(sector)) * fixef(slope.as.outcome)["exp_centered:sectorPublic"],
  overall = fixef(slope.as.outcome)["exp_centered"] +
    dummy(sector)[, 1] * fixef(slope.as.outcome)[["exp_centered:sectorPublic"]],
  prediction = overall + {
    ranef(slope.as.outcome)$company |>
      select("exp_centered") |>
      pull()
  },
  overall_to_prediction = {
    ranef(slope.as.outcome)$company |>
      select("exp_centered") |>
      pull()
  },
  overall_nosector_to_prediction = prediction - overall_nosector,
  observed = coef(slope.as.outcome.lm) |> select("exp_centered") |> pull()
) |>
  arrange(sector)




# reorder companies
plot_df_sao <- plot_df_sao |>
  mutate(company = factor(company, c(unique(company))))


ggplot(aes(x = company, y = prediction),
  data = plot_df_sao
) +
  geom_segment(aes(xend = company, yend = overall), alpha = .3) +
  geom_point(aes(color = sector, size = abs(overall_to_prediction))) + # Color mapped to abs(residuals)
  geom_point(aes(y = observed)) +
  scale_color_manual(values = c("#0077BB", "#33BBEE")) +
  guides(size = FALSE) + # Color legend removed
  labs(color = "Sector") + # Change title of color legend
  geom_smooth(aes(x = company, y = overall, color = sector, group = sector), method = "lm", se = FALSE, alpha = 1) +
  ggthemes::theme_tufte() +
  xlab("") +
  ylab("Effect (slope) of experience on salary") +
  ggtitle("Level-2 Slope-Residuen (Slope-as-Outcome Modell)") +
  theme(
    axis.text.x = element_text(angle = 35, hjust = 1),
    legend.position = c(.9, .9)
  )

2.4 Übung

Wir verwenden wieder die selben Daten wie in Übung 1.

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

name : Name der Schülerin / desSchülers

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

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

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.

Um diese herunterzuladen und aufzubereiten können wir also auf den Codechunk von letzter Woche zurückgreifen:

# Einlesen
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"))
# Aufbereiten
school_df <- school_df |> mutate(
  id = as.factor(id),
  class = as.factor(class)
)

Die Daten kann man sich etwa so vorstellen (subsample, nur 15 Klassen für die Visualisierung):

Code

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) +
  theme_tufte()

\(~\)

1. Kontexteffekt-Modell

1. Fitten Sie ein Kontexteffekt-Modell, welches die Klassenmittelwerte von grade und motivation als Prädiktoren im Modell hat (zusätzlich zu den Level-1-Effekten dieser Variablen).

Weil wir ja schon, dass es signifikante Random Slope Varianzen gibt, sollten diese auch Teil des Modells sein.

Beachten Sie, dass die Level-1-Prädiktoren der Motivation (motivation) und der Vornote (grade) in diesem Fall als gruppenzentrierte Variablen ins Modell genommen werden müssen. Diese müssen - wie auch die Gruppenmittelwerte dieser Variablen - zunächst berechnet werden.

Lösung

school_df <- school_df |>
  group_by(class) |>
  mutate(
    motivation_groupmean = mean(motivation),
    motivation_centered = motivation - motivation_groupmean,
    grade_groupmean = mean(grade),
    grade_centered = grade - grade_groupmean
  )

context_model <- lmer(
  formula = pisa ~ motivation_groupmean + motivation_centered + grade_groupmean + grade_centered +
    (motivation_centered + grade_centered | class),
  data = school_df
)

2. Interpretieren Sie alle Fixed Effects dieses Modells.

Lösung

summary(context_model)

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

REML criterion at convergence: 13579.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2030 -0.6152  0.0243  0.6313  2.7130 

Random effects:
 Groups   Name                Variance Std.Dev. Corr     
 class    (Intercept)         12951.6  113.80            
          motivation_centered   519.1   22.78   0.50     
          grade_centered        272.7   16.51   0.57 0.23
 Residual                      6382.2   79.89            
Number of obs: 1150, groups:  class, 50

Fixed effects:
                     Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)           427.267    305.723  48.972   1.398   0.1685    
motivation_groupmean   -9.419     58.681  48.446  -0.161   0.8731    
motivation_centered     7.012      3.944  48.697   1.778   0.0816 .  
grade_groupmean        39.239     46.639  49.152   0.841   0.4042    
grade_centered         14.957      2.927  47.024   5.110 5.78e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) mtvtn_g mtvtn_c grd_gr
mtvtn_grpmn -0.880                       
mtvtn_cntrd  0.049 -0.032                
grade_grpmn -0.537  0.076  -0.001        
grade_cntrd  0.013  0.010   0.162   0.005

Insgesamt wurde nur einer der festen Effekte signifikant, nämlich grade_centered. motivation_groupmean, motivation_centered & grade_groupmean scheinen in diesem Modell keinen Einfluss auf den pisa- Wert zu haben.

Auf die Vornote als Prädiktor bezogen bedeutet das: Höhere durchschnittliche Vornoten (grade_groupmean) der Schüler:innen einer Klasse sind nicht mit höheren (durchschnittlichen) pisa-Werten einer Klasse assoziiert. Die Vornote eines Schülers im Vergleich zu Mitschülern aus derselben Klasse (grade_centered) hat dagegen sehr wohl einen Einfluss auf die pisa Werte. Diesen Effekt kennen wir aber bereits aus dem Random-Coefficents-Modell aus der Übung zu Kapitel 1. So haben Schüler:innen derselben Klasse durchschnittlich einen um 14.957 Punkte höheren pisa Wert für jede (ganze) Note, die Sie besser sind als ihre Mitschüler:innen.

Die Motivation scheint hingegen weder auf Klassenebene noch auf Individualebene einen Einfluss auf die pisa-Werte zu haben.

Zu guter Letzt: Der Intercept wäre hier der vorhergesagte pisa-Wert einer Person mit Klassen-durchschnittlicher Motivation und Vornote (motivation_centered = 0 und grade_centered = 0) in einer Klasse, deren Mittelwert sowohl von Motivation und Vornote gleich 0 ist (motivation_groupmean = 0 und grade_groupmean = 0). Da letzteres insbesondere in Bezug auf grade_groupmean unmöglich ist (Notenskala von 1-6), ist eine Interpretation des Intercepts hier nicht wirklich sinnvoll.

2. Intercept-as-Outcome-Modell

1. Welche der Variablen im Datenframe bietet sich an, als Level-2-Prädiktor (auf Klassenebene) eingesetzt zu werden?

Lösung

Die Variable teacher_salary ist auf Klassen-Ebene (also Level-2), weil dieser Wert pro Klasse immer konstant bleibt und verschiedene Klassen unterschiedliche Werte haben. Wir können mit dieser Variable überprüfen, ob der Lohn eines Lehrers/einer Lehrerin mit dem Pisa-Resultat der Schüler zusammenhängt. Auch hier muss man mit der Frage der Kausalität aufpassen: Eine ziemlich offensichtliche konfundierende Variable könnte die Erfahrung des Lehrers/der Lehrerin sein (da erfahrenere Lehrer:innen oft mehr verdienen als weniger erfahrene).

2. Berechnen Sie das Intercept-as-Outcome-Modell mit random Slopes. Tipp: Dieses Modell hat nun 3 Prädiktoren im Modell.

Lösung

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

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

REML criterion at convergence: 13592

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.2012 -0.6111  0.0230  0.6365  2.7374 

Random effects:
 Groups   Name        Variance Std.Dev. Corr       
 class    (Intercept) 10129.1  100.64              
          motivation    514.9   22.69   -0.57      
          grade         274.8   16.58   -0.11  0.23
 Residual              6382.4   79.89              
Number of obs: 1150, groups:  class, 50

Fixed effects:
               Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)     367.899     97.584  50.205   3.770 0.000431 ***
grade            14.965      2.933  47.163   5.102 5.92e-06 ***
motivation        6.967      3.927  48.742   1.774 0.082270 .  
teacher_salary    9.314     14.283  48.816   0.652 0.517385    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) grade  motvtn
grade       -0.046              
motivation  -0.124  0.166       
teachr_slry -0.982 -0.002 -0.006
optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.00624123 (tol = 0.002, component 1)

Der Effekt des Level-2-Prädiktors teacher_salary ist positiv, aber nicht signifikant. Zwar würde der vorhergesagte Pisawert rein deskriptiv für die vorliegende Stichprobe mit zusätzlichen 1000 Franken Lehrergehalt um \(\hat\gamma_{01} =\) 9.314 ansteigen, dieser Effekt darf aber aufgrund der fehlenden Signifikanz nicht weiter interpretiert werden.

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