# Test for interaction effect when contrasting specific timepoints in mixed-model

I'm having a mixed model (lmer in R) with 5 repeated measurements (time is categorical), and 6 groups, and several covariates. I'm mainly interested in seeing whether there is a significant difference in 'trajectory' between groups comparing timepoint 3 and 5. Can I somehow test for the significance of the interaction specifically contrasting these two timepoints?

When I use an LRT (anova) to test for the interaction variable, of course I test for the overall significance. From the model summary, I only get the interaction between two contrasts (e.g. group 1-2 at timepoint 3-5). When I use the emmeans package and contrasts, i.e.

PCSmeans <-emmeans(lmerFitPCS, ~ time*group6, adjust = "holm")


I only get the simple contrasts for time or group (e.g. comparing group 1 and 2 at timepoint 3, or comparing group 2 at timepoint 3 and 5). but I would like to get the overall significance of the differences in groups between timepoint 3 and 5.

If I use only a subset of the dataset and use the same model, with only data from timepoints 3 and 5, I can use the LRT test to test the significance I want. However, my model diagnostics (residual plots etc) when I run the model with the subset don't look very good, while the overall model (including all 5 timepoints) does seem to meet the assumptions.

Edit, to make question clear based on @EdM 's comment: I'm looking for a test whether any of the groups differ between those two time points.

What can I do here?

• With the time:group interaction there isn't a unique "overall" measure of the difference between time points 3 and 5, as the outcome values at both time points depend on which group is involved. Do you want an average over the individuals in your study, or an average over the 6 groups (weighting them equally regardless of group size), or a test whether any of the groups differs between those time points? Please edit your question to provide that information, as comments are easy to overlook and can be deleted.
– EdM
May 17, 2022 at 11:59
• What if you add , at = list(time = c(3,5)) to the emmeans call? May 17, 2022 at 13:54
• @RussLenth Thanks for your suggestion. I added the ''at'' term to the emmeans call, but I still only get the simple contrasts (e.g. comparison of group 1 at time point 3 vs 5, or comparison of all groups with each other at one timepoint e.g. group 1-2, group 1-3 etc. at one time point). I'm looking for an overall test testing whether there are any differences between groups across timepoint 3 vs 5. May 17, 2022 at 14:20
• Look at joint_tests May 17, 2022 at 15:15

As I understand it, the goal is to obtain an omnibus test of the interaction of two factors, after excluding some levels of one of them. Here is a parallel example using the warpbreaks data. This uses lm, but the same techniques work for an lmer model.

> require(emmeans)

> warp.lm <- lm(breaks ~ wool * tension, data = warpbreaks)
>
> #---- Estimated marginal means with only two levels of 'tension' included
> (emm <- emmeans(warp.lm, ~ wool * tension,
+                 at = list(tension = c("L", "M"))))
wool tension emmean   SE df lower.CL upper.CL
A    L         44.6 3.65 48     37.2     51.9
B    L         28.2 3.65 48     20.9     35.6
A    M         24.0 3.65 48     16.7     31.3
B    M         28.8 3.65 48     21.4     36.1

Confidence level used: 0.95

> #---- Associated joint tests
> joint_tests(emm)
model term   df1 df2 F.ratio p.value
wool           1  48   2.510  0.1197
tension        1  48   7.519  0.0086
wool:tension   1  48   8.378  0.0057


The above is similar to a type-III ANOVA table, but it applies only to the subsetted factor levels (even though it uses the error estimate from the full model).

It may be that you want these three tests combined into one. That is possible too, using version 1.7.4 or newer of emmeans. We just combine the two factors into one consisting of all combinations of them:

> (emm1 <- comb_facs(emm, c("wool", "tension")))
wool.tension emmean   SE df lower.CL upper.CL
A:L            44.6 3.65 48     37.2     51.9
B:L            28.2 3.65 48     20.9     35.6
A:M            24.0 3.65 48     16.7     31.3
B:M            28.8 3.65 48     21.4     36.1

Confidence level used: 0.95

> joint_tests(emm1)
model term   df1 df2 F.ratio p.value
wool.tension   3  48   6.136  0.0013


