Bei einer mehrstufigen Stichprobenziehung handelt es sich nicht wie sonst üblich um eine einfache Zufallsstichprobe, bei der alle Elemente der Population die gleiche Wahrscheinlichkeit haben, gezogen zu werden. Bei beispielsweise drei Hierarchieebenen einer Untersuchung mit Schülern findet zuerst eine Zufallsauswahl von Schulen (Level-3), dann eine Zufallsauswahl von Schulklassen (Level-2) in diesen Schulen, dann eine Zufallsauswahl von Schülern (Level-1) aus den Schulklassen statt. Dies führt zu einer verschachtelten oder „genesteten“ Datenstruktur, bei der Stichprobenziehung handelt es sich nicht mehr um zufällige Ziehungen aus einer Population, sondern es gibt Abhängigkeiten: der Wert einer Person auf einem Merkmal hängt (auch) davon ab, aus welcher Level-2-Einheit (z.B. Schulklasse) sie stammt.
Beispiel: Schüler einer Klasse sind sich bezüglich der gezeigten Aggressivität ähnlicher als Schüler aus unterschiedlichen Klassen, da es gruppenspezifische Einflüsse wie das Klassenklima gibt, die sich auf das allgemeine Aggressivitätsniveau auswirken. Ähnlichkeiten können sich auch auf die ganze Schule beziehen, z.B. durch eine sozioökonomische Selektion von Schülern oder durch das Führungsklima an einer Schule (vgl. Eid, Gollwitzer & Schmitt, 2013).
Die Verletzung der Unabhängigkeitsannahme führt zu interpretatorischen und statistischen Problemen.
In HLM-Modellen werden die Untersuchungseinheiten einer höheren Ebene (z.B. die Level-2-Einheiten “Schulklassen”) als Zufallsauswahl aus einer Population von Untersuchungseinheiten (z.B. Schulklassen, Kulturen) konzeptualisiert. So interessiert also nicht mehr (wie z.B. in einer ANOVA) der Effekt einer einzelnen Klasse/Kultur (Abweichung vom Gesamtmittelwert), sondern nur die Varianz der Effekte. Dasselbe gilt für die Abweichung eines Regressionsgewichts einer Level-2-Einheit vom Gesamtregressionsgewicht.
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- oder Residualvarianz), desto grösser sind die Ähnlichkeiten innerhalb der Level-2-Einheiten im Vergleich zu zwischen den Level-2-Einheiten, da die Level-1-Varianz - also die Varianz innerhalb der Gruppen (separat berechnet für alle Gruppen und dann gepoolt/gemittelt) - dann relativ klein ist.
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.u.), 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 notwendig.
Bei Modellen OHNE Level-2-Prädiktoren handelt es sich um normale Regressionsmodelle, die die Mehrebenenstruktur der Daten berücksichtigen. Man unterscheidet hier zwischen Modellen, bei denen kein zufälliger Effekt für den/die Slope(s) der Level-1-Prädiktoren spezifiziert wird (d.h. ein Random-Effect wird nur für den Intercept definiert, daher: “Random-Intercept Modelle”) und solchen, bei denen zusätzlich zur Varianz des random Intercepts auch noch ein bzw. mehrere Varianzen für den/die random Slope/s geschätzt werden (“Random-Coefficients Modelle”).
Bei Modellen MIT Level-2-Prädiktoren unterscheidet man zwischen solchen, bei denen nur der Intercept durch ein bzw. mehrere Level-2 Merkmale vorhergesagt wird und solchen, bei denen sowohl die Intercept-Varianz als auch die Slope-Varianz durch einen bzw. mehrere Level-2-Prädiktoren vorhergesagt wird. Erstere bezeichnet man als “Intercept-as-Outcome-Modelle”, letztere als “Slope-as-Outcome-Modelle”.
Bei der ML-Schätzung werden alle Parameter (Fixed und Random) mittels Maximum-Likelihod geschätzt, bei REML nur die Random Effects (die Fixed Effects werden mittels GLM geschätzt). Für Modellvergleiche (Likelihood-Ratio-Tests) von Modellen, die sich nur im Random-Part unterscheiden, genügt REML, sonst benötigt man ML. Da REML die besseren (ungebiasten) Parameterschätzungen liefert, ist dieses zu bevorzugen. Ausserdem werden Modellvergleiche üblicherweise nur für den Signifikanztest von Random Parameters gemacht, da die Fixed Effects mittels T-Tests getestet werden können.
Mayer und Trommsdorff (2012) haben den Zusammenhang zwischen der Religiosität (religiosity) Jugendlicher und der von Jugendlichen angegebenen Wichtigkeit traditioneller Familienwerte (family) in 17 verschiedenen Kulturen (n = 4719) untersucht. Dabei wurde nur die Intensität der persönlichen Religiosität erfragt. Zusätzlich liegen Angaben zum Anteil der Religionszugehörigkeit (prop_rel) der Jugendlichen in den verschiedenen Kulturen vor (Level-2-Merkmal, z.B. Schweiz = 0,74: 74 % der Jugendlichen gehören einer Religion an).
Folgende Hypothesen sollen mit hierarchischen linearen Modellen überprüft werden:
Zusatzfrage: Wie gross ist der Anteil dieser Level-2-Varianz an der Gesamtvarianz der abhängigen Variablen (Intraklassenkorrelation)?
Vorbereitung der Daten: Importieren des Datensatzes relfam3.csv und speichern als .RData-File.
relfam3<-read.csv("relfam3.csv", header = TRUE)
save(relfam3, file="relfam3.RData")
load("relfam3.Rdata")
names(relfam3)## [1] "X" "culture" "religiosity" "family" "prop_rel"tail(relfam3)## X culture religiosity family prop_rel
## 4714 4714 15 1 3.8 0.16
## 4715 4715 15 1 4.0 0.16
## 4716 4716 15 1 4.0 0.16
## 4717 4717 15 1 4.0 0.16
## 4718 4718 15 1 4.0 0.16
## 4719 4719 15 1 4.0 0.16Packages nlme, lme4 und lmerTest installieren:
install.packages(nlme)
install.packages(lme4)
install.packages(lmerTest)
library(nlme)
library(lme4)
library(lmerTest)Dieses Modell enthält keine Prädiktorvariable, nur der Intercept wird als Kombination eines festen (durchschnittlichen) Effekts und einer zufälligen Abweichung (Level-2-Residuum) von diesem modelliert:
Level-1-Modell: \(\boldsymbol {y_{mi}=\beta_{0i}+\epsilon_{mi}}\)
Level-2-Modell: \(\boldsymbol {\beta_{0i}=\gamma_{00} + \upsilon_{0i}}\)
Gesamtmodell: \(\boldsymbol {y_{mi}=\gamma_{00} + \upsilon_{0i}+\epsilon_{mi}}\)
Mit lme-Funktion (Package nlme):
intercept.only <- lme(family ~ 1, random = ~1 | culture, data=relfam3, method = "REML")
summary(intercept.only)## Linear mixed-effects model fit by REML
## Data: relfam3
## AIC BIC logLik
## 7491.616 7510.993 -3742.808
##
## Random effects:
## Formula: ~1 | culture
## (Intercept) Residual
## StdDev: 0.3330947 0.5302261
##
## Fixed effects: family ~ 1
## Value Std.Error DF t-value p-value
## (Intercept) 4.236973 0.08118324 4702 52.19024 0
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -6.7676816 -0.5352317 0.1943606 0.6650780 3.1466657
##
## Number of Observations: 4719
## Number of Groups: 17VarCorr(intercept.only)## culture = pdLogChol(1)
## Variance StdDev
## (Intercept) 0.1109521 0.3330947
## Residual 0.2811397 0.5302261Mit lmer-Funktion (Package lme4 bzw. lmerTest):
intercept.only2 <- lmer(family ~ 1 + (1 | culture), relfam3, REML = TRUE)
summary(intercept.only2)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [merModLmerTest]
## Formula: family ~ 1 + (1 | culture)
## Data: relfam3
##
## REML criterion at convergence: 7485.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.7677 -0.5352 0.1944 0.6651 3.1467
##
## Random effects:
## Groups Name Variance Std.Dev.
## culture (Intercept) 0.1110 0.3331
## Residual 0.2811 0.5302
## Number of obs: 4719, groups: culture, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.23697 0.08118 15.94500 52.19 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Intraklassenkorrelation: \(\boldsymbol {\rho = \frac{\sigma^2_{Level-2}}{\sigma^2_{Level-2}+ \sigma^2_{Level-1}}= \frac {0.1109521}{0.1109521+0.2811397}}\)= 0.2829748
-> ca. 28% der Gesamtvarianz sind auf Level-2-(Kultur-)Unterschiede zurückzuführen.
Wir kennen nun also die Level-2-Varianz, aber wir haben keinen Signifikanztest für den zufälligen Effekt. Dieser kann über einen Modellvergleich (Likelihood-Ratio-Test) durchgeführt werden. Dazu müssen wir zusätzlich ein Modell spezifizieren, dass keinen Random-Effekt enthält, d.h. ein “normales” lineares Modell ohne Prädiktor. Das könnte man mit lm machen, aber wir machen es hier mit gls aus dem nlme-Package, da dieses auch REML beherrscht. Mittels anova oder rand (aus lmerTest) können wir dann die beiden Modelle mit einem LR-Test vergleichen.
rand(intercept.only2)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## culture 1200 1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1no.random2 <- gls(family ~ 1, data=relfam3, method = "REML")
summary(no.random2)## Generalized least squares fit by REML
## Model: family ~ 1
## Data: relfam3
## AIC BIC logLik
## 8689.353 8702.271 -4342.677
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 4.272202 0.00883475 483.568 0
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -5.3916403 -0.4485105 0.2105735 0.8696575 1.1991994
##
## Residual standard error: 0.6069029
## Degrees of freedom: 4719 total; 4718 residualanova(intercept.only, no.random2)## Model df AIC BIC logLik Test L.Ratio
## intercept.only 1 3 7491.616 7510.993 -3742.808
## no.random2 2 2 8689.353 8702.271 -4342.677 1 vs 2 1199.737
## p-value
## intercept.only
## no.random2 <.0001Dieses Modell enthält einen Level-1-Prädiktor (X), der jedoch als fester Effekt konzeptualisiert wird (kein Random Slope, nur Random Intercept). Dieses Modell wird benötigt, um den Anteil der durch den Level-1-Prädiktor erklärten Level-1-Residualvarianz zu ermitteln (vgl. Eid et al., 2013, s.u.).
Level-1-Modell:
\(\boldsymbol {y_{mi}=\beta_{0i} + \beta_{1i} \cdot x_{mi} + \epsilon_{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}+\epsilon_{mi}}\)
Mit lme-Funktion (Package nlme):
random.intercept <- lme(family ~ religiosity, random = ~1 | culture, data=relfam3,
method = "REML")
summary(random.intercept)## Linear mixed-effects model fit by REML
## Data: relfam3
## AIC BIC logLik
## 7271.415 7297.251 -3631.708
##
## Random effects:
## Formula: ~1 | culture
## (Intercept) Residual
## StdDev: 0.271298 0.5178104
##
## Fixed effects: family ~ religiosity
## Value Std.Error DF t-value p-value
## (Intercept) 3.921284 0.06936454 4701 56.53155 0
## religiosity 0.106045 0.00688554 4701 15.40118 0
## Correlation:
## (Intr)
## religiosity -0.296
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -7.1797255 -0.5169863 0.1532851 0.6844567 2.9233746
##
## Number of Observations: 4719
## Number of Groups: 17VarCorr(random.intercept)## culture = pdLogChol(1)
## Variance StdDev
## (Intercept) 0.07360259 0.2712980
## Residual 0.26812762 0.5178104Mit lmer-Funktion (Package lme4 bzw. lmerTest):
random.intercept2 <- lmer(family ~ religiosity + (1 | culture), relfam3, REML = TRUE)
summary(random.intercept2)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [merModLmerTest]
## Formula: family ~ religiosity + (1 | culture)
## Data: relfam3
##
## REML criterion at convergence: 7263.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.1797 -0.5170 0.1533 0.6845 2.9234
##
## Random effects:
## Groups Name Variance Std.Dev.
## culture (Intercept) 0.0736 0.2713
## Residual 0.2681 0.5178
## Number of obs: 4719, groups: culture, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.921e+00 6.936e-02 1.900e+01 56.53 <2e-16 ***
## religiosity 1.060e-01 6.886e-03 4.596e+03 15.40 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## religiosity -0.296Wieviel Level-1-Varianz wurde durch religiosity erklärt?
\(\boldsymbol {R^2_{Level-1} = \frac {\sigma^2_{\epsilon_1} - \sigma^2_{\epsilon_2}}{\sigma^2_{\epsilon_1}}= \frac {0.2811397-0.26812762}{0.2811397}}\)= 0.0462833
Welcher Anteil der Gesamtvarianz (Level-1- & Level-2-Varianz) wurde durch religiosity erklärt?
\(\boldsymbol {R^2_{Gesamt} = \frac {(\sigma^2_{\upsilon_01} + \sigma^2_{\epsilon_1}) - (\sigma^2_{\upsilon_02} + \sigma^2_{\epsilon_2})}{\sigma^2_{\upsilon_01} + \sigma^2_{\epsilon_1}}= \frac {(0.1109521 + 0.2811397) - (0.07360259 + 0.26812762)}{0.1109521 + 0.2811397}}\)= 0.1284434
Dieses Modell enthält sowohl einen Random Intercept als auch einen Random Slope. Mit Hilfe dieses Modells kann die zufällige Varianzkomponente des Level-1-Slopes berechnet, also die Frage beantwortet werden, in welchem Ausmass sich der Effekt des Level-1-Prädiktors (X) zwischen den Level-2-Einheiten unterscheidet.
Level-1-Modell:
\(\boldsymbol {y_{mi}=\beta_{0i} + \beta_{1i} \cdot x_{mi} + \epsilon_{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} + \epsilon_{mi}}\)
Mit lme-Funktion (Package nlme):
random.coefficients <- lme(family ~ religiosity, random = ~religiosity | culture,
data=relfam3, method = "REML")
summary(random.coefficients)## Linear mixed-effects model fit by REML
## Data: relfam3
## AIC BIC logLik
## 7253.803 7292.556 -3620.901
##
## Random effects:
## Formula: ~religiosity | culture
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.34203231 (Intr)
## religiosity 0.04790866 -0.705
## Residual 0.51563602
##
## Fixed effects: family ~ religiosity
## Value Std.Error DF t-value p-value
## (Intercept) 3.902325 0.08672852 4701 44.99471 0
## religiosity 0.106681 0.01376499 4701 7.75016 0
## Correlation:
## (Intr)
## religiosity -0.711
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -7.0766333 -0.5198621 0.1497904 0.6807766 2.8625981
##
## Number of Observations: 4719
## Number of Groups: 17VarCorr(random.coefficients)## culture = pdLogChol(religiosity)
## Variance StdDev Corr
## (Intercept) 0.11698610 0.34203231 (Intr)
## religiosity 0.00229524 0.04790866 -0.705
## Residual 0.26588051 0.51563602Mit lmer-Funktion (Package lme4 bzw. lmerTest):
random.coefficients2 <- lmer(family ~ religiosity + (religiosity | culture), relfam3, REML = TRUE)
summary(random.coefficients2)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [merModLmerTest]
## Formula: family ~ religiosity + (religiosity | culture)
## Data: relfam3
##
## REML criterion at convergence: 7241.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.0766 -0.5199 0.1498 0.6808 2.8626
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## culture (Intercept) 0.116986 0.34203
## religiosity 0.002295 0.04791 -0.70
## Residual 0.265881 0.51564
## Number of obs: 4719, groups: culture, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.90232 0.08673 16.08900 44.99 < 2e-16 ***
## religiosity 0.10668 0.01377 14.92700 7.75 1.31e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## religiosity -0.711Wenn wir dieses Modell mittels LR-Test mit dem Random Intercept Modell vergleichen, werden zwei Parameter gleichzeitig getestet: die Level-2-Varianz des Slopes (\(\boldsymbol {\sigma^2_{\upsilon_1}}\)) und die Level-2-Kovarianz/Korrelation zwischen Intercept und Slope (\(\boldsymbol {\sigma_{\upsilon0\upsilon1}}\)). Im Folgenden der Vergleich mittels rand und anschliessend mit anova (einmal in Bezug auf die mit lme geschätzten Modelle und einmal in in Bezug auf die mit lmer geschätzten Modelle:
rand(random.coefficients2)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## religiosity:culture 21.6 2 2e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(random.coefficients, random.intercept)## Model df AIC BIC logLik Test L.Ratio
## random.coefficients 1 6 7253.803 7292.556 -3620.901
## random.intercept 2 4 7271.415 7297.251 -3631.708 1 vs 2 21.61255
## p-value
## random.coefficients
## random.intercept <.0001anova(random.coefficients2, random.intercept2, refit = FALSE)## Data: relfam3
## Models:
## ..1: family ~ religiosity + (1 | culture)
## object: family ~ religiosity + (religiosity | culture)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 4 7271.4 7297.3 -3631.7 7263.4
## object 6 7253.8 7292.6 -3620.9 7241.8 21.613 2 2.027e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Bei der anova der beiden mit lmer geschätzten Modelle (random.coefficients2 vs. random.intercept2) 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 (s.o. 1.4). 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.
Das Ergebnis zeigt einen hochsignifikanten LR-Test, beide Parameter zusammen unterscheiden sich also signifikant von 0. Der p-Wert ist allerdings nicht ganz korrekt, da einer der beiden Parameter - \(\boldsymbol {\sigma^2_{\upsilon_1}}\) - einen bounded parameter space aufweist (kann nicht kleiner als 0 werden). Der korrekte p-Wert ist daher noch kleiner als der ausgegebene, da die korrekte Testverteilung eine Mischverteilung aus einer \(\boldsymbol {\chi^2_{df=1}}\) und einer \(\boldsymbol {\chi^2_{df=2}}\) ist. Die exakte Bestimmung des p-Werts ist wegen der Mischung von bounded (\(\boldsymbol {\sigma^2_{\upsilon_1}}\)) und unbounded (\(\boldsymbol {\sigma_{\upsilon0\upsilon1}}\)) Parametern schwierig (vgl. West, Welch, & Galecki, 2015).
Dieser Test beantwortet jedoch nicht die hauptsächlich interessierende Frage, ob die Level-2-Varianz des Slopes als solche signifikant ist. Dazu brauchen wir einen zusätzlichen Modellvergleich zwischen einem Modell, das nur die Level-2-Varianz des Slopes (\(\boldsymbol {\sigma^2_{\upsilon_1}}\)) zusätzlich im Vergleich zum Random-Intercept-Modell beinhaltet und daher gleichzeitig die Level-2-Kovarianz/Korrelation zwischen Intercept und Slope (\(\boldsymbol {\sigma_{\upsilon0\upsilon1}}\)) auf 0 restringiert (d.h. ein Modell mit einer Varianzkomponenten-Kovarianzstruktur anstatt eines mit einer unstrukturierten Kovarianzstruktur).
Ein solches Modell lässt sich (hier ausschliesslich mit lmer gezeigt) wie folgt spezifizieren (Bates, 2005):
random.coefficients.varcomp <- lmer(family ~ religiosity + (1 | culture) +
(religiosity - 1 | culture), relfam3, REML = TRUE)
summary(random.coefficients.varcomp)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [merModLmerTest]
## Formula:
## family ~ religiosity + (1 | culture) + (religiosity - 1 | culture)
## Data: relfam3
##
## REML criterion at convergence: 7249
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.1106 -0.5210 0.1457 0.6798 2.9881
##
## Random effects:
## Groups Name Variance Std.Dev.
## culture (Intercept) 0.095759 0.30945
## culture.1 religiosity 0.001766 0.04203
## Residual 0.266006 0.51576
## Number of obs: 4719, groups: culture, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.90378 0.07884 16.54500 49.516 < 2e-16 ***
## religiosity 0.10469 0.01253 14.95700 8.357 5.12e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## religiosity -0.161Den Modellvergleich für den Test der Signifikanz der Level-2-Varianz des Slopes (\(\boldsymbol {\sigma^2_{\upsilon_1}}\)) erhalten wir über die anova-Funktion. Wir können so auch gleich alle drei genesteten Modelle (Random-Intercept, Random-Coefficients mit Varianzkomponeten-Kovarianzstruktur, Random-Coefficients mit unstrukturierter Kovarianzstruktur) miteinander vergleichen, dann erhalten wir auch einen separaten Signifikanztest für die Kovarianz \(\boldsymbol {\sigma_{\upsilon0\upsilon1}}\):
anova(random.coefficients2, random.coefficients.varcomp, random.intercept2, refit = FALSE)## Data: relfam3
## Models:
## ..2: family ~ religiosity + (1 | culture)
## ..1: family ~ religiosity + (1 | culture) + (religiosity - 1 | culture)
## object: family ~ religiosity + (religiosity | culture)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..2 4 7271.4 7297.3 -3631.7 7263.4
## ..1 5 7259.0 7291.2 -3624.5 7249.0 14.463 1 0.0001429 ***
## object 6 7253.8 7292.6 -3620.9 7241.8 7.149 1 0.0075005 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hier verwenden wir wieder refit = FALSE (s.o.), um den LR-Test für die REML-geschätzten Modelle zu erhalten. Beide Vergleiche sind hochsignifikant, also sind beide (Ko-)Varianz-Parameter (\(\boldsymbol {\sigma^2_{\upsilon_1}}\) und \(\boldsymbol {\sigma_{\upsilon0\upsilon1}}\)) auch einzeln signifikant. Der p-Wert des zweiten Vergleichs ist korrekt (p = 0.0075), da es sich bei \(\sigma_{\upsilon0\upsilon1}\) um einen Parameter mit unbounded parameter space handelt. Der p-Wert des ersten Vergleichs (mit bounded parameter space: \(\boldsymbol {\sigma^2_{\upsilon_1}}\)) ist muss allerdings noch halbiert werden, da die Testverteilung (Mischverteilung: \(\boldsymbol {\frac {1}{2} \chi^2_{df=0} + \frac {1}{2}\chi^2_{df=1}}\)) mit dem Test einer gerichteten Alternativhypothese \(\boldsymbol {(H_1:\sigma^2_{\upsilon_1}>0)}\) gleichzusetzen ist (p = 0.00007).
Dies ist ein Modell mit Random Slope und mit Level-2-Haupteffekt (Prädiktion des Intercepts), aber ohne Prädiktor für den Slope (ohne Cross-level-Interaktion). Dieses Modell wird als Vergleichsmodell zu Modell 5 (s.u.) benötigt wird, um die Varianzreduktion durch den Cross-level-Effekt (Prädiktion des Slopes) zu quantifizieren.
Level-1-Modell:
\(\boldsymbol {y_{mi}=\beta_{0i} + \beta_{1i} \cdot x_{mi} + \epsilon_{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} + \epsilon_{mi}}\)
Mit lme-Funktion (Package nlme):
intercept.outcome <- lme(family ~ religiosity + prop_rel, random = ~religiosity | culture,
data=relfam3, method = "REML")
summary(intercept.outcome)## Linear mixed-effects model fit by REML
## Data: relfam3
## AIC BIC logLik
## 7249.254 7294.465 -3617.627
##
## Random effects:
## Formula: ~religiosity | culture
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.40393068 (Intr)
## religiosity 0.05130595 -0.924
## Residual 0.51561079
##
## Fixed effects: family ~ religiosity + prop_rel
## Value Std.Error DF t-value p-value
## (Intercept) 3.457511 0.15829647 4701 21.841995 0.000
## religiosity 0.103698 0.01446386 4701 7.169454 0.000
## prop_rel 0.625810 0.16804673 15 3.724026 0.002
## Correlation:
## (Intr) rlgsty
## religiosity -0.497
## prop_rel -0.769 -0.089
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -7.0965427 -0.5219249 0.1610212 0.6618710 2.8164437
##
## Number of Observations: 4719
## Number of Groups: 17VarCorr(intercept.outcome)## culture = pdLogChol(religiosity)
## Variance StdDev Corr
## (Intercept) 0.1631600 0.40393068 (Intr)
## religiosity 0.0026323 0.05130595 -0.924
## Residual 0.2658545 0.51561079Mit lmer-Funktion (Package lme4 bzw. lmerTest):
intercept.outcome2 <- lmer(family ~ religiosity + prop_rel + (religiosity | culture),
relfam3, REML = T)
summary(intercept.outcome2)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [merModLmerTest]
## Formula: family ~ religiosity + prop_rel + (religiosity | culture)
## Data: relfam3
##
## REML criterion at convergence: 7235.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.0966 -0.5219 0.1610 0.6619 2.8164
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## culture (Intercept) 0.163191 0.4040
## religiosity 0.002631 0.0513 -0.92
## Residual 0.265855 0.5156
## Number of obs: 4719, groups: culture, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.45736 0.15828 23.52800 21.843 < 2e-16 ***
## religiosity 0.10370 0.01446 15.03500 7.170 3.18e-06 ***
## prop_rel 0.62604 0.16800 17.28400 3.726 0.00164 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) rlgsty
## religiosity -0.497
## prop_rel -0.769 -0.089Level-1-Modell:
\(\boldsymbol {y_{mi}=\beta_{0i} + \beta_{1i} \cdot x_{mi} + \epsilon_{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} + \epsilon_{mi}}\)
slope.outcome <- lme(family ~ religiosity + prop_rel + religiosity*prop_rel,
random = ~religiosity | culture, data=relfam3, method = "REML")
summary(slope.outcome)## Linear mixed-effects model fit by REML
## Data: relfam3
## AIC BIC logLik
## 7249.125 7300.793 -3616.562
##
## Random effects:
## Formula: ~religiosity | culture
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.34564038 (Intr)
## religiosity 0.03956303 -0.891
## Residual 0.51557816
##
## Fixed effects: family ~ religiosity + prop_rel + religiosity * prop_rel
## Value Std.Error DF t-value p-value
## (Intercept) 3.930033 0.21940080 4700 17.912574 0.0000
## religiosity 0.021205 0.03187788 4700 0.665209 0.5059
## prop_rel -0.044307 0.28433547 15 -0.155826 0.8782
## religiosity:prop_rel 0.114005 0.04072977 4700 2.799059 0.0051
## Correlation:
## (Intr) rlgsty prp_rl
## religiosity -0.776
## prop_rel -0.917 0.719
## religiosity:prop_rel 0.727 -0.926 -0.802
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -7.1380611 -0.5293355 0.1671173 0.6635869 2.8783206
##
## Number of Observations: 4719
## Number of Groups: 17VarCorr(slope.outcome)## culture = pdLogChol(religiosity)
## Variance StdDev Corr
## (Intercept) 0.119467275 0.34564038 (Intr)
## religiosity 0.001565233 0.03956303 -0.891
## Residual 0.265820840 0.51557816Mit lmer-Funktion (Package lme4 bzw. lmerTest):
slope.outcome2 <- lmer(family ~ religiosity + prop_rel + religiosity*prop_rel +
(religiosity | culture), relfam3, REML = TRUE)
summary(slope.outcome2)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [merModLmerTest]
## Formula:
## family ~ religiosity + prop_rel + religiosity * prop_rel + (religiosity |
## culture)
## Data: relfam3
##
## REML criterion at convergence: 7233.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.1379 -0.5293 0.1668 0.6635 2.8785
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## culture (Intercept) 0.119453 0.34562
## religiosity 0.001567 0.03958 -0.89
## Residual 0.265819 0.51558
## Number of obs: 4719, groups: culture, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.93006 0.21939 13.84700 17.914 5.67e-11 ***
## religiosity 0.02120 0.03189 16.03700 0.665 0.5156
## prop_rel -0.04438 0.28432 14.56700 -0.156 0.8781
## religiosity:prop_rel 0.11402 0.04074 15.63100 2.798 0.0131 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) rlgsty prp_rl
## religiosity -0.776
## prop_rel -0.917 0.718
## rlgsty:prp_ 0.727 -0.926 -0.801Wieviel Level-2-Varianz des religiosity-Slopes wurde durch prop_rel erklärt (vgl. Eid et al., 2013)?
\(\boldsymbol {R^2_{Cross-level} = \frac {\sigma^2_{\upsilon_14} - \sigma^2_{\upsilon_15}}{\sigma^2_{\upsilon_14}}= \frac {0.0026323-0.001565233}{0.0026323}}\)= 0.4053744
Um zu überprüfen, ob die Varianz des Level-2 Slopes noch signifikant ist, nachdem wir prop_rel in das Modell aufgenommen haben, brauchen wir ein Vergleichsmodell zum slope.outcome2-Modell, das sich nur darin von diesem unterscheidet, dass es den Random-Slope Parameter \(\boldsymbol {\sigma^2_{\upsilon_1}}\) nicht enthält. Da die rand-Funktion wieder beide Parameter (\(\boldsymbol {\sigma^2_{\upsilon_1}}\) und \(\boldsymbol {\sigma_{\upsilon0\upsilon1}}\)) auf einmal testet, müssen wir zwei weitere Modelle spezifizieren: eines ohne Random-Slope Parameter (aber mit Cross-level-Interaktion) und eines mit Random-Slope Parameter, aber mit einer Varianzkomponenten-Kovarianzstruktur (d.h. mit \(\boldsymbol {\sigma^2_{\upsilon_1}}\) aber ohne \(\boldsymbol {\sigma_{\upsilon0\upsilon1}}\)), und diese beiden Modelle müssen dann mit unserem slope.outcome2-Modell mittels anova verglichen werden.
Zunächst aber der Gesamtvergleich zum simultanen Test beider Parameter mittels rand:
rand(slope.outcome2)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## religiosity:culture 14.7 2 6e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Beide Parameter zusammen sind nach wie vor hoch signifikant (tatsächlicher p-Wert noch kleiner wegen bounded parameter space bei \(\boldsymbol {\sigma^2_{\upsilon_1}}\)).
slope.outcome.intercept <- lmer(family ~ religiosity + prop_rel + religiosity*prop_rel +
(1 | culture), relfam3, REML = TRUE)
summary(slope.outcome.intercept)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [merModLmerTest]
## Formula: family ~ religiosity + prop_rel + religiosity * prop_rel + (1 |
## culture)
## Data: relfam3
##
## REML criterion at convergence: 7247.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.2591 -0.5220 0.1748 0.6724 3.0639
##
## Random effects:
## Groups Name Variance Std.Dev.
## culture (Intercept) 0.06837 0.2615
## Residual 0.26704 0.5168
## Number of obs: 4719, groups: culture, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.95201 0.16776 16.00000 23.557 5.95e-14 ***
## religiosity 0.01806 0.02000 4714.00000 0.903 0.367
## prop_rel -0.07882 0.21815 17.00000 -0.361 0.722
## religiosity:prop_rel 0.11729 0.02518 4699.00000 4.658 3.28e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) rlgsty prp_rl
## religiosity -0.200
## prop_rel -0.916 0.203
## rlgsty:prp_ 0.204 -0.938 -0.255slope.outcome.varcomp <- lmer(family ~ religiosity + prop_rel + religiosity*prop_rel +
(1 | culture) + (religiosity - 1 | culture), relfam3,
REML = TRUE)
summary(slope.outcome.varcomp)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [merModLmerTest]
## Formula: family ~ religiosity + prop_rel + religiosity * prop_rel + (1 |
## culture) + (religiosity - 1 | culture)
## Data: relfam3
##
## REML criterion at convergence: 7244.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.1789 -0.5247 0.1692 0.6669 3.0563
##
## Random effects:
## Groups Name Variance Std.Dev.
## culture (Intercept) 0.0854605 0.29234
## culture.1 religiosity 0.0006695 0.02587
## Residual 0.2661269 0.51587
## Number of obs: 4719, groups: culture, 17
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.94883 0.18670 13.72800 21.150 7.17e-12 ***
## religiosity 0.01802 0.02595 18.53300 0.694 0.4960
## prop_rel -0.07762 0.24263 14.59300 -0.320 0.7536
## religiosity:prop_rel 0.11855 0.03294 17.71300 3.599 0.0021 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) rlgsty prp_rl
## religiosity -0.142
## prop_rel -0.916 0.145
## rlgsty:prp_ 0.146 -0.930 -0.185anova(slope.outcome2, slope.outcome.varcomp, slope.outcome.intercept, refit = FALSE)## Data: relfam3
## Models:
## ..2: family ~ religiosity + prop_rel + religiosity * prop_rel + (1 |
## ..2: culture)
## ..1: family ~ religiosity + prop_rel + religiosity * prop_rel + (1 |
## ..1: culture) + (religiosity - 1 | culture)
## object: family ~ religiosity + prop_rel + religiosity * prop_rel + (religiosity |
## object: culture)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..2 6 7259.9 7298.6 -3623.9 7247.9
## ..1 7 7258.6 7303.8 -3622.3 7244.6 3.3031 1 0.0691472 .
## object 8 7249.1 7300.8 -3616.6 7233.1 11.4374 1 0.0007198 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Der Vergleich zeigt, dass der p-Wert für \(\boldsymbol {\sigma^2_{\upsilon_1}}\) nach Halbierung (wegen bounded parameter space, s.o.) nach wie vor signifikant ist (p = 0.0346), und dass auch die Kovarianz \(\boldsymbol {\sigma_{\upsilon0\upsilon1}}\) nach wie vor signifikant ist (p = 0.0007).
Bates, D. (2013). Fitting linear mixed models in R. Using the lme4 package. R Newsletter, 5,, 27-30.
Eid, M., Gollwitzer, M., & Schmidt, M. (2013). Statistik und Forschungsmethoden (3. Aufl.). Weinheim: Beltz.
Mayer, B., & Trommsdorff, G. (2012). Cross-cultural perspectives on adolescents’ religiosity and family orientation. In G. Trommsdorff & X. Chen (Eds.), Values, religion, and culture in adolescent development (pp. 341-369). New York: Cambridge University Press.
West, B. T., Welch, K. B., & Galecki, A. T. (2015). Linear mixed models. A practical guide using statstical software (2nd ed.). Boca Raton: CRC Press.