library(lme4)
library(tidyverse)
library(ggplot2)21 Introduction to Mixed Models
Exercises in R
22 Introduction
Mixed effects models are powerful statistical tools that allow us to analyze data with hierarchical or grouped structures. In these exercises, we will explore the use of mixed models in R.
23 Setup
You will be using the lme4 package for fitting mixed models and you need to install this package. In the examples we will use tidyverse for data manipulation and ggplot2 for visualization, but if you prefer other alternatives, that is of course fine.
24 Random intercept and slope
In these exercises we cover both models with random intercept and random slope using the lme4 package in R.
24.1 Drug effect on blood pressure
A study investigates how a drug affects blood pressure. 200 patients are recruited at 10 hospitals. We suspect that the average blood pressure differs between hospitals, mainly due to different routines for measurement, but we expect the drug effect to be consistent.
## Generate the data
databp <- data.frame(hospital=factor(rep(1:10, each=20)),
drug=rep(c("placebo", "treatment"), times=100),
bp=c(123.9,114.6,119.8,113.3,115.0,121.7,120.2,102.9,121.3,110.4,112.4,111.7,112.6,109.1,114.6,104.3,121.9,113.5,112.1,119.0,121.2,112.6,123.6,118.5,123.2,117.5,121.8,113.8,117.5,112.2,115.6,113.0,112.8,124.9,125.1,108.5,117.1,111.7,123.0,113.7,127.5,121.1,126.0,128.1,125.1,128.8,118.5,124.2,126.9,122.3,128.1,118.7,124.6,116.1,120.9,122.8,128.5,121.5,130.8,131.5,117.8,103.7,125.3,111.7,116.8,120.4,118.9,109.2,121.2,114.6,120.3,117.2,118.4,118.5,119.2,116.9,125.8,117.5,118.7,121.0,125.5,118.3,121.7,112.4,127.3,112.5,131.5,123.2,119.3,110.4,117.0,116.8,119.3,113.8,115.8,115.3,116.6,107.2,118.6,120.1,124.0,124.9,118.8,121.6,129.5,123.4,127.4,118.7,122.6,116.7,127.4,117.1,124.4,120.6,136.1,118.6,128.0,122.3,122.1,121.5,129.1,119.1,122.0,114.7,111.6,122.5,114.5,120.5,131.4,109.6,125.4,115.5,114.0,109.3,113.8,114.2,114.5,120.3,132.3,110.4,118.9,113.8,116.6,104.9,114.3,108.5,117.8,108.1,119.8,108.1,120.2,104.7,108.6,126.1,112.9,111.4,118.1,107.5,117.5,111.8,116.2,112.6,117.1,122.9,113.5,106.8,117.4,113.8,119.4,110.0,111.9,118.6,115.5,107.9,116.1,111.3,122.8,112.7,121.0,109.8,119.3,111.6,118.7,108.7,111.7,123.2,121.2,107.0,115.2,107.3,129.2,119.8,116.9,115.9,116.1,110.8,114.3,110.2,126.5,112.9))Visualize the data using e.g. boxplots or jitter plots, both overall and by hospital.
Solution
## Blood pressure vs drug treatment, treat each measurement as independent
databp %>%
ggplot(aes(x = drug, y = bp)) +
geom_boxplot() +
geom_jitter(width = 0.2, alpha = 0.5) +
theme_bw()
## Blood pressure per hospital and drug treatment
databp |>
ggplot(aes(x = hospital, y = bp, color = drug)) +
geom_boxplot() +
geom_jitter(width = 0.2, alpha = 0.5) +
theme_bw()
Based on the plots, do you think there is variability in blood pressure across hospitals that has to be taken care of?
Tip
Yes, the plots suggest that there is variability in blood pressure across hospitals. The boxplots show different medians for each hospital, indicating that the baseline blood pressure differs between hospitals. This variability should be accounted for in our analysis.
If we ignore any grouping structure, we can fit an ordinary linear regression model.
Fit a linear regression model to the data with blood pressure as the response variable and drug treatment as the independent variable.
Solution
lm_model <- lm(bp ~ drug, data = databp)
summary(lm_model)
Call:
lm(formula = bp ~ drug, data = databp)
Residuals:
Min 1Q Median 3Q Max
-12.329 -4.302 -1.016 4.603 16.271
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 120.4020 0.5835 206.338 < 2e-16 ***
drugtreatment -5.1730 0.8252 -6.269 2.24e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.835 on 198 degrees of freedom
Multiple R-squared: 0.1656, Adjusted R-squared: 0.1614
F-statistic: 39.3 on 1 and 198 DF, p-value: 2.236e-09According to the linear regression model, is the average blood pressure higher in the treatment group compared to the placebo group? By how much does the drug change the blood pressure (on average)?
Tip
The average blood pressure is lower in the treatment group compared to the placebo group. The estimated effect of the drug treatment on blood pressure is approximately -5.17, indicating that patients treated with the drug have an average blood pressure that is 5.17 units lower than those treated with the placebo.
Is this a valid model given the data structure? Why/why not?
Tip
No, not really. The model treats all observations as independent, but patients within the same hospital may be more similar to each other than patients from different hospitals. This violates the independence assumption and can lead to incorrect standard errors and inference.
To account for the variability in blood pressure across hospitals, fit a mixed model with a random intercept for hospital. Random intercepts allow each hospital to have its own baseline blood pressure.
Solution
model1 <- lmer(bp ~ drug + (1 | hospital), data = databp)
summary(model1)In addition to the summary function, try VarCorr and fixef. The function VarCorr gives the variance and correlations in a model, including the the variance of the random intercept for hospital. fixef extracts the estimated fixed effects of the model.
Code
print(VarCorr(model1), comp="Variance") Groups Name Variance
hospital (Intercept) 11.897
Residual 23.234 fixef(model1) (Intercept) drugtreatment
120.402 -5.173 Based on the mixed model results, what is the estimated effect of the drug treatment on blood pressure?
Solution
The estimated effect of the drug treatment on blood pressure is -5.17.
Based on the mixed model, how much variability in blood pressure is due to differences between hospitals?
To estimate the between-hospital variability, we can use the variance component from VarCorr. The proportion of total variance due to differences between hospitals, called the intraclass correlation coefficient (ICC), can be calculated as:
vc <- as.data.frame(VarCorr(model1))
icc <- vc$vcov[1] / (vc$vcov[1] + vc$vcov[2])Based on the variance due to the difference between hospitals, do you think it is appropriate to use an ordinary linear regression model instead of a mixed model? Why or why not?
Tip
The estimated between-hospital variability is substantial, indicating that there are significant differences in blood pressure between hospitals. This suggests that the ordinary linear regression is not appropriate.
The random intercepts for each hospital can be computed using the ranef function. This gives us the deviations from the overall mean (as reported by the fixed effect intercept) for each hospital.
Compute the random effects using ranef. Note that ranef gives the deviations from the overall mean for each hospital.
Code
ranef(model1)$hospital
(Intercept)
1 -2.8246697
2 -0.4104221
3 6.1809293
4 -0.1462214
5 0.2865211
6 4.9829159
7 0.3821799
8 -3.9498001
9 -2.6880141
10 -1.8134188
with conditional variances for "hospital" What would the average blood pressure be at hospital 3 for a patient treated with the placebo?
Solution
fixef(model1)[1] + ranef(model1)$hospital[3,1](Intercept)
126.5829 24.2 Sleep study
The sleepstudy dataset from the lme4 package contains data on reaction times of subjects measured over several days of sleep deprivation (not allowed more than 3h of sleep).
Days 0 and 1 are adaptation and training and will be excluded. Day 2 is baseline, but for simplicity we will rename it to 0.
We will initially work with a subset of the data from only 5 subjects.
##Remove the first two days (before baseline) and adjust Days to be 0 at baseline.
sleepall <- sleepstudy |> filter(Days>=2) |> mutate(Days=Days-2)
##Select a subset of 5 individuals
sleep <- sleepstudy |> filter(Subject %in% c(333, 334, 369, 349, 372), Days>=2) |> mutate(Days=Days-2)Investigate the dataset. What columns are available? Plot the data, e.g. using ggplot and add a linear regression line to visualize the relationship between reaction time and days of sleep deprivation.
Code
head(sleep, 12)
##Plot Reaction time vs days of sleep deprivation, treat all measurements as independent
sleep |> ggplot(aes(x=Days, y=Reaction)) + geom_point() + geom_smooth(method="lm", se=FALSE) + theme_bw() +
labs(x = "Days of sleep deprivation",
y = "Reaction time (ms)")`geom_smooth()` using formula = 'y ~ x'##Plot Reaction time vs days of sleep deprivation, treating each subject as a separate group
sleep |> ggplot(aes(x=Days, y=Reaction, color=Subject)) + geom_point() + geom_smooth(method="lm", se=FALSE) +
theme_bw() + theme(legend.position = "none") +
labs(x = "Days of sleep deprivation",
y = "Reaction time (ms)")`geom_smooth()` using formula = 'y ~ x' Reaction Days Subject
1 276.7693 0 333
2 299.8097 1 333
3 297.1710 2 333
4 338.1665 3 333
5 332.0265 4 333
6 348.8399 5 333
7 333.3600 6 333
8 362.0428 7 333
9 243.3647 0 334
10 254.6723 1 334
11 279.0244 2 334
12 284.1912 3 334
Fit a linear model to the data, ignoring the individual differences.
Code
lm_sleep <- lm(Reaction ~ Days, data = sleep)
summary(lm_sleep)
Call:
lm(formula = Reaction ~ Days, data = sleep)
Residuals:
Min 1Q Median 3Q Max
-39.385 -14.384 0.434 15.952 35.621
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 261.976 5.905 44.36 < 2e-16 ***
Days 14.060 1.412 9.96 3.82e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 20.46 on 38 degrees of freedom
Multiple R-squared: 0.723, Adjusted R-squared: 0.7157
F-statistic: 99.2 on 1 and 38 DF, p-value: 3.816e-12According to this model, how many milliseconds does reaction time increase per day of sleep deprivation?
Solution
The model estimates that reaction time increases by approximately 14.060 milliseconds per day of sleep deprivation.
Fit a mixed model with a random intercept for each subject to account for individual differences in reaction time.
Code
mm1_sleep <- lmer(Reaction ~ Days + (1 | Subject), data = sleep)
summary(mm1_sleep)Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (1 | Subject)
Data: sleep
REML criterion at convergence: 325.3
Scaled residuals:
Min 1Q Median 3Q Max
-1.79871 -0.71049 -0.07413 0.64819 1.91336
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 278.9 16.70
Residual 183.6 13.55
Number of obs: 40, groups: Subject, 5
Fixed effects:
Estimate Std. Error t value
(Intercept) 261.9758 8.4307 31.07
Days 14.0599 0.9351 15.04
Correlation of Fixed Effects:
(Intr)
Days -0.388According to this models, how many milliseconds does reaction time increase per day of sleep deprivation?
Solution
The model estimates that reaction time increases by approximately 14.060 milliseconds per day of sleep deprivation, which is the same as in the linear model.
How much variability in reaction time is explained by individual differences?
Tip
The estimated variance of the random intercept for Subject is 278.9, which indicates substantial between-subject variability in baseline reaction time.
Now, include the full dataset and redo the analyses, i.e. visualize and fit a mixed model with random intercept.
Code
sleepall |> ggplot(aes(x=Days, y=Reaction, color=Subject)) + geom_point() + geom_smooth(method="lm", se=FALSE) +
theme_bw() + theme(legend.position = "none") +
labs(x = "Days of sleep deprivation",
y = "Reaction time (ms)")`geom_smooth()` using formula = 'y ~ x'
mm_sleepall <- lmer(Reaction ~ Days + (1 | Subject), data = sleepall)
summary(mm_sleepall)Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (1 | Subject)
Data: sleepall
REML criterion at convergence: 1430
Scaled residuals:
Min 1Q Median 3Q Max
-3.6261 -0.4450 0.0474 0.5199 4.1378
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 1746.9 41.80
Residual 913.1 30.22
Number of obs: 144, groups: Subject, 18
Fixed effects:
Estimate Std. Error t value
(Intercept) 267.967 10.871 24.65
Days 11.435 1.099 10.40
Correlation of Fixed Effects:
(Intr)
Days -0.354What is the estimated fixed effect of sleep deprivation on reaction time in the full dataset?
What is the estimated variability in reaction time explained by subject differences?
The plot suggests that also the slope varies between subjects, i.e. some subjects are more affected by sleep deprivation than others. Investigate this further by including a random slope for Days within Subject.
Code
mm2_sleepall <- lmer(Reaction ~ Days + (Days | Subject), data = sleepall)
summary(mm2_sleepall)Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (Days | Subject)
Data: sleepall
REML criterion at convergence: 1404.1
Scaled residuals:
Min 1Q Median 3Q Max
-4.0157 -0.3541 0.0069 0.4681 5.0732
Random effects:
Groups Name Variance Std.Dev. Corr
Subject (Intercept) 958.35 30.957
Days 45.78 6.766 0.18
Residual 651.60 25.526
Number of obs: 144, groups: Subject, 18
Fixed effects:
Estimate Std. Error t value
(Intercept) 267.967 8.266 32.418
Days 11.435 1.845 6.197
Correlation of Fixed Effects:
(Intr)
Days -0.062Based on the model with random slope, what is the average effect of sleep deprivation on reaction time? How does it compare to the model with only a random intercept?
Compute the random effects for the model with random slope using ranef(mm2_sleepall). What is the estimated effect of sleep deprivation on reaction time for subject 309?
24.3 Orthodontic measurement over time
The Orthodont dataset from the nlme package contains orthodontic measurements of patients over time. The dataset includes measurements of the distance from the pituitary to the pterygomaxillary fissure (mm).
Use this dataset to study how the distance changes with age, accounting for individual differences in growth patterns. Also, take any relevant covariates into account.
library(nlme)
data("Orthodont", package = "nlme")Start by looking at the data table.
summary(Orthodont) distance age Subject Sex
Min. :16.50 Min. : 8.0 M16 : 4 Male :64
1st Qu.:22.00 1st Qu.: 9.5 M05 : 4 Female:44
Median :23.75 Median :11.0 M02 : 4
Mean :24.02 Mean :11.0 M11 : 4
3rd Qu.:26.00 3rd Qu.:12.5 M07 : 4
Max. :31.50 Max. :14.0 M08 : 4
(Other):84 Visualize the data.
Code
Orthodont |>
ggplot(aes(x = age, y = distance, group = Subject, color = Sex)) +
geom_line() +
geom_point() +
xlab("Age (years)") +
ylab("Distance (mm)") +
theme_bw()
What is the grouping structure in the data?
What covariates could be relevant to include in the model?
First, fit a linear model to the data, ignoring the grouping structure, but include any relevant covariates.
Then, fit a mixed model taking the grouping structure into account. Use a random intercept for Subject to account for individual differences in growth patterns.
Code
The groups would be the individual subjects, and the covariates could be age and Sex.
lmortho <- lm(distance ~ age + Sex, data=Orthodont)
summary(lmortho)
Call:
lm(formula = distance ~ age + Sex, data = Orthodont)
Residuals:
Min 1Q Median 3Q Max
-5.9882 -1.4882 -0.0586 1.1916 5.3711
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.70671 1.11221 15.920 < 2e-16 ***
age 0.66019 0.09776 6.753 8.25e-10 ***
SexFemale -2.32102 0.44489 -5.217 9.20e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.272 on 105 degrees of freedom
Multiple R-squared: 0.4095, Adjusted R-squared: 0.3983
F-statistic: 36.41 on 2 and 105 DF, p-value: 9.726e-13mmortho <- lmer(distance ~ age + Sex + (1 | Subject), data=Orthodont)
summary(mmortho)Linear mixed model fit by REML ['lmerMod']
Formula: distance ~ age + Sex + (1 | Subject)
Data: Orthodont
REML criterion at convergence: 437.5
Scaled residuals:
Min 1Q Median 3Q Max
-3.7489 -0.5503 -0.0252 0.4534 3.6575
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 3.267 1.807
Residual 2.049 1.432
Number of obs: 108, groups: Subject, 27
Fixed effects:
Estimate Std. Error t value
(Intercept) 17.70671 0.83392 21.233
age 0.66019 0.06161 10.716
SexFemale -2.32102 0.76142 -3.048
Correlation of Fixed Effects:
(Intr) age
age -0.813
SexFemale -0.372 0.000What is the estimated effect of age on distance in the linear model? How does it compare to the mixed model?
What is the estimated variability in distance due to differences between subjects? Does this suggest that the random intercept is needed in the model?