Building Intuitions about Mixed-Effects Models
Martin Zettersten
25 February, 2025
This code is adapted from Violet Brown’s introduction to mixed-effects modeling:
Brown VA. An Introduction to Linear Mixed-Effects Modeling in R. Advances in Methods and Practices in Psychological Science. 2021;4(1). doi:10.1177/2515245920960351
if (!("knitr" %in% installed.packages())) install.packages("knitr")
if (!("lme4" %in% installed.packages())) install.packages("lme4")
if (!("lmerTest" %in% installed.packages())) install.packages("lmerTest")
if (!("tidyverse" %in% installed.packages())) install.packages("tidyverse")
if (!("here" %in% installed.packages())) install.packages("here")Load libraries and set up paths
library(knitr)
library(tidyverse)
library(lme4)
library(lmerTest)
library(here)
data_path <- here("data")Fixed-effects only, random intercepts, and random slopes plots
Load data
figuredata <- read_csv(here(data_path,"figure_data.csv"))## Rows: 16 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (3): PID, xvar, yvar
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Make PID a factor
figuredata$PID <- as.factor(figuredata$PID)Fixed-effects only regression plot
Build regression model and view the summary output to look at the residuals
Intercept only model
ols.mod_intercept <- lm(yvar ~ 1, data = figuredata)
summary(ols.mod_intercept)##
## Call:
## lm(formula = yvar ~ 1, data = figuredata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -865.6 -421.9 -165.6 396.9 1034.4
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1165.6 143.3 8.134 7.01e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 573.2 on 15 degrees of freedomUncentered predictor (simple regression)
ols.mod <- lm(yvar ~ xvar, data = figuredata)
summary(ols.mod)##
## Call:
## lm(formula = yvar ~ xvar, data = figuredata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -721.36 -233.80 3.13 361.59 561.08
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 582.92 180.88 3.223 0.00614 **
## xvar 119.84 30.64 3.911 0.00157 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 410.2 on 14 degrees of freedom
## Multiple R-squared: 0.5221, Adjusted R-squared: 0.4879
## F-statistic: 15.29 on 1 and 14 DF, p-value: 0.001568Centered predictor (what changes?)
figuredata <- figuredata %>%
mutate(
xvar_c = xvar - mean(xvar)
)
ols.mod_c <- lm(yvar ~ xvar_c, data = figuredata)
summary(ols.mod_c)##
## Call:
## lm(formula = yvar ~ xvar_c, data = figuredata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -721.36 -233.80 3.13 361.59 561.08
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1165.62 102.55 11.367 1.87e-08 ***
## xvar_c 119.84 30.64 3.911 0.00157 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 410.2 on 14 degrees of freedom
## Multiple R-squared: 0.5221, Adjusted R-squared: 0.4879
## F-statistic: 15.29 on 1 and 14 DF, p-value: 0.001568figuredata <- figuredata %>%
mutate(
yvar_adj = yvar - 700
)
ols.mod_c <- lm(yvar_adj ~ xvar_c, data = figuredata)
summary(ols.mod_c)##
## Call:
## lm(formula = yvar_adj ~ xvar_c, data = figuredata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -721.36 -233.80 3.13 361.59 561.08
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 465.63 102.55 4.541 0.000462 ***
## xvar_c 119.84 30.64 3.911 0.001568 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 410.2 on 14 degrees of freedom
## Multiple R-squared: 0.5221, Adjusted R-squared: 0.4879
## F-statistic: 15.29 on 1 and 14 DF, p-value: 0.001568Plot
Build a fixed effects only plot
ggplot(figuredata, aes(x = xvar, y = yvar)) +
stat_smooth(method = lm, se = FALSE, linetype = "solid",
color = "black", linewidth = .6) +
geom_point(aes(shape = PID), size = 3.25, color = "grey70") +
scale_shape_manual(values = c(15, 16, 17, 18)) +
geom_segment(aes(x = xvar, xend = xvar,
y = yvar, yend = fitted(ols.mod)),
color = "grey70") +
scale_y_continuous(expand = c(0, 0), breaks = c(0, 750, 1500, 2250, 3000),
limits = c(0, 2600)) +
scale_x_continuous(expand = c(0, 0), breaks = c(0, 2, 4, 6, 8, 10),
limits = c(-0.5, 10.5)) +
theme(panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
legend.position = "none",
axis.text = element_text(size = 14),
axis.title = element_text(size = 14)) +
labs (x = "Word Difficulty", y = "Response Time") ## `geom_smooth()` using formula = 'y ~ x'
Random intercepts plot
Build the model with random intercepts and view the summary output to look at the residuals
random_intercepts.mod <- lmer(yvar ~ 1 + xvar + (1|PID), data = figuredata)
summary(random_intercepts.mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: yvar ~ 1 + xvar + (1 | PID)
## Data: figuredata
##
## REML criterion at convergence: 210.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6361 -0.6381 0.3516 0.6020 1.3016
##
## Random effects:
## Groups Name Variance Std.Dev.
## PID (Intercept) 108442 329.3
## Residual 75324 274.5
## Number of obs: 16, groups: PID, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 576.362 204.356 5.033 2.820 0.036817 *
## xvar 121.185 20.507 11.001 5.909 0.000102 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## xvar -0.488Extract the fixed effects estimates for the intercept and slope
model_intercept <- as.numeric(fixef(random_intercepts.mod)[1])
model_slope <- as.numeric(fixef(random_intercepts.mod)[2])Extract the individual participant intercepts for this model and add it to the data frame
figuredata$intercepts <- rep(coef(random_intercepts.mod)$PID[,1], each = 4)Build random intercepts plot
ggplot(figuredata, aes(x = xvar, y = yvar)) +
geom_abline(slope = model_slope, intercept = model_intercept,
linetype = "solid", color = "black", size = 1) +
geom_abline(mapping = aes(slope = model_slope, intercept = intercepts,color=PID),
linetype = "dashed", size = .4) +
geom_point(aes(shape = PID,color=PID), size = 3.25) +
scale_shape_manual(values = c(15, 16, 17, 18)) +
geom_segment(aes(x = xvar, xend = xvar,
y = yvar, yend = fitted(random_intercepts.mod),color=PID)) +
scale_y_continuous(expand = c(0, 0), breaks = c(0, 500, 1000, 1500, 2000, 2500),
limits = c(0, 2600)) +
scale_x_continuous(expand = c(0, 0), breaks = c(0, 2, 4, 6, 8, 10),
limits = c(-0.5, 10.5)) +
theme(panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
legend.position = "none",
axis.text = element_text(size = 14),
axis.title = element_text(size = 14)) +
labs (x = "Word Difficulty", y = "Response Time") ## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Random intercepts and slopes plot
Build the model with random intercepts and slopes and view the summary output to look at the residuals
random_slopes.mod <- lmer(yvar ~ 1 + xvar + (1 + xvar|PID), data = figuredata)
summary(random_slopes.mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: yvar ~ 1 + xvar + (1 + xvar | PID)
## Data: figuredata
##
## REML criterion at convergence: 191.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.0716 -0.4797 -0.1726 0.6937 1.0952
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## PID (Intercept) 61192 247.37
## xvar 5854 76.51 -0.40
## Residual 5638 75.08
## Number of obs: 16, groups: PID, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 561.105 128.055 2.990 4.382 0.0222 *
## xvar 124.636 38.667 2.997 3.223 0.0485 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## xvar -0.415Extract the individual participant intercepts and slopes from this model and add them to the data frame
figuredata$intercepts2 <- rep(coef(random_slopes.mod)$PID[,1], each = 4)
figuredata$slopes <- rep(coef(random_slopes.mod)$PID[,2], each = 4)Build plot
ggplot(figuredata, aes(x = xvar, y = yvar)) +
geom_abline(slope = model_slope, intercept = model_intercept,
linetype = "solid", color = "black", size = 1) +
geom_abline(mapping = aes(slope = slopes,
intercept = intercepts2, linetype = PID,color=PID),
linetype = "dashed", size = .4) +
geom_point(aes(shape = PID,color=PID), size = 3.25) +
scale_shape_manual(values = c(15, 16, 17, 18)) +
geom_segment(aes(x = xvar, xend = xvar,
y = yvar, yend = fitted(random_slopes.mod),color=PID)) +
scale_y_continuous(expand = c(0, 0), breaks = c(0, 750, 1500, 2250),
limits = c(0, 2600)) +
scale_x_continuous(expand = c(0, 0), breaks = c(0, 2, 4, 6, 8, 10),
limits = c(-0.5, 10.5)) +
theme(panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(colour = "black", fill = NA),
legend.position = "none",
axis.text = element_text(size = 14),
axis.title = element_text(size = 14)) +
labs (x = "Word Difficulty", y = "Response Time") 