Szenario & Aufgaben:

Pädagogische Psycholog_innen wollen die Wirksamkeit eines neuen Lernkompetenz-Trainings in einer breit angelegten Studie testen. 500 Schüler_innen aus 25 verschiedenen Schulen (20 pro Schule) werden zufällig der Trainingsgruppe oder einer Kontrollgruppe zugewiesen. Die Schüler_innen in der Trainingsgruppe lernen über ein Schulhalbjahr hinweg Strategien zum selbstständigen und strukturierten Nachbereiten des Unterrichtsstoffs. Die Kontrollgruppe erhält kein Training. Nach Ablauf des Zeitraums werden alle Schüler_innen einem standardisierten Test unterzogen, der das im letzten Schulhalbjahr vermittelte Wissen abprüft.

Die Forscher_innen wollen wissen, ob das Training zu einer besseren Leistung im Test führt. Sie gehen jedoch auch davon aus, dass es Leistungsunterschiede zwischen Schulen geben könnte, und dass das Training möglicherweise nicht an allen Schulen gleich stark wirkt. Sie entscheiden sich, eine Multilevel-Analyse durchzuführen.

Daten einlesen

data <- read.csv("school_training.csv",sep=";")

Zusätzliche Faktorvariablen für school- und training-Variablen erstellen

Diese sind evtl. für Grafiken nötig (je nach Grafik). Inhaltlich gesehen sind beide Variablen Faktoren, die Berechnungen funktionieren aber auch mit den numerischen Versionen.

data$school_id_faktor <- as.factor(data$school_id)
data$training_faktor <- as.factor(data$training)

Pakete laden:

library(Rcmdr)
library(nlme) #für Berechnen des Modells
library(lmtest) #für likelihood ratio tests
library(rcompanion) #für Effektgrößen

Welche Komponenten könnte das Modell beinhalten?

Statistische Hypothesen:

Siehe Formelsammlung:

Grafik

Streudiagramm mit Linien für Schulen:

(zur besseren Ansicht kann die Abbildung in einem separaten Fenster betrachtet werden, dafür im ungeknitteten Markdown rechts unter dem ausgeführten chunk das linke von den drei Symbolen anklicken.)

scatterplot(score~training | school_id_faktor, regLine=TRUE, smooth=FALSE, 
  boxplots=FALSE, by.groups=TRUE, data=data)
## Warning in scatterplot.default(X[, 2], X[, 1], groups = X[, 3], xlab = xlab, : number of groups exceeds number of available colors
##   colors are recycled

Eine andere Möglichkeit ist ein Plot über Mittelwerte, nach den Gruppen Training und Schule aufgeteilt.

Schritt 1: Nullmodell gegen Modell mit random intercept für Schule

Null <- gls(score ~ 1, data = data, method="ML")

summary(Null)
## Generalized least squares fit by maximum likelihood
##   Model: score ~ 1 
##   Data: data 
##        AIC      BIC    logLik
##   4438.412 4446.841 -2217.206
## 
## Coefficients:
##                Value Std.Error  t-value p-value
## (Intercept) 122.3538 0.9131733 133.9875       0
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.93419014 -0.69630750  0.05864086  0.69103267  2.11023754 
## 
## Residual standard error: 20.39875 
## Degrees of freedom: 500 total; 499 residual
RI <- lme(score ~ 1, data= data, random = ~1|school_id, method="ML")

summary(RI)
## Linear mixed-effects model fit by maximum likelihood
##   Data: data 
##       AIC      BIC    logLik
##   4044.41 4057.054 -2019.205
## 
## Random effects:
##  Formula: ~1 | school_id
##         (Intercept) Residual
## StdDev:    16.06302 12.57331
## 
## Fixed effects:  score ~ 1 
##                Value Std.Error  DF  t-value p-value
## (Intercept) 122.3538  3.264708 475 37.47772       0
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.78455259 -0.72031748 -0.04285167  0.72870697  3.09735166 
## 
## Number of Observations: 500
## Number of Groups: 25
lrtest(Null,RI)
## Warning in modelUpdate(objects[[i - 1]], objects[[i]]): original model was of
## class "gls", updated model is of class "lme"
## Likelihood ratio test
## 
## Model 1: score ~ 1
## Model 2: score ~ 1
##   #Df  LogLik Df Chisq Pr(>Chisq)    
## 1   2 -2217.2                        
## 2   3 -2019.2  1   396  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modell mit random intercept wird beibehalten.

Schritt 2: Modell mit zusätzlichem Prädiktor training (fixed effect) gegen Modell mit nur random intercept

RI_pred <- lme(score ~ training, data = data, random = ~1|school_id, method="ML")

