Willkommen bei Stats by Randolph. Hier geht es zur Hauptseite mit weiteren Informationen und Inhalten.
Autor:innen dieser Seite: An den Inhalten dieser Seite haben mitgearbeitet: Markus Janczyk und Valentin Koob. Der Inhalt dieser Seite wird in der Lehre in den Studiengängen Psychologie von der AG Psychologische Forschungsmethoden und Kognitive Psychologie an der Universität Bremen verwendet, steht aber allen Interessierten zur Verfügung. Rückmeldungen/Fehler/Vorschläge können gesendet werden an randolph@uni-bremen.
Versionshistory:
Wir beginnen mit der Simulation eines Datensatzes, so wie wir es im Haupttext auch getan haben:
set.seed(1) # um die Simulation reproduzierbar zu machen
n_subjects <- 50 # Wieviele Versuchspersonen?
n_items <- 50 # Wieviele Items?
conditions <- c(-0.5, 0.5) # Wir benennen die beiden Bedingungen mit -0.5 und 0.5
# damit der fixed effect nachher den generellen RT Unterschied
# zwischen beiden Bedingungen darstellt. Die beiden Bedingungen
# könnten auch, wie im Beispiel, kongruente und inkongruente
# Zuordnungen von Stimuli zu Tasten sein
sigma_residual <- 20 # SD der Residuen
# fixed effects
beta_0 <- 800 # Intercept
beta_1 <- 100 # Slope
# Varianzkomponenten subject random factor
SD_0s <- 30 # Intercept
SD_1s <- 200 # Slope
COV_01s <- 4000 # Kovarianz Intercept-Slope
# Varianzkomponenten item random factor
SD_0i <- 20 # Intercept
SD_1i <- 40 # Slope
COV_01i <- 200 # Kovarianz Intercept-Slope
# Kovarianzmatrix für subject random factor
CovMatrixSubject <- matrix(c(SD_0s^2, COV_01s,
COV_01s, SD_1s^2),
nrow = 2, ncol = 2,
byrow = TRUE)
# Kovarianzmatrix für item random factor
CovMatrixItems <- matrix(c(SD_0i^2, COV_01i,
COV_01i, SD_1i^2),
nrow = 2, ncol = 2,
byrow = TRUE)
# Wir beginnen mit Intercept und Slope für den subject random factor
# simulieren der Daten gemäß der oben spezifizirten Kovarianzmatrix und mit einem
# Erwartungswert 0
randomEffects_Subjects <- MASS::mvrnorm(n = n_subjects,
mu = c(0, 0),
Sigma = CovMatrixSubject,
empirical = TRUE)
# Die simulierten Werte und die Labels für die Versuchspersonen in ein DataFrame...
randomEffects_Subjects <- data.frame(Subject = c(1:n_subjects),
round(randomEffects_Subjects))
# ...und passende Namen geben
names(randomEffects_Subjects)[2:3] <- c("Intercept_Subject", "Slope_Subject")
# Das gleiche machen wir dann für den item random factor
randomEffects_Items <- MASS::mvrnorm(n = n_items,
mu = c(0, 0),
Sigma = CovMatrixItems,
empirical = TRUE)
randomEffects_Items <- data.frame(Item = c(1:n_items),
round(randomEffects_Items))
names(randomEffects_Items)[2:3] <- c("Intercept_Item", "Slope_Item")
# Zunächst legen wir den DataFrame und entsprechende Variablen für Subject,
# Condition und Item an.
data <- data.frame(Subject = rep(c(1:n_subjects), each = 2*n_items),
Condition = c(-0.5, 0.5),
Item = rep(rep(c(1:n_items), each = 2), times = n_subjects))
# Dann fügen wir die Werte für die fixed effects hinzu, die ja für alle Versuchspersonen
# und Item gelten sollen.
data$beta_0 <- beta_0
data$beta_1 <- beta_1
# Nun fügen wir dem DataFrame die Intercepts und Slopes des item random factors hinzu...
data <- merge(data, randomEffects_Items,
by = "Item")
# ...und danach das Ganze noch für den subject random factor.
data <- merge(data, randomEffects_Subjects,
by = "Subject")
# Als letzter Schritt wird eine Variable mit den Residuen hinzugefügt.
data$residual <- round(rnorm(n = 2*n_subjects*n_items,
mean = 0,
sd = sigma_residual))
data$RT <- data$beta_0 + data$Condition*data$beta_1 +
data$Intercept_Item + data$Condition*data$Slope_Item +
data$Intercept_Subject + data$Condition*data$Slope_Subject +
data$residualInsgesamt sieht der Datensatz nun wie folgt aus:
car::some(data) Auf diese Daten wenden wir nun das Mixed-Model an:
M1 <- lmerTest::lmer(RT ~ Condition + (Condition | Subject) + (Condition | Item),
data = data)
summary(M1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RT ~ Condition + (Condition | Subject) + (Condition | Item)
## Data: data
##
## REML criterion at convergence: 45458.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4572 -0.6370 -0.0087 0.6635 3.7238
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 915.2 30.25
## Condition 39991.2 199.98 0.66
## Item (Intercept) 409.9 20.25
## Condition 1661.7 40.76 0.24
## Residual 420.4 20.50
## Number of obs: 5000, groups: Subject, 50; Item, 50
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 799.849 5.156 85.135 155.123 < 2e-16 ***
## Condition 98.755 28.869 53.040 3.421 0.00121 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Condition 0.565Wir sehen, dass hier eine Korrelation der fixed effects, also des Intercepts und des Slopes von \(0.565\) ausgegeben wird.
Diese Korrelation mag auf den ersten Blick komisch erscheinen, da die
fixed effects ja einzelne Werte, die auch von summary()
ausgegeben werden. Allerdings werden ja auch Standardfehler der fixed
effects ausgegeben, was also schon einmal bedeutet, dass die Schätzungen
nicht immer konstant sein können. Insbesondere würden wir, und das ist
auch der Hauptpunkt der Interpretation bereits, andere Schätzungen der
fixed effects erhalten, wenn wir das Mixed-Model an einer anderen
Zufallsstichprobe schätzen. Mit jeder solchen weiteren Schätzung
erhalten wir aber wieder einen Intercept und einen Slope und die
genannte Korrelation ist nichts anderes als die Schätzung der
Korrelation dieser Werte, wenn wir immer und immer wieder das Modell auf
eine neue Zufallsstichprobe anwenden.
Dies können wir uns anhand einer kleinen Simulation anschauen:
set.seed(1)
# Wir simulieren 100 Datensätze gemäß der Stuktur von M1. Das Resultat ist
# ein DataFrame mit 100 Spalten (für jede Simulation) und 5000 Zeilen (so wie
# im oben simulierten Ausgangsdatensatz)
simulierteWerte <- simulate(M1,
nsim = 100)
# Die folgende Funktion schätzt mit den simulierten Daten (eine Spalte von
# simulierteWerte) das Mixed-Model
LMM_Auswertung <- function(response)lme4::lmer(response ~ Condition +
(Condition | Subject) +
(Condition | Item),
data = data)
# Anwendung der Funktion auf das DataFrame (dauert ein wenig...)
ergebnisse <- lapply(X = sim,
FUN = LMM_Auswertung)
# In einem nächsten Schritt extrahieren wir die fixed effects der Mixed-Models...
fixedEffects <- lapply(X = ergebnisse,
FUN = lme4::fixef)
# ...und erstellen daraus eine Matrix mit den Werten für Intercept und Slope
fixedEffects <- do.call(rbind, fixedEffects)
head(fixedEffects) # Anzeigen der ersten sechs FälleNun können wir die Korrelation der beiden Spalten der Matrix berechnen und wir sehen, dass der erhaltene Wert in die Richtung des oben angezeigten Wertes geht:
cor(fixedEffects[,1], fixedEffects[,2])## [1] 0.5104151Tatsächlich ist eine Korrelation der fixed effects kein Spezifikum
von Mixed-Models, sondern die gleiche Logik gilt auch z.B. bei
Standardregressionen im Rahmen des ALM. Typischerweise werden sie dort
aber nicht von summary() ausgegeben (und bei Mixed-Models
könnte ihre Ausgabe mit corr = FALSE auch unterdrückt
werden).
Zur Illustration schauen wir uns die Situation im ALM daher auch kurz an. Wir simulieren daher 100 Werte, die als Prädiktor \(X\) und als Kriterium \(Y\) dienen sollen und mit \(r=0.5\) korreliert sind:
# Kovarianzmatrix die r = 0.5 impliziert
CovMatrix <- matrix(c(40, 20,
20, 40),
nrow = 2, ncol = 2,
byrow = TRUE)
# Simulation der Daten
simulierteDaten <- round(MASS::mvrnorm(n = 100,
mu = c(0, 0),
Sigma = CovMatrix,
empirical = TRUE))
simulierteDaten <- data.frame(simulierteDaten)
names(simulierteDaten) <- c("X","Y")
cor(simulierteDaten$X, simulierteDaten$Y) # tatsächliche Korrelation## [1] 0.4958333Im nächsten Schritt führen wir eine ganz normale Regression mit
lm() durch und extrahieren den geschätzte Wert der
Korrelation der beiden Parameter des Modells, also des Intercepts und
des Slopes (oder \(\hat{\beta_0}\) bzw.
\(\hat{\beta_1}\)):
ALM_Regression <- lm(Y ~ X,
data = simulierteDaten)
cov2cor(vcov(ALM_Regression))[1, 2] # geschätzte Korrelation der beiden Parameter## [1] 0.004814349Die Logik ist nun die gleiche wie oben beschrieben, daher fassen wir
uns hier etwas kürzer. Da die Berechnung mit lm() nicht
annähernd solange dauert wie im Fall eines Mixed-Models, simulieren wir
hier direkt 10000 Datensätze:
set.seed(1)
# Daten (Y) simulieren
simulierteWerte_LM <- simulate(ALM_Regression,
nsim = 10000)
# Funktion zur Auswertung
LM_Auswertung <- function(response)lm(response ~ X,
data = simulierteDaten)
ergebnisse <- lapply(X = simulierteWerte_LM,
FUN = LM_Auswertung)
# Koeffizienten extrahieren...
Koeffizienten <- lapply(X = ergebnisse,
FUN = coef)
# ...und in eine Matrix packen
Koeffizienten <- do.call(rbind, Koeffizienten)
head(Koeffizienten)## (Intercept) X
## sim_1 0.5705466 0.5809363
## sim_2 -0.2455490 0.4622317
## sim_3 0.1293851 0.5008570
## sim_4 0.2495444 0.4575279
## sim_5 -0.2529207 0.4613596
## sim_6 -0.2808386 0.5250193Nun können wir uns die Korrelation ausrechnen und sehen, dass die simulierte Korrelation in die Richtung der vorab aus dem Modell extrahierten Schätzung geht:, also quasi keine Korrelation der Parameter vorliegt:
cor(Koeffizienten[,1], Koeffizienten[,2])## [1] -0.00078855