Multiple pairwise comparisons for categorical predictors

Dale Barr (@datacmdr) recently had a nice blog post about coding categorical predictors, which reminded me to share my thoughts about multiple pairwise comparisons for categorical predictors in growth curve analysis. As Dale pointed out in his post, the R default is to treat the reference level of a factor as a baseline and to estimate parameters for each of the remaining levels. This will give pairwise comparisons for each of the other levels with the baseline, but not among those other levels. Here's a simple example using the ChickWeight data set (part of the datasets package). As a reminder, this data set is from an experiment on the effect of diet on early growth of chicks. There were 50 chicks, each fed one of 4 diets, and their weights were measured up to 12 times over the first 3 weeks after they were born.
Here are the data:

ggplot(ChickWeight, aes(Time, weight, color = Diet))
 + stat_summary(fun.y = mean, geom = "line")
 + stat_summary(fun.data = mean_se, geom = "pointrange")

Now let's fit the model:

m <- lmer(weight ~ Time * Diet + (1 | Chick),
    data = ChickWeight, REML = F)
print(m, corr = F)

## Linear mixed model fit by maximum likelihood 
## Formula: weight ~ Time * Diet + (1 | Chick) 
##    Data: ChickWeight 
##   AIC  BIC logLik deviance REMLdev
##  5508 5552  -2744     5488    5467
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Chick    (Intercept) 498      22.3    
##  Residual             638      25.3    
## Number of obs: 578, groups: Chick, 50
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   31.508      5.911    5.33
## Time           6.713      0.257   26.09
## Diet2         -2.874     10.191   -0.28
## Diet3        -13.258     10.191   -1.30
## Diet4         -0.398     10.200   -0.04
## Time:Diet2     1.896      0.427    4.44
## Time:Diet3     4.710      0.427   11.04
## Time:Diet4     2.949      0.432    6.82

The Intercept parameter (31.5) is the weight of the "Diet 1" chicks at birth, and the Time parameter (6.7) is their (linear) rate of growth. The Diet2, Diet3, and Diet4 parameters are the birth weights for the chicks in those diet groups relative to Diet 1; the Time:Diet2, Time:Diet3, and Time:Diet4 parameters are the relative slopes -- these are the ones that capture the different effects of diet on rate of growth. But that's rate of growth relative to Diet 1, what if we want to check for differences between Diet 2 and Diet 3?
One option is to re-level the factor so that Diet 2 becomes the reference and the parameters are estimated relative to that:

ChickWeight$DietReleveled <- relevel(ChickWeight$Diet, "2")
m.2 <- lmer(weight ~ Time * DietReleveled + (1 | Chick),
    data = ChickWeight, REML = F)
print(m.2, corr = F)

## Linear mixed model fit by maximum likelihood 
## Formula: weight ~ Time * DietReleveled + (1 | Chick) 
##    Data: ChickWeight 
##   AIC  BIC logLik deviance REMLdev
##  5508 5552  -2744     5488    5467
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Chick    (Intercept) 498      22.3    
##  Residual             638      25.3    
## Number of obs: 578, groups: Chick, 50
## 
## Fixed effects:
##                     Estimate Std. Error t value
## (Intercept)           28.634      8.302    3.45
## Time                   8.609      0.340   25.29
## DietReleveled1         2.874     10.191    0.28
## DietReleveled3       -10.383     11.740   -0.88
## DietReleveled4         2.476     11.748    0.21
## Time:DietReleveled1   -1.896      0.427   -4.44
## Time:DietReleveled3    2.814      0.481    5.84
## Time:DietReleveled4    1.053      0.486    2.17

