library(lme4)
library(tidyverse)
library(ggplot2)21 Exercises: Mixed Effects Models
Mixed-effects models are powerful statistical tools for analyzing data with hierarchical or grouped structures. In these exercises, we will explore the use of mixed models in R.
You will use the lme4 package for fitting mixed models, so make sure it is installed. In the examples we also use tidyverse and ggplot2 for data manipulation and visualization, but you are welcome to use other alternatives.
21.1 Random intercepts and random slopes
In these exercises we cover both models with random intercepts and random slopes using the lme4 package in R.
Exercise 21.1 (Drug effect on blood pressure) A study investigates how a drug affects blood pressure. A total of 200 patients are recruited from 10 hospitals. We suspect that average blood pressure differs between hospitals, mainly due to different measurement routines, but we expect the drug effect to be consistent across hospitals.
## 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, for example, boxplots or jitter plots, both overall and by hospital. Based on the plots, do you think there is variability in blood pressure across hospitals that should be accounted for?
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()
databp |>
ggplot(aes(x = hospital, y = bp, color = drug)) +
geom_boxplot(outlier.shape = NA) +
geom_jitter(width = 0.2, alpha = 0.5) +
theme_bw()
Yes, the plots suggest that there is variability in blood pressure across hospitals. The hospital-specific boxplots show different baseline levels per hospital, which indicates that the grouping structure should be taken into account.
- If we ignore the grouping structure, we can fit an ordinary linear regression model. Fit a linear regression model with blood pressure as the response and drug treatment as the independent variable. According to this model, is the average blood pressure higher in the treatment group or in the placebo group? By how much does the drug change blood pressure on average?
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-09The average blood pressure is lower in the treatment group than in the placebo group. The estimated treatment effect is approximately -5.17, meaning that the treatment reduces blood pressure by about 5.17 units on average.
- Is this linear model valid given the data structure? Why or why not?
Solution
No, not really. The model treats all observations as independent, but patients from the same hospital may be more similar to each other than patients from different hospitals. This violates the independence assumption and may lead to incorrect standard errors and inference.
- To account for 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. Investigate the model using the functions
summary,VarCorr, andfixef.
Solution
model1 <- lmer(bp ~ drug + (1 | hospital), data = databp)
summary(model1)Linear mixed model fit by REML ['lmerMod']
Formula: bp ~ drug + (1 | hospital)
Data: databp
REML criterion at convergence: 1215.7
Scaled residuals:
Min 1Q Median 3Q Max
-2.36147 -0.65012 -0.08265 0.56163 3.07473
Random effects:
Groups Name Variance Std.Dev.
hospital (Intercept) 11.90 3.449
Residual 23.23 4.820
Number of obs: 200, groups: hospital, 10
Fixed effects:
Estimate Std. Error t value
(Intercept) 120.4020 1.1925 100.968
drugtreatment -5.1730 0.6817 -7.589
Correlation of Fixed Effects:
(Intr)
drugtretmnt -0.286VarCorr(model1) Groups Name Std.Dev.
hospital (Intercept) 3.4492
Residual 4.8202 fixef(model1) (Intercept) drugtreatment
120.402 -5.173 VarCorr() gives the variance components, including the variance of the random intercept for hospital. fixef() extracts the estimated fixed effects.
- Based on the mixed model results, what is the estimated effect of the drug treatment on blood pressure?
Solution
The estimated fixed effect of treatment is about -5.17, so the drug is estimated to lower blood pressure by approximately 5.17 units on average.
- Based on the mixed model, how much variability in blood pressure is due to differences between hospitals?
Solution
We can quantify the between-hospital variability using the variance components and compute the intraclass correlation coefficient (ICC):
vc <- as.data.frame(VarCorr(model1))
vc grp var1 var2 vcov sdcor
1 hospital (Intercept) <NA> 11.89669 3.449158
2 Residual <NA> <NA> 23.23435 4.820202icc <- vc$vcov[1] / (vc$vcov[1] + vc$vcov[2])
icc[1] 0.3386376The ICC gives the proportion of total variance attributable to differences between hospitals.
- Based on the estimated between-hospital variability, do you think it is appropriate to use an ordinary linear regression model instead of a mixed model? Why or why not?
Solution
The estimated between-hospital variability is substantial, and the ICC is clearly above 0. This means that a non-negligible part of the total variation is due to hospital-level differences. Hence, it is more appropriate to use a mixed model, since ordinary linear regression ignores this dependence structure.
- The random intercepts for each hospital can be obtained using the
ranef()function. This gives the deviation of each hospital from the overall intercept.
Compute the random effects using ranef(). What is the estimated average blood pressure at hospital 3 for a patient receiving the placebo?
Solution
## Note that `ranef` gives the deviations from the overall mean for each hospital.
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" ## The average blood pressure at a specific hospital is the sum of the fixed and random effects
fixef(model1)[1] + ranef(model1)$hospital[3,1](Intercept)
126.5829 This gives the estimated mean blood pressure for placebo patients at hospital 3.
Exercise 21.2 (Sleep study) The sleepstudy dataset from the lme4 package contains reaction times for subjects measured over several days of sleep deprivation.
Days 0 and 1 correspond to adaptation and training, and will be excluded. Day 2 is the baseline day, but for simplicity we rename it to 0.
We will initially work with a subset containing 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
ggplot2and add a linear regression line to visualize the relationship between reaction time and days of sleep deprivation, both overall and separately per subject.
Solution
## Print the first 12 rows
head(sleep, 12) 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##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'
- Fit a linear model to the data, ignoring individual differences.
According to this model, by how many milliseconds does reaction time increase per day of sleep deprivation?
Solution
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-12The estimated slope is approximately 14.06, so reaction time increases by about 14.06 milliseconds per day of sleep deprivation.
- Fit a mixed model with a random intercept for each subject to account for individual differences in baseline reaction time.
According to this model, by how many milliseconds does reaction time increase per day of sleep deprivation?
Solution
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.388The fixed effect of Days gives the average increase in reaction time per day. In this subset, it is approximately the same as in the ordinary linear model, around 14.06 ms/day.
- How much variability in reaction time is due to differences between subjects?
Solution
VarCorr(mm1_sleep) Groups Name Std.Dev.
Subject (Intercept) 16.699
Residual 13.551 The random-intercept variance from VarCorr(mm1_sleep) gives the between-subject variability in baseline reaction times. You can also compute an ICC if you want the proportion of total variance attributable to subject differences.
- Now include the full dataset and repeat the analysis: visualize the data and fit a mixed model with a random intercept for subject.
What is the estimated fixed effect of sleep deprivation on reaction time in the full dataset?
What is the estimated variability in reaction time due to subject differences?
Solution
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.354VarCorr(mm_sleepall) Groups Name Std.Dev.
Subject (Intercept) 41.796
Residual 30.217 The fixed effect of Days gives the average increase in reaction time per day in the full dataset. The random-intercept variance gives the variability in baseline reaction time due to differences between subjects.
- The plot suggests that the slope may also vary between subjects, i.e., some subjects are more affected by sleep deprivation than others. Investigate this further by including a random slope for
DayswithinSubject.
Solution
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.062VarCorr(mm2_sleepall) Groups Name Std.Dev. Corr
Subject (Intercept) 30.9573
Days 6.7659 0.178
Residual 25.5264 Based on the model with a random slope, what is the average effect of sleep deprivation on reaction time? How does it compare with the model with only a random intercept?
Compute the random effects for the random-slope model using ranef(mm2_sleepall). What is the estimated effect of sleep deprivation on reaction time for subject 309?
Solution
ranef(mm2_sleepall)$Subject (Intercept) Days
308 24.4992891 8.6020000
309 -59.3723102 -8.1277534
310 -39.4762764 -7.4292365
330 1.3500428 -2.3845976
331 18.4576169 -3.7477340
332 30.5270040 -4.8936899
333 13.3682027 0.2888639
334 -18.1583020 3.8436686
335 -16.9737887 -12.0702333
337 44.5850842 10.1760837
349 -26.6839022 2.1946699
350 -5.9657957 8.1758613
351 -5.5710355 -2.3718494
352 46.6347253 -0.5616377
369 0.9616395 1.7385130
370 -18.5216778 5.6317534
371 -7.3431320 0.2729282
372 17.6826159 0.6623897fixef(mm2_sleepall)["Days"] + ranef(mm2_sleepall)$Subject["309", "Days"] Days
3.307675 The subject-specific slope for subject 309 is the sum of: - the fixed slope for Days, and - the random slope deviation for subject 309.
Exercise 21.3 (Orthodontic measurements over time) The Orthodont dataset from the nlme package contains orthodontic measurements on patients over time. The dataset includes measurements of the distance from the pituitary to the pterygomaxillary fissure (mm).
library(nlme)
data("Orthodont", package = "nlme")Use this dataset to study how the distance changes with age, while accounting for individual differences in growth patterns. Also, take any relevant covariates into account.
- Start by looking at the data. What columns are available?
Solution
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 head(Orthodont)Grouped Data: distance ~ age | Subject
distance age Subject Sex
1 26.0 8 M01 Male
2 25.0 10 M01 Male
3 29.0 12 M01 Male
4 31.0 14 M01 Male
5 21.5 8 M02 Male
6 22.5 10 M02 Male- We can visualize the data useing
ggplot2.
Orthodont |>
ggplot(aes(x = age, y = distance, group = Subject, color = Sex)) +
geom_line() +
geom_point() +
xlab("Age (years)") +
ylab("Distance (mm)") +
theme_bw()
- Based on the plot and what you now know of the available variables, what grouping structure could be relevant to consider?
Solution
The groups would be the individual subjects.
What covariates could be relevant to include in the model?
Solution
The covariates could be age and Sex.
- 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.
Solution
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.000VarCorr(mmortho) Groups Name Std.Dev.
Subject (Intercept) 1.8074
Residual 1.4316 - What 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?
Hint
The estimated effect of age is given by the coefficient of age in each model. In many repeated-measures settings, the estimated slope is similar in the linear and mixed models, but the mixed model gives more appropriate standard errors because it accounts for the dependence within subject.
The estimated variability due to differences between subjects is given by the random-intercept variance in VarCorr(mmortho). If this variance is clearly larger than 0, that suggests that a random intercept is needed.
These exercises focus on more complex grouped structures with more than one grouping level as well as statistical tests for mixed models.
21.2 Nested design
Exercise 21.4 (Mice in litters) We simulate data for mice nested within litters, which are themselves nested within breeding lines. Our aim is to evaluate the effect of a drug on mouse weight while accounting for the hierarchical grouping structure.
## Data simulation
set.seed(42)
# Parameters
n_lines <- 5
n_litters_per_line <- 4
n_mice_per_litter <- 6
# Create data
nested_data <- expand.grid(
line = factor(1:n_lines),
litter = factor(1:n_litters_per_line),
mouse = 1:n_mice_per_litter
) %>%
mutate(
litter_id = factor(paste(line, litter, sep=".")),
drug = sample(c("control", "treated"), n(), replace = TRUE),
line_effect = rnorm(n_lines, 0, 1)[as.integer(line)],
litter_effect = rnorm(n_lines * n_litters_per_line, 0, 0.8)[as.integer(litter_id)],
residual = rnorm(n(), 0, 1),
drug_effect = ifelse(drug == "treated", -2, 0),
weight = 25 + drug_effect + line_effect + litter_effect + residual
)Visualize the data, e.g. using boxplots and jitter plots, both overall and showing the line and litter.
Solution
ggplot(nested_data, aes(x = drug, y = weight)) +
geom_boxplot() +
geom_jitter() +
theme_bw()
ggplot(nested_data, aes(x = litter_id, y = weight, color = drug, shape=line)) +
geom_jitter() +
labs(title = "Weight by Litter", x = "Line.Litter") +
theme_bw()
Fit three different models:
- A linear regression model ignoring the grouping structure.
- A linear mixed model with a random intercept for litter nested within line.
- A linear mixed model with a random intercept and random slope for drug treatment for litters nested within lines.
Solution
# Model 1: Ignore nesting
m0 <- lm(weight ~ drug, data = nested_data)
# Model 2: Random intercept for litter (nested within line)
m1 <- lmer(weight ~ drug + (1 | line/litter_id), data = nested_data)
# Compare to this notation
m1_alt <- lmer(weight ~ drug + (1 | line) + (1 | litter_id:line), data = nested_data)
# or
m1_alt2 <- lmer(weight ~ drug + (1 | line) + (1 | litter_id), data = nested_data)
# Model 3: Add random slope for drug per litter within line
m2 <- lmer(weight ~ drug + (1 + drug | line/litter), data = nested_data)boundary (singular) fit: see help('isSingular')Use summary() to study the different models and anova() to compare the linear mixed models with different random-effects structures.
Solution
summary(m0)
Call:
lm(formula = weight ~ drug, data = nested_data)
Residuals:
Min 1Q Median 3Q Max
-3.1747 -0.8468 -0.0684 0.8763 4.3433
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 25.4423 0.1947 130.65 < 2e-16 ***
drugtreated -2.3272 0.2606 -8.93 6.75e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.418 on 118 degrees of freedom
Multiple R-squared: 0.4032, Adjusted R-squared: 0.3982
F-statistic: 79.74 on 1 and 118 DF, p-value: 6.745e-15summary(m1)Linear mixed model fit by REML ['lmerMod']
Formula: weight ~ drug + (1 | line/litter_id)
Data: nested_data
REML criterion at convergence: 365.4
Scaled residuals:
Min 1Q Median 3Q Max
-2.06105 -0.62779 -0.05877 0.59189 2.72304
Random effects:
Groups Name Variance Std.Dev.
litter_id:line (Intercept) 0.4566 0.6757
line (Intercept) 0.8129 0.9016
Residual 0.9092 0.9535
Number of obs: 120, groups: litter_id:line, 20; line, 5
Fixed effects:
Estimate Std. Error t value
(Intercept) 25.4609 0.4520 56.33
drugtreated -2.3605 0.1904 -12.40
Correlation of Fixed Effects:
(Intr)
drugtreated -0.235summary(m2)Linear mixed model fit by REML ['lmerMod']
Formula: weight ~ drug + (1 + drug | line/litter)
Data: nested_data
REML criterion at convergence: 360
Scaled residuals:
Min 1Q Median 3Q Max
-2.2294 -0.5673 -0.0806 0.6434 2.4779
Random effects:
Groups Name Variance Std.Dev. Corr
litter:line (Intercept) 0.15865 0.3983
drugtreated 0.20968 0.4579 1.00
line (Intercept) 1.17439 1.0837
drugtreated 0.08797 0.2966 -1.00
Residual 0.84284 0.9181
Number of obs: 120, groups: litter:line, 20; line, 5
Fixed effects:
Estimate Std. Error t value
(Intercept) 25.4883 0.5096 50.01
drugtreated -2.3550 0.2471 -9.53
Correlation of Fixed Effects:
(Intr)
drugtreated -0.580
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')anova(m1, m2, refit = FALSE)Data: nested_data
Models:
m1: weight ~ drug + (1 | line/litter_id)
m2: weight ~ drug + (1 + drug | line/litter)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
m1 5 375.44 389.37 -182.72 365.44
m2 9 377.96 403.04 -179.98 359.96 5.4815 4 0.2414- Why is model 1 more appropriate than model 0?
- For the random intercept model 1, does line or litter account for more variation in weight?
- Does adding a random slope for
drugimprove the model fit? - Based on the best model, what is the estimated fixed effect of the drug treatment? Interpret it.
Solution
Model
m1accounts for the nested structure—ignoring it can underestimate standard errors and inflate Type I error.The
VarCorr(m1)output shows variance estimates for line and litter. These tell us how much of the total variance in weight is due to those grouping levels.
# Variance components for model m1
print(VarCorr(m1), comp="Variance") Groups Name Variance
litter_id:line (Intercept) 0.45659
line (Intercept) 0.81293
Residual 0.90917 Compare AIC and p-value from the ANOVA. If model
m2doesn’t improve AIC much or the p-value is large, the added complexity is not justified.From
summary(m1), the fixed effect ofdrugtreatedis approximately-2.3, indicating that treatment lowers weight by 2.3 g on average.
21.3 Lab, technician and sample
We will use the faraway package to analyze a dataset with a hierarchical structure involving laboratory, technician and sample. The dataset eggs contains measurements of fat content in egg powder. Four samples, 2 G samples and 2 H samples, were analyzed by six laboratories, each with two technicians. The samples were each divided into two parts and fat content measured in each. The goal is to assess the effect of lab and technician on the fat content.
# Load the faraway package and the eggs dataset
library(faraway)
Attaching package: 'faraway'The following object is masked from 'package:lme4':
toenaildata("eggs", package = "faraway")
# Display the first few rows of the dataset
head(eggs) Fat Lab Technician Sample
1 0.62 I one G
2 0.55 I one G
3 0.34 I one H
4 0.24 I one H
5 0.80 I two G
6 0.68 I two GVisualize the data.
Solution
eggs |> ggplot(aes(x= Lab, y = Fat, color = Technician, shape=Sample)) +
geom_jitter() +
theme_bw()
What is the grouping structure in the data, and how is it nested? Using this information, fit a linear mixed model.
Solution
Replicate measurements are made on samples, samples are handled by technicians, and technicians are nested within laboratories. Thus, the grouping structure is hierarchical, with technician nested within lab, and sample nested within technician within lab.
##These two models are equivalent, try it!
mfat <- lmer(Fat ~ 1 + (1|Lab/Technician/Sample), data = eggs)
mfat <- lmer(Fat ~ 1 + (1|Lab) + (1|Lab:Technician) + (1|Lab:Technician:Sample), data = eggs)
summary(mfat)Linear mixed model fit by REML ['lmerMod']
Formula:
Fat ~ 1 + (1 | Lab) + (1 | Lab:Technician) + (1 | Lab:Technician:Sample)
Data: eggs
REML criterion at convergence: -64.2
Scaled residuals:
Min 1Q Median 3Q Max
-2.04098 -0.46576 0.00927 0.59713 1.54276
Random effects:
Groups Name Variance Std.Dev.
Lab:Technician:Sample (Intercept) 0.003065 0.05536
Lab:Technician (Intercept) 0.006980 0.08355
Lab (Intercept) 0.005920 0.07694
Residual 0.007196 0.08483
Number of obs: 48, groups:
Lab:Technician:Sample, 24; Lab:Technician, 12; Lab, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.38750 0.04296 9.019VarCorr(mfat) Groups Name Std.Dev.
Lab:Technician:Sample (Intercept) 0.055359
Lab:Technician (Intercept) 0.083548
Lab (Intercept) 0.076941
Residual 0.084828Examine the model output and interpret the random effects. At what level is the most variation in fat content?
Hint
Use VarCorr(mfat) to compare the variance components at the lab, technician-within-lab, and sample-within-technician-within-lab levels. The level with the largest estimated variance contributes the most to the total variability in fat content.
21.4 Model comparisons
Exercise 21.5 (Orthodont) Let’s get back to the Orthodont dataset in the nlme package, which contains measurements of 27 children’s teeth growth over time.
library(nlme)
data("Orthodont", package = "nlme")
summary(Orthodont)When assessing the change in distance with increasing age, it is important to take the grouping structure into account. Each child is measured multiple times, so we need to account for the correlation of measurements within subjects. Add a random intercept for subject and fit a linear mixed model to assess the association between distance and age.
mo1 <- lmer(distance ~ age + (1|Subject), data = Orthodont)
summary(mo1)Linear mixed model fit by REML ['lmerMod']
Formula: distance ~ age + (1 | Subject)
Data: Orthodont
REML criterion at convergence: 447
Scaled residuals:
Min 1Q Median 3Q Max
-3.6645 -0.5351 -0.0129 0.4874 3.7218
Random effects:
Groups Name Variance Std.Dev.
Subject (Intercept) 4.472 2.115
Residual 2.049 1.432
Number of obs: 108, groups: Subject, 27
Fixed effects:
Estimate Std. Error t value
(Intercept) 16.76111 0.80240 20.89
age 0.66019 0.06161 10.72
Correlation of Fixed Effects:
(Intr)
age -0.845There are both boys and girls in the dataset. Your PI wants to estimate the effect of sex on the distance. How could you do this? With a random intercept per subject, build models with different sets of fixed effects. Start with only intercept and age, then add also sex and finally also the interaction between age and sex.
mof2 <- lmer(distance ~ age + Sex + (1|Subject), data = Orthodont)
mof3 <- lmer(distance ~ age*Sex + (1|Subject), data = Orthodont)Compare the models using a likelihood ratio test. When comparing models with different fixed effects, the comparison should be based on ML rather than REML. Use the anova function and remember to set refit = TRUE, as anova will refit the models with ML.
What do the results tell you? Which model would you choose?
Solution
anova(mo1, mof2, mof3, refit = TRUE)refitting model(s) with ML (instead of REML)Data: Orthodont
Models:
mo1: distance ~ age + (1 | Subject)
mof2: distance ~ age + Sex + (1 | Subject)
mof3: distance ~ age * Sex + (1 | Subject)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
mo1 4 451.39 462.12 -221.69 443.39
mof2 5 444.86 458.27 -217.43 434.86 8.5331 1 0.003488 **
mof3 6 440.64 456.73 -214.32 428.64 6.2174 1 0.012650 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Choose the simplest model that adequately describes the data. If adding Sex improves the fit, keep mof2. If the interaction age:Sex also improves the fit, keep mof3. Otherwise, prefer the simpler model.
21.5 Computing p-values for mixed models
In this section you will use the lmerTest package to compute p-values for mixed models using Satterthwaite’s approximation. This package extends the lme4 package and provides functions to compute p-values for fixed effects in linear mixed models.
Go back to the exercises “Drug effect on blood pressure” and “Sleep study” and use the lmerTest package to compute p-values for the fixed effects in the models you fitted.
Solution
library(lmerTest)
Attaching package: 'lmerTest'The following object is masked from 'package:lme4':
lmerThe following object is masked from 'package:stats':
stepmodel1_test <- lmer(bp ~ drug + (1 | hospital), data = databp)
summary(model1_test)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: bp ~ drug + (1 | hospital)
Data: databp
REML criterion at convergence: 1215.7
Scaled residuals:
Min 1Q Median 3Q Max
-2.36147 -0.65012 -0.08265 0.56163 3.07473
Random effects:
Groups Name Variance Std.Dev.
hospital (Intercept) 11.90 3.449
Residual 23.23 4.820
Number of obs: 200, groups: hospital, 10
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 120.4020 1.1925 10.6685 100.968 < 2e-16 ***
drugtreatment -5.1730 0.6817 189.0000 -7.589 1.43e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
drugtretmnt -0.286mm1_sleep_test <- lmer(Reaction ~ Days + (1 | Subject), data = sleep)
summary(mm1_sleep_test)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
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 df t value Pr(>|t|)
(Intercept) 261.9758 8.4307 5.5251 31.07 2.04e-07 ***
Days 14.0599 0.9351 34.0000 15.04 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
Days -0.388mm_sleepall_test <- lmer(Reaction ~ Days + (1 | Subject), data = sleepall)
summary(mm_sleepall_test)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
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 df t value Pr(>|t|)
(Intercept) 267.967 10.871 22.152 24.65 <2e-16 ***
Days 11.435 1.099 125.000 10.40 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
Days -0.354The summary() output from lmerTest includes approximate p-values for the fixed effects.
21.6 Logistic mixed model
Here we will use the ohio dataset in the faraway package. This dataset contains information on 537 children and their wheezing status, a binary variable indicating if the child has a pulmonary problem or not. Each child has multiple observations (one for each year from ages 7 to 10, where age 9 is coded as 0). There is also an indicator variable for whether the mother of the child is a smoker.
Take a look at the first two children:
library(faraway)
data("ohio", package = "faraway")
head(ohio, 8) resp id age smoke
1 0 0 -2 0
2 0 0 -1 0
3 0 0 0 0
4 0 0 1 0
5 0 1 -2 0
6 0 1 -1 0
7 0 1 0 0
8 0 1 1 0What would the grouping structure be in this case?
Solution
The grouping variable is id, since each child has repeated observations over time.
If we ignore the grouping structure, we can fit a logistic regression model using glm():
mohio_glm <- glm(resp ~ age + smoke, family = binomial, data = ohio)
summary(mohio_glm)
Call:
glm(formula = resp ~ age + smoke, family = binomial, data = ohio)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.88373 0.08384 -22.467 <2e-16 ***
age -0.11341 0.05408 -2.097 0.0360 *
smoke 0.27214 0.12347 2.204 0.0275 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1829.1 on 2147 degrees of freedom
Residual deviance: 1819.9 on 2145 degrees of freedom
AIC: 1825.9
Number of Fisher Scoring iterations: 4To compute the odds ratios (OR), we exponentiate the coefficients:
exp(coef(mohio_glm))(Intercept) age smoke
0.1520213 0.8927821 1.3127689 According to this model, do the odds of wheezing increase or decrease with age? What about smoking? Are any of the effects statistically significant?
Hint
The sign of the coefficient for age tells us whether the odds of wheezing increase or decrease with age. The coefficient for smoke tells us how maternal smoking is associated with the odds of wheezing. Statistical significance can be assessed from the p-values in summary(mohio_glm), and the odds ratios are obtained from exp(coef(mohio_glm)).
The above model is however incorrect as it ignores the fact that each child has multiple observations. To take this grouping into account, we can use a mixed model with a random intercept for each child.
mohio <- glmer(resp ~ age + smoke + (1|id), family=binomial, data=ohio)
summary(mohio)Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: resp ~ age + smoke + (1 | id)
Data: ohio
AIC BIC logLik -2*log(L) df.resid
1597.9 1620.6 -794.9 1589.9 2144
Scaled residuals:
Min 1Q Median 3Q Max
-1.4027 -0.1802 -0.1577 -0.1321 2.5176
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 5.49 2.343
Number of obs: 2148, groups: id, 537
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.37395 0.27496 -12.271 <2e-16 ***
age -0.17676 0.06797 -2.601 0.0093 **
smoke 0.41478 0.28705 1.445 0.1485
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) age
age 0.227
smoke -0.419 -0.010We compute the ORs:
exp(fixef(mohio))(Intercept) age smoke
0.03425393 0.83797708 1.51403850 According to this model, do the odds of wheezing increase or decrease with age? What about smoking? What does the model output suggest about the age and smoking effects?
Hint
The sign of the fixed effects tells us whether the odds of wheezing increase or decrease with age and with maternal smoking. The exponentiated fixed effects, exp(fixef(mohio)), give the corresponding odds ratios. Comparing this model with the ordinary logistic regression allows us to see how accounting for repeated observations within child affects the estimates and inference.