Mixed model multiple comparisons for interaction between continuous and categorical predictor

Question

I would like to use lme4 to fit a mixed effects regression and multcomp to compute the pairwise comparisons. I have a complex data set with multiple continuous and categorical predictors, but my question can be demonstrated using the built-in ChickWeight data set as an example:

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

Time is continuous and Diet is categorical (4 levels) and there are multiple Chicks per diet. All the chicks started out at about the same weight, but their diets (may) affect their growth rate, so the Diet intercepts should be (more or less) the same, but the slopes might be different. I can get the pairwise comparisons for the intercept effect of Diet like this:

summary(glht(m, linfct=mcp(Diet = "Tukey")))

and, indeed, they are not significantly different, but how can I do the analogous test for the Time:Diet effect? Just putting the interaction term into mcp produces an error:

summary(glht(m, linfct=mcp('Time:Diet' = "Tukey")))
Error in summary(glht(m, linfct = mcp(`Time:Diet` = "Tukey"))) : 
  error in evaluating the argument 'object' in selecting a method for function
 'summary': Error in mcp2matrix(model, linfct = linfct) : 
Variable(s) ‘Time:Diet’ have been specified in ‘linfct’ but cannot be found in ‘model’! 

Answer 1

By default, lmer treats the reference level of a categorical predictor as the baseline and estimates parameters for the other levels. So you get some pairwise comparisons in the default output and you can get the others by using relevel to define a new reference level and re-fitting the model. This has the advantage of letting you use model comparisons or MCMC to get the p-values, but does not correct for multiple comparisons (though you could apply your own correction afterwards).

To use multcomp, you need to define the contrast matrix. Each row in the contrast matrix represents weights for the effects you get in the default model output, starting with the Intercept. So if you want an effect that's already included in the basic output, you just put a "1" in the position corresponding to that effect. Since the parameter estimates are relative to a common reference level, you can get comparisons between any two other levels by setting the weight of one to "-1" and of the other "1". Here's how that would work for the Time:Diet terms in the ChickWeight example:

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))
summary(glht(m, contrast.matrix))

Caveat emptor: This approach seems to use the normal approximation to get p-values, which is somewhat anti-conservative, and then applies some correction for multiple comparisons. The upshot is that this method gives you easy access to as many pairwise parameter estimates and standard errors as you want, but the p-values may or may not be what you want.

(Thanks to Scott Jackson from r-ling-lang-L for help with this)