# Load libraries
library(lme4)
library(lmerTest)
library(ggplot2)
library(tidyr)
library(dplyr)
library(gtsummary)Workshop: Linear Mixed Models in R
1 Quick Background
- Linear Mixed Models (LMMs): LMMs extend linear regression to account for clustered (correlated) data, such as repeated measures on subjects. They include fixed effects (e.g., treatment, time) for population-level trends and random effects (e.g., random intercepts or slopes per subject) to model within-subject correlation and heterogeneity.
- Key Components:
- Fixed Effects: Coefficients for predictors like treatment (
trt) and time. - Random Effects: Typically
(1|id)for random intercepts (subject-specific baseline), or(time|id)for random slopes (individual time trends).
- Fixed Effects: Coefficients for predictors like treatment (
- Assumptions: Linearity, normality of residuals/random effects, homoscedasticity; handles missing data under MAR.
- Here: We’ll analyze strength gains over time by treatment group using the
lme4package. Outcome: strength (y); Predictor: time (days); Group:trt(1=increase reps, 2=increase weight); Clustering: subjectid.
The dataset has 37 subjects (16 in trt=1, 21 in trt=2), with up to 7 measurements (assuming time points at days 0,2,4,6,8,10,12 for illustration, based on 7 columns). Some missings (denoted .).
Note: Run in RStudio. Save the provided stren.txt file in your working directory.
2 Install and Load Packages
We need lme4 for LMMs, lmerTest for p-values, ggplot2 for plots, and tidyr/dplyr for data reshaping.
- What this does:
lme4::lmer()fits LMMs;lmerTestadds ANOVA-style tests. - Expected output: Libraries loaded without errors.
3 Load and Prepare the Dataset
Read the wide-format data and reshape to long format for analysis.
# Read wide data (columns: id, trt, y1, y2, y3, y4, y5, y6, y7)
stren_wide <- read.table("https://www.dropbox.com/s/f1n20of73hrispo/stren.txt?dl=1",
header=F, na.strings = ".",
col.names = c("id", "trt", "y1", "y2",
"y3", "y4", "y5", "y6", "y7"),
colClasses = c("integer",
"integer",
rep("numeric", 7)))
# Reshape to long: add time variable (assuming days 0,2,4,6,8,10,12)
stren_long <- stren_wide %>%
pivot_longer(cols = y1:y7, names_to = "time_str", values_to = "strength") %>%
mutate(time = case_when(
time_str == "y1" ~ 0,
time_str == "y2" ~ 2,
time_str == "y3" ~ 4,
time_str == "y4" ~ 6,
time_str == "y5" ~ 8,
time_str == "y6" ~ 10,
time_str == "y7" ~ 12
)) %>%
select(id, trt, time, strength) %>%
filter(!is.na(strength)) # Remove missings for now
# View summary
summary(stren_long) id trt time strength
Min. : 1.00 Min. :1.000 Min. : 0.000 Min. :74.00
1st Qu.:10.00 1st Qu.:1.000 1st Qu.: 2.000 1st Qu.:79.00
Median :19.00 Median :2.000 Median : 6.000 Median :82.00
Mean :18.88 Mean :1.561 Mean : 5.707 Mean :81.51
3rd Qu.:28.50 3rd Qu.:2.000 3rd Qu.: 9.000 3rd Qu.:84.00
Max. :37.00 Max. :2.000 Max. :12.000 Max. :91.00 table(stren_long$trt, stren_long$time)
0 2 4 6 8 10 12
1 16 15 16 15 15 13 15
2 21 21 20 21 19 17 15- What this does: Reads wide data, pivots to long (one row per measurement), defines
timein days (adjusted to 7 points), drops NAs. - Expected output: ~200-250 rows (some missings). Means: strength ~80-82 increasing over time; ~16 subjects trt=1, ~21 trt=2; measurements per time ~30-37.
4 Summarise and visualise the Data
Visualize trajectories with a spaghetti plot.
tbl_summary(stren_wide)| Characteristic | N = 371 |
|---|---|
| id | 19 (10, 28) |
| trt | |
| 1 | 16 (43%) |
| 2 | 21 (57%) |
| y1 | 80.00 (79.00, 83.00) |
| y2 | 81.0 (79.0, 84.0) |
| Unknown | 1 |
| y3 | 81.0 (79.0, 84.0) |
| Unknown | 1 |
| y4 | 82.0 (80.0, 84.0) |
| Unknown | 1 |
| y5 | 82.0 (80.0, 84.0) |
| Unknown | 3 |
| y6 | 82.5 (79.0, 83.0) |
| Unknown | 7 |
| y7 | 82.0 (80.0, 85.0) |
| Unknown | 7 |
| 1 Median (Q1, Q3); n (%) | |
# Spaghetti plot
ggplot(stren_long,
aes(x = time,
y = strength,
color = factor(trt),
group = id)) +
geom_line(alpha = 0.6) +
stat_summary(fun=mean,
geom="line",
show.legend = F,
lwd=3,
aes(group=trt)) +
# geom_smooth(aes(group = trt), method = "lm", se = TRUE) + #linear trend
labs(title = "Strength Over Time by Treatment Group",
x = "Time (Days)", y = "Strength",
color = "Treatment (1=Reps, 2=Weight)") +
theme_minimal()
- What this does: Plots individual trajectories (lines per id) and group trends (smoothed lines).
- Expected output: Plot shows upward trends; trt=2 may increase faster (steeper slope). Some lines missing points due to dropouts.
5 Fit a Basic Random Intercept Model
Model strength as function of time and trt, with random intercept per subject.
# Basic LMM: fixed time + trt, random intercept
model_ri <- lmer(strength ~ time + trt + (1 | id), data = stren_long)
# Compute predicted values (including random effects)
stren_long$pred_ri <- predict(model_ri)
# Spaghetti plot of predicted values
ggplot(stren_long, aes(x = time,
y = pred_ri,
color = factor(trt),
group = id)) +
geom_line(alpha = 0.7) +
geom_smooth(aes(group = trt), method = "lm", se = TRUE) + #linear trend
labs(title = "Predicted Strength Trajectories by Treatment Group",
x = "Time (Days)",
y = "Predicted Strength",
color = "Treatment (1=Reps, 2=Weight)") +
theme_minimal()
# Summary
summary(model_ri)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: strength ~ time + trt + (1 | id)
Data: stren_long
REML criterion at convergence: 877.1
Scaled residuals:
Min 1Q Median 3Q Max
-2.6996 -0.5657 0.0124 0.5911 3.0591
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 10.646 3.263
Residual 1.228 1.108
Number of obs: 239, groups: id, 37
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 78.63195 1.79934 35.22012 43.700 < 2e-16 ***
time 0.13805 0.01834 201.13186 7.528 1.7e-12 ***
trt 1.38358 1.09248 34.96059 1.266 0.214
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) time
time -0.062
trt -0.952 0.005# In table format
tbl_regression(model_ri)| Characteristic | Beta | 95% CI | p-value |
|---|---|---|---|
| time | 0.14 | 0.10, 0.17 | <0.001 |
| trt | 1.4 | -0.83, 3.6 | 0.2 |
| Abbreviation: CI = Confidence Interval | |||
#Useful if there were categorical variables with >2 categories
#as it would give us the omnibus test (global p-value)
anova(model_ri) # Type III tests with lmerTestType III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
time 69.564 69.564 1 201.132 56.6667 1.7e-12 ***
trt 1.969 1.969 1 34.961 1.6039 0.2137
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- What this does: Fixed effects: intercept, time slope, trt effect; random: subject-specific deviation from intercept.
- Expected output: Positive time coef ; trt coef ( positive for trt=2); random SD (between-subject variability).
6 Compare to Ordinary Linear Model (For Illustration)
Ignore clustering with lm() on long data.
# Ignore clustering: lm on long data
model_lm <- lm(strength ~ time + trt, data = stren_long)
summary(model_lm)
Call:
lm(formula = strength ~ time + trt, data = stren_long)
Residuals:
Min 1Q Median 3Q Max
-7.4748 -2.4141 0.0017 2.2142 9.0859
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 78.76442 0.78871 99.865 < 2e-16 ***
time 0.10984 0.05484 2.003 0.04633 *
trt 1.35517 0.43638 3.105 0.00213 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.347 on 236 degrees of freedom
Multiple R-squared: 0.05326, Adjusted R-squared: 0.04524
F-statistic: 6.639 on 2 and 236 DF, p-value: 0.001567# Compare coefficients/SEs
coef(summary(model_lm)) Estimate Std. Error t value Pr(>|t|)
(Intercept) 78.7644219 0.78870963 99.864918 4.609448e-195
time 0.1098411 0.05484071 2.002911 4.633064e-02
trt 1.3551741 0.43637883 3.105499 2.131917e-03coef(summary(model_ri)) Estimate Std. Error df t value Pr(>|t|)
(Intercept) 78.6319495 1.79934434 35.22012 43.700335 2.821369e-32
time 0.1380451 0.01833821 201.13186 7.527731 1.699711e-12
trt 1.3835818 1.09247558 34.96059 1.266465 2.137201e-01- Expected: lm SEs smaller (underestimates uncertainty); biased p-values. LMM SEs larger, more conservative.
7 Add Treatment-Time Interaction
Allow different slopes by treatment.
# LMM with interaction
model_int <- lmer(strength ~ time * trt + (1 | id), data = stren_long)
# Compute predicted values (including random effects)
stren_long$pred_int <- predict(model_int)
# Spaghetti plot of predicted values
ggplot(stren_long, aes(x = time,
y = pred_int,
color = factor(trt),
group = id)) +
geom_line(alpha = 0.7) +
geom_smooth(aes(group = trt), method = "lm", se = TRUE) + #linear trend
labs(title = "Predicted Strength Trajectories by Treatment Group",
x = "Time (Days)",
y = "Predicted Strength",
color = "Treatment (1=Reps, 2=Weight)") +
theme_minimal()
# Summary
summary(model_int)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: strength ~ time * trt + (1 | id)
Data: stren_long
REML criterion at convergence: 881.2
Scaled residuals:
Min 1Q Median 3Q Max
-2.5979 -0.5463 0.0335 0.6079 3.0065
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 10.64 3.263
Residual 1.23 1.109
Number of obs: 239, groups: id, 37
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 78.90221 1.82853 37.55998 43.151 <2e-16 ***
time 0.09088 0.05975 200.07996 1.521 0.130
trt 1.20998 1.11229 37.56560 1.088 0.284
time:trt 0.03056 0.03684 200.11978 0.829 0.408
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) time trt
time -0.188
trt -0.954 0.181
time:trt 0.178 -0.952 -0.188# In table format
tbl_regression(model_int)| Characteristic | Beta | 95% CI | p-value |
|---|---|---|---|
| time | 0.09 | -0.03, 0.21 | 0.13 |
| trt | 1.2 | -1.0, 3.5 | 0.3 |
| time * trt | 0.03 | -0.04, 0.10 | 0.4 |
| Abbreviation: CI = Confidence Interval | |||
anova(model_int)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
time 2.84453 2.84453 1 200.080 2.3135 0.1298
trt 1.45499 1.45499 1 37.566 1.1834 0.2836
time:trt 0.84596 0.84596 1 200.120 0.6880 0.4078- What this does: Tests if time slope differs by trt (interaction term).
- Expected output: Interaction p<0.05? (trt=2 faster gains); main effects adjusted.
8 Include Random Slopes
Allow subject-specific time slopes.
# LMM with random slope for time
model_rs <- lmer(strength ~ time * trt + (time | id), data = stren_long)
# Summary
summary(model_rs)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: strength ~ time * trt + (time | id)
Data: stren_long
REML criterion at convergence: 818.5
Scaled residuals:
Min 1Q Median 3Q Max
-1.9432 -0.6199 -0.0612 0.5352 3.2512
Random effects:
Groups Name Variance Std.Dev. Corr
id (Intercept) 9.95311 3.1549
time 0.03433 0.1853 -0.03
Residual 0.66469 0.8153
Number of obs: 239, groups: id, 37
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 79.00099 1.74864 34.97864 45.179 <2e-16 ***
time 0.06501 0.11045 33.35236 0.589 0.560
trt 1.13141 1.06373 34.98892 1.064 0.295
time:trt 0.05198 0.06745 33.80024 0.771 0.446
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) time trt
time -0.084
trt -0.954 0.081
time:trt 0.081 -0.953 -0.085# In table format
tbl_regression(model_int)| Characteristic | Beta | 95% CI | p-value |
|---|---|---|---|
| time | 0.09 | -0.03, 0.21 | 0.13 |
| trt | 1.2 | -1.0, 3.5 | 0.3 |
| time * trt | 0.03 | -0.04, 0.10 | 0.4 |
| Abbreviation: CI = Confidence Interval | |||
anova(model_rs)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
time 0.23026 0.23026 1 33.352 0.3464 0.5601
trt 0.75197 0.75197 1 34.989 1.1313 0.2948
time:trt 0.39478 0.39478 1 33.800 0.5939 0.4463# Compute predicted values (including random effects)
stren_long$pred_rs <- predict(model_rs)
# Spaghetti plot of predicted values
ggplot(stren_long, aes(x = time,
y = pred_rs,
color = factor(trt),
group = id)) +
geom_line(alpha = 0.7) +
geom_smooth(aes(group = trt), method = "lm", se = TRUE) + #linear trend
labs(title = "Predicted Strength Trajectories by Treatment Group",
x = "Time (Days)",
y = "Predicted Strength",
color = "Treatment (1=Reps, 2=Weight)") +
theme_minimal()
# Compare models (AIC)
AIC(model_ri, model_int, model_rs) df AIC
model_ri 5 887.1339
model_int 6 893.2104
model_rs 8 834.5395- What this does: Random
(time | id)adds subject-specific slopes; compare fit via AIC (lower better). - Expected output: Random slope SD; better fit (lower AIC) than intercept-only.
9 Model Diagnostics
Check assumptions.
# Residuals plot
plot(model_rs) # Residuals vs fitted
# QQ plot for normality
qqnorm(resid(model_rs))
qqline(resid(model_rs))
# Histogram of residuals
hist(resid(model_rs), main = "Residuals Distribution")
- What this does: Plots check linearity/homoscedasticity; QQ for normality.
- Expected output: Points scatter around line (no patterns); QQ straight; residuals ~normal.