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.
data <- read.csv("school_training.csv",sep=";")
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
Siehe Formelsammlung:
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.
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.
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.
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.
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 :)