We take off from where we left, and add an interaction term to the mixed model, to run customised linear contrasts. As usual, we are using the sleepstudy dataset:

{r}
# install.packages("lmerTest")
# install.packages("knitr")
# install.packages("multcomp")
# install.packages("dplyr")
# install.packages("ggplot2")
data("sleepstudy")
# Load packages
packs <- c("lmerTest", "knitr", "multcomp", "dplyr", "ggplot2")
lapply(packs, require, character.only = TRUE) 

kable(head(sleepstudy))

We are generating an additional factor (i.e., a categorical variable), lunch, with levels 0 or 1, recording whether the participants had lunch or not before the experiment:

set.seed(88)
lunch <- sample(c(0,1), replace=TRUE, size=18)
sleepstudy$lunch <- factor(lunch)
kable(head(sleepstudy))

At this point, we have 2 levels of variation, days of sleep deprivation and whether the participant had a full belly or not. We are going to investigate whether having lunch or not change the effect of Days on Reaction, by adding an interaction between the factors.

Before fitting the model, we are also changing the format of Days to factor:

# Change the data type of Days from number (continous) to factor (categorical)
sleepstudy$Days <- factor(as.character(sleepstudy$Days), 
                          levels = unique(sleepstudy$Days))
levels(sleepstudy$Days)
# Fit the model
m <- lmer(Reaction ~ Days*lunch - 1 + ( 1 | Subject), data=sleepstudy)
kable(coef(summary(m)))

The above output shows the main effects of Days and lunch and their interactions. We should interpret the main effect of Days without lunch (lunch = 0). In addition to the main effect of lunch, each day will receive an additional coefficient given by the interaction.

The output shows that the main effect of lunch and the interactions are not significant. But we are interested in the differences between all the days and the first day of the experiment, Day 1, so we want to compare all the other days to Day 1. We need again to build a contrasts matrix.

Building the Contrasts Matrix

First, we need to compute all the combinations of the levels of Days and lunch:

group <- paste0(sleepstudy$Days, sleepstudy$lunch)
group <- aggregate(model.matrix(m) ~ group, FUN=mean)
# group
rownames(group) <- group$group
(group <- group[,-1])
kable(head(group)) #displaying only the first 6 rows

The first column specifies the interaction between Days and lunch for each row: we have Day 0 with lunch 0 and 1 (00 and 01), day 1 with lunch 0 and 1 (10 and 11), and so on. The other columns contains zeros and ones specifying the factors combination (1s mark the active level - Days 0 and lunch 1, on row 2) in the extended model matrix.

Now we can specify which comparisons we want to compute (all days vs Day 1):

contrasts <- rbind(group["11",] - group["01",],
                   group["11",] - group["21",],
                   group["11",] - group["31",],
                   group["11",] - group["41",],
                   group["11",] - group["51",],
                   group["11",] - group["61",],
                   group["11",] - group["71",],
                   group["11",] - group["81",],
                   group["11",] - group["91",])
kable(head(contrasts))

Now the matrix specifies couples of contrasts that we want to estimate from the model.

Extracting the contrasts from the model

For extracting the contrasts from the model, we are going to use the glht function from the multcomp package.

# Transform into Matrix
contrast.matrix <- rbind("Lunch Day 0 vs Day 1"=as.numeric(contrasts[1,]),
                         "Lunch Day 2 vs Day 1"=as.numeric(contrasts[3,]),
                         "Lunch Day 3 vs Day 1"=as.numeric(contrasts[4,]),
                         "Lunch Day 4 vs Day 1"=as.numeric(contrasts[5,]),
                         "Lunch Day 5 vs Day 1"=as.numeric(contrasts[6,]),
                         "Lunch Day 6 vs Day 1"=as.numeric(contrasts[7,]),
                         "Lunch Day 7 vs Day 1"=as.numeric(contrasts[8,]),
                         "Lunch Day 8 vs Day 1"=as.numeric(contrasts[8,]),
                         "Lunch Day 9 vs Day 1"=as.numeric(contrasts[8,])
                         )

comparisons <- summary(glht(m, contrast.matrix))

comparisons

From this output, we can extract that the effect of Lunch is critical on Days 4-6, compared to Day 1.

Let’s format the results of the contrasts test into a table:

pq <- comparisons$test

mtests <- data.frame(pq$coefficients, pq$sigma, pq$tstat, pq$pvalues)
error <- attr(pq$pvalues, "error")
pname <- switch(comparisons$alternativ, 
                less = paste("Pr(<", 
                             ifelse(comparisons$df ==0, "z", "t"), ")", 
                             sep = ""), 
                 greater = paste("Pr(>", 
                                 ifelse(comparisons$df == 0, "z", "t"), ")", 
                                 sep = ""), 
                 two.sided = paste("Pr(>|", 
                                   ifelse(comparisons$df == 0, "z", "t"), "|)"
                                   , sep = ""))                                   
colnames(mtests) <- c("Estimate", "Std.Error", 
                      ifelse(comparisons$df ==0, "zvalue", "tvalue"), pname)

mtests <- mtests %>%
  tibble::rownames_to_column("Comparison") %>%
  mutate(`Pr(>|z|)`=ifelse(`Pr(>|z|)`< 0.001, 
                           "< 0.001", 
                           ifelse(`Pr(>|z|)` < 0.01, 
                                  "< 0.01",
                                  ifelse(`Pr(>|z|)` < 0.05, 
                                         "< 0.05",
                                         paste(round(`Pr(>|z|)`,
                                                     4)))))) %>%
  mutate_if(is.numeric, funs(round(.,
                                   digits=2)))
kable(mtests)

This is one way of investigating interactions with mixed models. Remember that, the more comparisons one does, the more there is the chance of incurring in Type I error, and there is the need of correcting the p-values (that the function glht does automatically).