Gemischte Modelle
Auf dieser Seite findest du Übungen zu Gemischten Modellen in R.
Der gesamte Code, den du hier findest, kann von dir 1:1 in R umgesetzt werden. Die größte Lernwirkung hat diese Einführung, wenn du möglichst alles einmal selbst in R durchführst.
01 Vorbereitungen
Bevor wir loslegen können, müssen wir die Packages laden, die wir für unser Vorhaben benötigen. Neben dem SfL Package laden wir zusätzliche Packages, mit denen wir Gemischte Modelle erstellen können:
install.packages("lme4")
install.packages("lmerTest")
install.packages("LMERConvenienceFunctions")library(SfL)
library(lme4)
library(lmerTest)
library(LMERConvenienceFunctions)Hinweis: Solltest du im späteren Verlauf feststellen, dass dir keine p-Werte angezeigt werden, führe folgenden Schritt zur manuellen Installation einer früheren Version des Packages durch:
require(devtools)
install_version("lme4", version = "1.1.23", repos = "http://cran.us.r-project.org")Führe diesen Schritt wirklich nur dann aus, wenn du später keine p-Werte angezeigt bekommst!
Nun laden wir die Datensätze, die wir nutzen werden:
data(data_a)
data(data_c)
data(data_s)
data(data_t)
data(data_v)02 Einführung
Die Nutzung von Gemischten Modellen ist im Vergleich zu anderen Methoden aufwendig. Aber lohnt sich der Aufwand?
Hierzu werfen wir einen Blick auf die Beispiel-Daten zu Kindern und ihrem Wachstum aus dem Theorie-Teil.
height_s <- lm(height ~ age, data_h) # simple
height_m <- lm(height ~ age + bsex, data_h) # multiple
height_l1 <- lmer(height ~ age + bsex + (1|name), data_h, REML=F) # mixed w/ random intercept
height_l2 <- lmer(height ~ age + bsex + (age|name), data_h, REML=F) # mixed w/ random intercept + random slopeHier haben wir nun also ein Simples Lineares Modell, ein Multiples Lineares Modell, ein Gemischtes Modell mit Random Intercept und ein Gemischtes Modell mit Random Intercecpt sowie Random Slope. Wie schneiden diese Modelle ab?
AIC(height_s, height_m, height_l1, height_l2) df AIC
height_s 3 284.4112
height_m 4 279.4706
height_l1 5 255.5833
height_l2 7 248.3431Wir wir sehen nimmt der AIC-Wert von oben nach unten ab. Demnach werden die Modelle mit zunehmender Komplexität besser, und das Gemischte Modell mit Random Intercept und Random Slope beschreibt die Daten am besten.
03 Vorbereitung der Daten
Wir möchten sDur anhand der folgenden Variablen modelln:
baseDur, speakingRate, googleFreq, age, speaker, typeOfS, folSeg, folType, item, transcription, biphoneProbSum, biphoneProbSumBin, monoMultilingual, pauseBin, preC
Hierzu müssen wir zunächst einige Dinge überprüfen.
Distribution Check
Wie wir bereits wissen, ist die Verteilung der Daten ein entscheidender Punkt. Daher checken wir die Verteilung der folgenden numerischen Variablen mit Shapiro-Wilk Tests:
sDur, baseDur, speakingRate, googleFreq
sDur
shapiro.test(data_s$sDur)
Shapiro-Wilk normality test
data: data_s$sDur
W = 0.93622, p-value = 2.761e-06Mit \(p=2.761e^{-6}<0.05\), ist sDur nicht normalverteilt.
data_s$sDurLog <- log(data_s$sDur)
shapiro.test(data_s$sDurLog)
Shapiro-Wilk normality test
data: data_s$sDurLog
W = 0.9927, p-value = 0.6438Mit \(p=0.6438>0.05\), ist sDurLog normalverteilt.
baseDur
shapiro.test(data_s$baseDur)
Shapiro-Wilk normality test
data: data_s$baseDur
W = 0.90932, p-value = 4.606e-08Mit \(p=4.606e^{-8}<0.05\), ist baseDur nicht normalverteilt.
data_s$baseDurLog <- log(data_s$baseDur)
shapiro.test(data_s$baseDurLog)
Shapiro-Wilk normality test
data: data_s$baseDurLog
W = 0.98401, p-value = 0.07966Mit \(p=0.07966>0.05\), ist baseDurLog normalverteilt.
speakingRate
shapiro.test(data_s$speakingRate)
Shapiro-Wilk normality test
data: data_s$speakingRate
W = 0.96027, p-value = 0.0002596Mit \(p=0.0003<0.05\), ist speakingRate nicht normalverteilt.
data_s$speakingRateLog <- log(data_s$speakingRate)
shapiro.test(data_s$speakingRateLog)
Shapiro-Wilk normality test
data: data_s$speakingRateLog
W = 0.99308, p-value = 0.6884Mit \(p=0.6884>0.05\), ist speakingRateLog normalverteilt.
googleFreq
shapiro.test(data_s$googleFreq)
Shapiro-Wilk normality test
data: data_s$googleFreq
W = 0.3621, p-value < 2.2e-16Mit \(p=2.2e^{-16}<0.05\), ist googleFreq nicht normalverteilt.
data_s$googleFreqLog <- log(data_s$googleFreq)
shapiro.test(data_s$googleFreqLog)
Shapiro-Wilk normality test
data: data_s$googleFreqLog
W = 0.97312, p-value = 0.004867Mit \(p=0.004867<0.05\), googleFreqLog nicht normalverteilt. Allerdings ist die Verteilung von googleFreqLog dennoch näher an einer Normalverteilung:
Correlation Check
Wie wir bereits wissen, ist Kollinearität ein Problem für Modelle mit vielen Prediktorvariablen. Daher testen wir zunächst die Korrelation der Variablen.
correlation_matrix(data = data_s,
variables = c("baseDurLog", "speakingRateLog", "googleFreqLog",
"age", "speaker", "typeOfS", "folSeg", "folType", "item", "transcription", "biphoneProbSum", "biphoneProbSumBin", "monoMultilingual", "pauseBin", "preC"))
Betrachten wir die Korrelationskoeffizienten näher, finden wir folgende problematische Fälle:
folSegvs.folTypebiphoneProbSumBinvs.biphoneProbSumbiphoneProbSumBinvs.itembiphoneProbSumBinvs.transcriptionsbiphoneProbSumvs.itembiphoneProbSumvs.transcriptionitemvs.transcription
Im nächsten Schritt nutzen wir die “Vergleich und Ausschluss”-Strategie zur Vermeidung von Kollinearitätsproblemen.
Vergleich und Ausschluss
Idee: Die stark korrelierten Variablen werden in simple Modelle eingesetzt und diese dann vergliechen. Die Variable mit der jeweils größeren Vorhersagekraft wird behalten.
Potentielles Problem: Der Verlust von Variablen kann einem Verlust an Vorhersagekraft entsprechen.
Vorgehen: Mit der predictor_competition() Funktion des SfL Packages kann herausgefunden werden, welche Variable der jeweils bessere Prediktor ist. Alternativ kann dieses Vorgehen auch manuell durchgeführt werden.
Da wir diesmal mit Gemischten Modellen arbeiten, spezifizieren wir im predictor_competition() Befehl auch einen Random Intercept.
Schritt 1
folSeg vs. folType
predictor_competition(data = data_s,
dependent = "sDurLog",
independent1 = "folSeg",
independent2 = "folType",
random.intercept = "speaker")Data: data
Models:
mdl2: sDurLog ~ folType + (1 | speaker)
mdl1: sDurLog ~ folSeg + (1 | speaker)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
mdl2 7 135.74 156.81 -60.869 121.74
mdl1 17 140.86 192.04 -53.431 106.86 14.877 10 0.1366Wie vom Output angedeutet, sind beide Variablen gleich gute Prediktoren. In diesem Fall können wir wählen, welchen wir für die weitere Analyse behalten möchten.
Wir behalten folType. Weitere Vergleiche mit folSeg entfallen.
Schritt 2
biphoneProbSumBin vs. biphoneProbSum
predictor_competition(data = data_s,
dependent = "sDurLog",
independent1 = "biphoneProbSumBin",
independent2 = "biphoneProbSum",
random.intercept = "speaker")Data: data
Models:
mdl1: sDurLog ~ biphoneProbSumBin + (1 | speaker)
mdl2: sDurLog ~ biphoneProbSum + (1 | speaker)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
mdl1 4 130.32 142.36 -61.161 122.32
mdl2 4 129.73 141.77 -60.863 121.73 0.5955 0 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Wie vom Output angedeutet, ist eine Variable ein besserer Prediktor als die andere Variable. Anhand des AIC-Wertes stellen wir fest, die bessere Prediktorvariable ist biphoneProbSum.
Wir behalten biphoneProbSum. Weitere Vergleiche mit biphoneProbSumBin entfallen.
Schritt 3
item vs. transcription
Achtung: Diese Variablen möchten wir als Random Effects einsetzen. Daher erfolgt der Vergleich mit etwas anderen Mitteln.
mdl_item <- lmer(sDurLog ~ 1 + (1|item), data = data_s, REML = F)
mdl_transcription <- lmer(sDurLog ~ 1 + (1|transcription), data = data_s, REML = F)
anova(mdl_item, mdl_transcription)Data: data_s
Models:
mdl_item: sDurLog ~ 1 + (1 | item)
mdl_transcription: sDurLog ~ 1 + (1 | transcription)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
mdl_item 3 162.51 171.55 -78.257 156.51
mdl_transcription 3 162.53 171.57 -78.268 156.53 0 0 1Wie vom Output angedeutet, sind beide Variablen gleich gute Prediktoren. In diesem Fall können wir wählen, welchen wir für die weitere Analyse behalten möchten.
Wir behalten transcription. Weitere Vergleiche mit item entfallen.
Schritt 4
biphoneProbSum vs. transcription
Achtung: Da wir eine der beiden Variablen als Fixed Effect und die andere als Random Effect nutzen möchten, müssen wir keine Kollinearitätsprobleme befürchten. Entsprechend muss keine der beiden Variablen aussortiert werden.
Interim Summary
Zu Beginn haben wir problematische Korrelationen zwischen folgenden Variablen gefunden:
folSegvs.folTypebiphoneProbSumBinvs.biphoneProbSumbiphoneProbSumBinvs.itembiphoneProbSumBinvs.transcriptionsbiphoneProbSumvs.itembiphoneProbSumvs.transcriptionitemvs.transcription
Nach dem Anwenden der “Vergleich und Auschluss”-Strategie sind folgende Variablen erhalten geblieben: folType, biphoneProbSum, und transcription.
04 Volles Modell
Mit den Variablen, die das Ausschlussverfahren überlebt haben, sowie mit den Variablen, die keine problematischen Korrelationen gezeigt haben, kann nun ein volles Modell erstellt werden:
mdl_full <- lmer(sDurLog ~
folType +
biphoneProbSum +
baseDurLog +
speakingRateLog +
googleFreqLog +
age +
typeOfS +
monoMultilingual +
pauseBin +
preC +
(1 | speaker) +
(1 | transcription),
data_s, REML=F)05 Bestes Modell
Mit der step() Funktion können wir erneut das beste Modell finden:
step(mdl_full)Backward reduced random-effect table:
Eliminated npar logLik AIC LRT Df Pr(>Chisq)
<none> 20 -6.3749 52.750
(1 | transcription) 1 19 -6.3752 50.750 0.0008 1 0.9781
(1 | speaker) 0 18 -15.5262 67.052 18.3020 1 1.885e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Backward reduced fixed-effect table:
Degrees of freedom method: Satterthwaite
Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
preC 1 0.16862 0.05621 3 124.653 1.1940 0.314902
age 2 0.04662 0.04662 1 28.033 0.9408 0.340359
googleFreqLog 3 0.11108 0.11108 1 125.054 2.2381 0.137166
folType 4 0.42808 0.10702 4 133.867 2.0869 0.085967 .
biphoneProbSum 0 0.40410 0.40410 1 131.726 7.3133 0.007749 **
baseDurLog 0 1.45686 1.45686 1 144.633 26.3656 8.971e-07 ***
speakingRateLog 0 0.25002 0.25002 1 146.146 4.5247 0.035087 *
typeOfS 0 1.21841 0.60921 2 136.646 11.0251 3.640e-05 ***
monoMultilingual 0 0.25425 0.25425 1 28.866 4.6012 0.040496 *
pauseBin 0 1.70080 1.70080 1 137.036 30.7803 1.440e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Model found:
sDurLog ~ biphoneProbSum + baseDurLog + speakingRateLog + typeOfS + monoMultilingual + pauseBin + (1 | speaker)Zuerst stellen wir fest, dass transcription als Random Intercept nicht benötigt wird.
Dann sehen wir, dass preC, age, googleFreqLog und folType keine signifikanten Prediktoren für sDurLog sind, und demnach aussortiert werden können.
Das beste Modell ist demnach folgendes:
mdl_best <- lmer(sDurLog ~
biphoneProbSum +
baseDurLog +
speakingRateLog +
typeOfS +
monoMultilingual +
pauseBin +
(1 | speaker),
data_s, REML=F)06 Diagnostik / Trimming
Gemischte Modelle folgenden den gleichen Assumptions, denen auch Simple und Multiple Lineare Regressionsmodelle folgen. Da wir bereits einige Male gesehen haben, wie man diese überprüft, sparen wir uns diesen Schritt hier.
Stattdessen schauen wir uns an, wie das finale Modell noch ein wenig besser gemacht werden kann: durch Residual Trimming.
Durch Residual Trimming wird der “Fit” des Modells verbessert, da Datenpunkte entfernt werden, die durch das Modell nur sehr schlecht erklärt werden können. Anders als das Aussortierten von Outliern vor dem Modell-Prozess ist dieses Vorgehen statistisch motiviert und entspricht modernen Konventionen.
Nachdem man, so wie im vorliegenden Fall, ein “bestes” Modell gefunden hat, kann man mit dem LMERConvenienceFunctions Package die Daten trimmen:
library(LMERConvenienceFunctions)
data_s2 <- romr.fnc(mdl_best, # das beste Modell
data_s, # der Datensatz
trim = 2.5) # die Standardabweichung mit welcher man trimmen möchte; 2.5 ist die Standardeinstellungn.removed = 1
percent.removed = 0.6666667 Mit diesem Befehl trimmt man die Daten, deren Residuen im Modell \(2.5\) oder mehr Standardabweichungen vom Durchschnitt der Gesamtresiduen entfernt liegen.
Im vorliegenden Fall wird lediglich \(1\) Datenpunkt getrimmt, d.h. \(0.67\%\) der gesamten Datenmenge.
Anschließend wird das getrimmte Datenset im Environment gesichert:
data_s2 <- data_s2$dataZuletzt wird das beste Modell mit dem getrimmten Datensatz erneut erstellt:
mdl_trim <- lmer(sDurLog ~ biphoneProbSum +
baseDurLog +
speakingRateLog +
typeOfS +
monoMultilingual +
pauseBin +
(1 | speaker),
data_s2, REML=F)Betrachten wir nun die AIC-Werte der beiden “besten” Modelle, sehen wir eine Verbesserung durch die getrimmten Daten:
AIC(mdl_best, mdl_trim) df AIC
mdl_best 10 47.30116
mdl_trim 10 39.3118407 Interpretation
Zu guter Letzt können wir uns nun an die Interpretation des Gemischten Modells begeben.
Signifikanz & Predictor Strength
Wie zuvor können wir mit der anova() Funktion die Signifikanz von Variablen überprüfen:
anova(mdl_trim)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
biphoneProbSum 0.37715 0.37715 1 129.745 7.2917 0.007852 **
baseDurLog 1.48311 1.48311 1 145.405 28.6737 3.272e-07 ***
speakingRateLog 0.27047 0.27047 1 144.081 5.2291 0.023669 *
typeOfS 1.32953 0.66477 2 134.213 12.8522 7.818e-06 ***
monoMultilingual 0.22935 0.22935 1 28.915 4.4341 0.044028 *
pauseBin 1.72387 1.72387 1 134.739 33.3284 5.116e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Offenbar sind einige Prediktoren stärker als andere, da die p-Werte insgesamt sehr unterschiedlich sind. Um ein genaues Bild von diesen Unterschieden zu erhalten, können wir die predictor_strength()-Funktion aus dem SfL Package nutzen:
strength <- predictor_strength(dependent = "sDurLog",
fixed = c("biphoneProbSum", "baseDurLog", "speakingRateLog", "typeOfS", "monoMultilingual", "pauseBin"),
random_str = c("( 1 | speaker)"),
data = data_s2)
strength[order(strength$R2m),] predictor R2m R2c
5 baseDurLog 0.3637483 0.6439048
3 typeOfS 0.4591623 0.6394286
1 pauseBin 0.4594154 0.6005595
6 biphoneProbSum 0.4979181 0.6737540
2 monoMultilingual 0.5009224 0.6798167
4 speakingRateLog 0.5147190 0.6686202Was bedeuten diese Werte?
R2m ist der s.g. “marginal R-squared”-Wert. Dieser Wert gibt an, wie viel Varianz in den Daten anhand der Fixed Effects im Modell erklärt wird. Je höher dieser Wert ist, desto mehr Varianz wird vom Modell und seinen Prediktoren erklärt.
R2c ist der s.g. “conditional R-squared”-Wert. Dieser Wert gibt an, wie viel Varianz in den Daten anhand der Fixed Effects + Random Effects im Modell erklärt wird. Je höher dieser Wert ist, desto mehr Varianz wird vom Modell und seinen Prediktoren erklärt.
Wenn wir nun also einen Fixed Effect Prediktor aus dem Modell berechnen und einen großen Verlust beim R2m-Wert feststellen, hat der entfernte Prediktor eine hohe Predictor Strength. Ist der Verlust gering, ist der Prediktor eher schwach.
Genau diese Idee setzt die predictor_strength()-Funktion um. Der Output ist oben bereits durch die letzte Code-Zeile nach Stärke sortiert worden. Wir sehen, dass baseDurLog der stärkste Prediktor ist, während speakingRateLog der schwächste Prediktor ist.
Koeffizienten
Mit der summary()-Funktionen können wir uns wie gewohnt die Koeffizienten des Modells anschauen:
summary(mdl_trim)Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: sDurLog ~ biphoneProbSum + baseDurLog + speakingRateLog + typeOfS + monoMultilingual + pauseBin + (1 | speaker)
Data: data_s2
AIC BIC logLik deviance df.resid
39.3 69.4 -9.7 19.3 139
Scaled residuals:
Min 1Q Median 3Q Max
-2.0987 -0.5820 -0.0124 0.6397 2.2855
Random effects:
Groups Name Variance Std.Dev.
speaker (Intercept) 0.02421 0.1556
Residual 0.05172 0.2274
Number of obs: 149, groups: speaker, 38
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -1.63149 0.16998 104.48189 -9.598 5.11e-16 ***
biphoneProbSum 8.04380 2.97884 129.74520 2.700 0.00785 **
baseDurLog 0.58748 0.10971 145.40485 5.355 3.27e-07 ***
speakingRateLog -0.20866 0.09125 144.08078 -2.287 0.02367 *
typeOfSnm 0.25994 0.05175 132.20258 5.023 1.61e-06 ***
typeOfSpl 0.10968 0.05327 135.24495 2.059 0.04140 *
monoMultilingualmonolingual 0.16531 0.07851 28.91453 2.106 0.04403 *
pauseBinpause 0.26015 0.04506 134.73859 5.773 5.12e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) bphnPS bsDrLg spknRL typOfSn typOfSp mnMltl
biphonPrbSm -0.134
baseDurLog 0.654 0.035
speakngRtLg -0.360 -0.074 0.294
typeOfSnm -0.120 -0.098 -0.097 -0.085
typeOfSpl -0.139 -0.013 -0.200 -0.237 0.526
mnMltlnglmn -0.354 0.010 -0.058 -0.075 -0.081 -0.018
pauseBinpas -0.339 -0.028 -0.252 0.067 -0.051 -0.094 0.051Die Interpretation der Koeffizienten funktioniert genau so wie bei Multiplen Linearen Modellen.
Tukey Contrasts
Tukey Contrasts werden genutzt um die Level von kategorischen Prediktoren zu untersuchen, d.h. um zu schauen, ob diese signifikant verschieden zueinander sind.
tukey(model = mdl_trim,
predictor = typeOfS)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lmer(formula = sDurLog ~ biphoneProbSum + baseDurLog + speakingRateLog +
typeOfS + monoMultilingual + pauseBin + (1 | speaker), data = data_s2,
REML = F)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
nm - is == 0 0.25994 0.05175 5.023 < 0.001 ***
pl - is == 0 0.10968 0.05327 2.059 0.09844 .
pl - nm == 0 -0.15025 0.05114 -2.938 0.00919 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)tukey(model = mdl_trim,
predictor = monoMultilingual)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lmer(formula = sDurLog ~ biphoneProbSum + baseDurLog + speakingRateLog +
typeOfS + monoMultilingual + pauseBin + (1 | speaker), data = data_s2,
REML = F)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
monolingual - bilingual == 0 0.16531 0.07851 2.106 0.0352 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)tukey(model = mdl_trim,
predictor = pauseBin)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lmer(formula = sDurLog ~ biphoneProbSum + baseDurLog + speakingRateLog +
typeOfS + monoMultilingual + pauseBin + (1 | speaker), data = data_s2,
REML = F)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
pause - no_pause == 0 0.26015 0.04506 5.773 7.78e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)08 Aufgaben
Nutze die anderen Datensätze um Gemischte Modelle zu erstellen. Zum Beispiel könntest du folgendes Modell erstellen:
Aus data_v:
abhängige Variable:
vowel_durPrediktoren:
subject,gender,item,voicing,vowel,prosodic_boundary
Führe alle Schritte zur Modellerstellung durch, die dir sinnvoll erscheinen.
Alternativ kannst du dich aber zunächst auch an der Visualisierung von Gemischten Modellen versuchen. Hierzu verwendet man ebenfalls das visreg Package, welches wir bereits zur Visualisierung Linearer Modelle kennengelernt haben.