That worked fine, though it would be kind of annoying and repetitive to have to cycle through all different re-levelings. Especially since we actually have all of the relevant information in the original parameter estimates, we just need to compare them directly to one another. For example, in the original model, the parameter estimates for the growth rate for Diet 2 (Time:Diet2) was 1.896 and for Diet 3 (Time:Diet3) it was 4.71. If we take the difference between those, we get 2.814, which is exactly the parameter estimate for Diet 3 relative to Diet 2 in the releveled model (Time:DietReleveled3). The multcomp package provides a way to do this sort of comparison using just the original model.
First, we need to set up a contrast matrix that defines the comparisons that we want to test. Each column in this matrix corresponds to a parameter estimate from the original model, in the order that they appear in the output. So the first column corresponds to "(Intercept)", the second column to "Time", the sixth column to "Time:Diet2", etc. Each row in the contrast matrix corresponds to a contrast that you want to test and the elements in the matrix are weights for that contrast. The simplest case is when the contrast you want to test corresponds to an estimated parameter: you put a 1 in that column and a 0 in all of the others. For example, to test "Time:Diet1 vs. Time:Diet2" you just put a 1 in the sixth column for the Time:Diet2 parameter. If you want to compare two estimated parameters to one another, you put a -1 in one of their columns and a 1 in the other. Here's the contrast matrix for testing all pairwise slope differences:

contrast.matrix <- rbind(
    `Time:Diet1 vs. Time:Diet2` = c(0, 0, 0, 0, 0, 1, 0, 0),
    `Time:Diet1 vs. Time:Diet3` = c(0, 0, 0, 0, 0, 0, 1, 0),
    `Time:Diet1 vs. Time:Diet4` = c(0, 0, 0, 0, 0, 0, 0, 1),
    `Time:Diet2 vs. Time:Diet3` = c(0, 0, 0, 0, 0, -1, 1, 0),
    `Time:Diet2 vs. Time:Diet4` = c(0, 0, 0, 0, 0, -1, 0, 1),
    `Time:Diet3 vs. Time:Diet4` = c(0, 0, 0, 0, 0, 0, -1, 1))

Now we can use this contrast matrix to test the pairwise comparisons using the glht function from the multcomp package:

library(multcomp)
comps <- glht(m, contrast.matrix)
summary(comps)

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lmer(formula = weight ~ Time * Diet + (1 | Chick),  
##     data = ChickWeight, REML = F)
## 
## Linear Hypotheses:
##                                Estimate Std. Error z value Pr(>|z|)    
## Time:Diet1 vs. Time:Diet2 == 0    1.896      0.427    4.44   <0 ---="" --="" -1.760="" -3.62="" .001="" 0.0016="" 0.001="" 0.01="" 0.05="" 0.1321="" 0.1="" 0.427="" 0.432="" 0.481="" 0.486="" 0="" 1.053="" 11.04="" 1="" 2.17="" 2.814="" 2.949="" 4.710="" 5.84="" 6.82="" codes:="" djusted="" method="" p="" pre="" reported="" signif.="" single-step="" time:diet1="" time:diet2="" time:diet3="" time:diet4="=" values="" vs.="">

This gives us the parameter estimates, standard errors, z-values for each of those comparisons, and adjusted p-values (i.e., corrected for multiple comparisons). If you don't like the default correction, you can change it or suppress it altogether by doing this:

summary(glht(m, contrast.matrix), test = adjusted("none"))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lmer(formula = weight ~ Time * Diet + (1 | Chick),  
##     data = ChickWeight, REML = F)
## 
## Linear Hypotheses:
##                                Estimate Std. Error z value Pr(>|z|)    
## Time:Diet1 vs. Time:Diet2 == 0    1.896      0.427    4.44  8.9e-06 ***
## Time:Diet1 vs. Time:Diet3 == 0    4.710      0.427   11.04  < 2e-16 ***
## Time:Diet1 vs. Time:Diet4 == 0    2.949      0.432    6.82  8.9e-12 ***
## Time:Diet2 vs. Time:Diet3 == 0    2.814      0.481    5.84  5.1e-09 ***
## Time:Diet2 vs. Time:Diet4 == 0    1.053      0.486    2.17  0.03031 *  
## Time:Diet3 vs. Time:Diet4 == 0   -1.760      0.486   -3.62  0.00029 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)

Thanks to Scott Jackson for helping me figure out how to use the contrast matrix.