22  Mixed models II

23 Introduction

These exercises focus on more complex grouped structures with more than one grouping level as well as statistical tests for mixed models.

## Load R packages
library(lme4) # To fit linear mixed models
library(tidyverse) # For data manipulation
library(ggplot2) # For visualization

24 Nested design

24.1 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.

Code

ggplot(nested_data, aes(x = drug, y = weight)) + 
  geom_boxplot() +
  geom_jitter() +
  theme_bw()

Figure 24.1: Weight by Drug Treatment

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()

Figure 24.2: Weight by Litter and Drug Treatment

Fit three different models:

1. A linear regression model ignoring the grouping structure.
2. A linear mixed model with a random intercept for litter nested within line.
3. A linear mixed model with a random intercept and random slope for drug treatment for litters nested within lines.

Code

# 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 line
m2 <- lmer(weight ~ drug + (1 + drug | line/litter_id), 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.

Code

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-15

summary(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.235

summary(m2)

Linear mixed model fit by REML ['lmerMod']
Formula: weight ~ drug + (1 + drug | line/litter_id)
   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_id: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_id:line, 20; line, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept)  25.4883     0.5096   50.02
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_id)
   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

1. Why is model 1 more appropriate than model 0?
2. For the random intercept model 1, does line or litter account for more variation in weight?
3. Does adding a random slope for drug improve the model fit?
4. Based on the best model, what is the estimated fixed effect of the drug treatment? Interpret it.

Solutions

1. Model m1 accounts for the nested structure—ignoring it can underestimate standard errors and inflate Type I error.

2. 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 

3. Compare AIC and p-value from the ANOVA. If model m2 doesn’t improve AIC much or the p-value is large, the added complexity is not justified.

4. From summary(m1), the fixed effect of drugtreated is approximately -2.3, indicating that treatment lowers weight by 2 g on average.

24.2 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':

    toenail

data("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      G

Visualize the data.

Code

eggs |> ggplot(aes(x= Lab, y = Fat, color = Technician, shape=Sample)) +
  geom_jitter() +
  theme_bw()

Fat content by Lab and Technician

What is the grouping structure in the data, and how is it nested? Using this information, fit a linear mixed model.

Tip

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.

There are multiple ways to specify the random effects structure, use the method that is least confusing to you!

Code

##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)

Examine the model output and interpret the random effects. At what level is the most variation in fat content?

25 Model comparisons

25.1 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)

    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              

When assessing the change of distance with increased 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.845

There 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 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 likelihood ratio test. When comparing models with different fixed effects, the comparison should be based on ML rather than REML. Use the anovafunction and remember to set `refit=TRUE as anova will the refit with ML.

What do the results tell you? Which model would you choose?

Code

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 ' ' 1

26 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.

27 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     0

What would the grouping structure be in this case?

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: 4

To 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, does the odds for wheezing increase or decrease with age? What about smoking? Are any of the effects statistically significant?

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.010

We compute the ORs:

exp(fixef(mohio))

(Intercept)         age       smoke 
 0.03425393  0.83797708  1.51403850 

According to this model, do the odds for wheezing increase or decrease with age? What about smoking? What does the model output suggest about the age and smokling effects?