However, the OP asks "I'm looking for a test whether any of the groups differ between those two time points." Viewing in this example wool in the role of groups and tension in the role of time, I wonder if what you need is this:

> (diffs <- update(contrast(emm, "consec", simple = "tension"), by = NULL))
contrast wool estimate   SE df t.ratio p.value
M - L    A     -20.556 5.16 48  -3.986  0.0005
M - L    B       0.556 5.16 48   0.108  0.9926

P value adjustment: mvt method for 2 tests

> #---- Joint test of both of these differences
> test(diffs, joint = TRUE)
df1 df2 F.ratio p.value
2  48   7.949  0.0010


This is like combining the tension and wool:tension entries from the first joint_tests output.

Or yet another idea, a multivariate test (Hotelling $$T^2$$) of these results:

> mvcontrast(emm, "consec", mult.name = "wool")
contrast T.square df1 df2 F.ratio p.value
M - L      15.898   2  47   7.783  0.0012



(Hmmm, there is only one test here, so that sidak adjustment really didn't do anything. I need to correct that.)

• Fantastic! I indeed meant the first solution you gave here. Thank you so much, it works. May 31, 2022 at 13:41

If for some reason you can't get the joint_tests() function in emmeans to do what you want, you can do this with a Wald test on multiple parameters. That's based on the variance-covariance matrix of the coefficient estimates (obtained via the vcov() function on the model object) and a matrix describing the combination of Q hypotheses and P parameters that you want to examine together.

In this case you would be doing that on a set of 6 differences, the difference between time_5 and time_3 for each group (for Q = 6 rows of the hypothesis matrix described in the above link). The null hypothesis would be that all of those 6 within-group time_5 - time_3 differences equal 0.

That gives "an overall test testing whether there are any differences between groups across timepoint 3 vs 5." The matrix representing your set of hypotheses has 6 rows (1 hypothesis for each group) and 30+ columns (the number of coefficients for a model with interactions between two categorical predictors, one with 5 levels and one with 6, plus extra coefficients for covariates).

Each row of that matrix then contains a representation of the linear combination of parameter values that you wish to test. For example, consider group $$j$$ at time 5. If the model used standard treatment coding for the categorical predictors, then the outcome estimate $$y_{jt_5}$$ is:*

$$y_{jt_5} = \beta_0 + \beta_j + \beta_{t_5} + \beta_{jt_5},$$

where $$\beta_0$$ is the intercept (for reference group and reference time), $$\beta_j$$ is the difference for group $$j$$ at the reference time, $$\beta_{t_5}$$ is the difference for the reference group at time $$t_5$$, and $$\beta_{jt_5}$$ is the extra difference represented by the interaction coefficient.

For group $$j$$ at time 3 you have:

$$y_{jt_3} = \beta_0 + \beta_j + \beta_{t_3} + \beta_{jt_3},$$

so the difference you seek to evaluate is:

$$y_{jt_5}- y_{jt_3}= (\beta_{t_5}-\beta_{t_3}) + (\beta_{jt_5}-\beta_{jt_3}).$$

So in each of the 6 rows there would be a +1 in the column for $$\beta_{t_5}$$ and a -1 in the column for $$\beta_{t_3}$$. For group $$j$$ (except for the reference group) there would be an additional +1 in the column for $$\beta_{jt_5}$$ and a -1 for $$\beta_{jt_3}$$. All other entries in the Q x P hypothesis matrix would be 0.

The wald.test() function in the R aod package might be adapted to this task.

*This formula assumes that other covariate values are 0. Provided that there are no interactions of covariates with group or time, then any extra contributions to outcome from non-zero covariate values would cancel when you take the difference between time_5 and time_3.

• Moreover, one can use (for the objects in my answer) vcov(emm) or vcov(diffs) to directly estimate the covariate matrix of the EMMs or differences thereof, thus saving having to work through the coefficients. May 19, 2022 at 3:44
• Thank you @EdM for taking your time again to answering me with such detail, it's much appreciated! It worked with the emmeans package but this is still good to know! May 31, 2022 at 13:48