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

Question

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")
pairs(PCSmeans, simple = "each", 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?

Answer 1

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

P value adjustment: sidak 

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

Answer 2

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.