summary(RI_pred)
## Linear mixed-effects model fit by maximum likelihood
##   Data: data 
##        AIC      BIC    logLik
##   3868.186 3885.045 -1930.093
## 
## Random effects:
##  Formula: ~1 | school_id
##         (Intercept) Residual
## StdDev:    16.13981 10.42256
## 
## Fixed effects:  score ~ training 
##                Value Std.Error  DF  t-value p-value
## (Intercept) 115.4992  3.301190 474 34.98715       0
## training     13.7092  0.934092 474 14.67650       0
##  Correlation: 
##          (Intr)
## training -0.141
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -2.678646687 -0.684574721  0.005923703  0.634436226  3.081870611 
## 
## Number of Observations: 500
## Number of Groups: 25
lrtest(RI,RI_pred)
## Likelihood ratio test
## 
## Model 1: score ~ 1
## Model 2: score ~ training
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1   3 -2019.2                         
## 2   4 -1930.1  1 178.22  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modell mit RI und fixed effect Prädiktor wird beibehalten.

Modell mit zusätzlich random slope für school gegen Modell mit nur RI für school und fixed effect für Training

RI_pred_RS <- lme(score ~ training, data = data, random = ~training|school_id, method="ML")

summary(RI_pred_RS)
## Linear mixed-effects model fit by maximum likelihood
##   Data: data 
##       AIC      BIC    logLik
##   3857.95 3883.237 -1922.975
## 
## Random effects:
##  Formula: ~training | school_id
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 17.915050 (Intr)
## training     3.816212 -0.935
## Residual    10.237053       
## 
## Fixed effects:  score ~ training 
##                Value Std.Error  DF  t-value p-value
## (Intercept) 115.4992  3.648341 474 31.65801       0
## training     13.7092  1.194413 474 11.47777       0
##  Correlation: 
##          (Intr)
## training -0.686
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.47226959 -0.66105966  0.04941544  0.64941923  3.10773367 
## 
## Number of Observations: 500
## Number of Groups: 25
lrtest(RI_pred,RI_pred_RS)
## Likelihood ratio test
## 
## Model 1: score ~ training
## Model 2: score ~ training
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1   4 -1930.1                         
## 2   6 -1923.0  2 14.237  0.0008101 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(RI_pred,RI_pred_RS)
##            Model df      AIC      BIC    logLik   Test L.Ratio p-value
## RI_pred        1  4 3868.186 3885.045 -1930.093                       
## RI_pred_RS     2  6 3857.950 3883.237 -1922.975 1 vs 2 14.2367   8e-04

Das Modell mit random intercept und random slope für school und fixed effect für Training beschreibt die Daten am besten.

Effektstärke

Achtung: Die Reihenfolge ist hier bedeutsam. Der Test bestimmt, wie viel Varianz das erste Modell im Verhältnis zum zweiten Modell erklärt Das erste Modell sollte also Modell mit unseren gewählten Prädiktoren sein (ansonsten bekommen wir heraus, wie viel Varianz das Nullmodell im Verhältnis zum Modell mit Prädiktoren enthält, also meistens etwas Negatives, da das Nullmodell ja normalerweise weniger erklärt).

nagelkerke(RI_pred_RS,Null)
## $Models
##                                                                        
## Model: "lme.formula, score ~ training, data, ~training | school_id, ML"
## Null:  "gls, score ~ 1, data, ML"                                      
## 
## $Pseudo.R.squared.for.model.vs.null
##                              Pseudo.R.squared
## McFadden                             0.132704
## Cox and Snell (ML)                   0.691775
## Nagelkerke (Cragg and Uhler)         0.691872
## 
## $Likelihood.ratio.test
##  Df.diff LogLik.diff  Chisq     p.value
##       -4     -294.23 588.46 4.8659e-126
## 
## $Number.of.observations
##           
## Model: 500
## Null:  500
## 
## $Messages
## [1] "Note: For models fit with REML, these statistics are based on refitting with ML"
## 
## $Warnings
## [1] "None"

Die negative Korrelation zeigt an, dass Schüler_innen an Schulen mit einem niedrigeren Intercept eine stärkere Steigung im Kriterium durch den Trainings-Faktor erfahren als solche an Schulen mit höherem intercept (zwischen den Trainings-Gruppen und den Kontrollgruppen). Schüler_innen an Schulen mit durchschnittlich etwas geringerer Leistung profitieren also vermutlich stärker vom Training als solche von Schulen mit einer bereits höheren durchschnittlichen Leistung. Die Forschungsgruppe könnte daher empfehlen, das Training vorrangig an den Schulen anzubieten, die mehr davon profitieren (und dadurch ggf. im Vergleich zu den anderen Schulen aufholen könnten).

Nichts :)