# post hoc test on multivariate repeated measures linear mixed effects model (lmer)

I'm running a multivariate linear mixed effects analysis on data with the following structure:

Where: Participants are nested within group, Time is a pre-post-post repeated measurement. The repeated measurement of Time is both a fixed effect and a random slope across the group/participant random effects nesting structure.

The variable consists of the three stacked DV's (i.e., c(DV1,DV2, DV3) for the analysis as per this post

The equation I created to analyze this repeated measures linear mixed-effects model is as follows:

model <- lmer(value ~ variable + variable:Time - 1 + (0 + Time | Group_ID/Participant_ID)


The output of the analysis was:

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: value ~ variable + variable:Time - 1 + (0 + Time | Group_ID/Participant_ID)
Data: Fixed_Data
Control: lmerControl(optCtrl = list(maxfun = 20000))

REML criterion at convergence: 12152

Scaled residuals:
Min      1Q  Median      3Q     Max
-6.7460 -0.3175 -0.0032  0.3513  4.5339

Random effects:
Groups                            Name  Variance Std.Dev. Corr
Participant_ID:Group_ID           Time1 1.6447   1.2825
Time2 0.7142   0.8451   1.00
Time3 1.6477   1.2836   1.00 1.00
Group_ID                          Time1 2.0583   1.4347
Time2 0.8759   0.9359   0.61
Time3 0.6493   0.8058   0.60 1.00
Residual                                7.0815   2.6611
Number of obs: 2451, groups:  Participant_ID:Group_ID, 280; Group_ID, 62

Fixed effects:
Estimate Std. Error       df t value Pr(>|t|)
DV1                        5.8524     0.2549 112.1553  22.962  < 2e-16 ***
DV2                        5.8836     0.2549 112.1553  23.084  < 2e-16 ***
DV3                        26.0502    0.2549 112.1553 102.209  < 2e-16 ***
DV1:Time2                  0.6434     0.2698 212.5869   2.385   0.0180 *
DV2:Time2                  0.6351     0.2698 212.5869   2.354   0.0195 *
DV3:Time2                  2.5075     0.2698 212.5869   9.293  < 2e-16 ***
DV1:Time3                  0.3340     0.2735 217.3497   1.221   0.2233
DV2:Time3                  0.2661     0.2735 217.3497   0.973   0.3318
DV3:Time3                  1.3563     0.2735 217.3497   4.959 1.43e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


The part I'm trying to figure out now is how to run a post-hoc test to determine the differences between the different time points for the fixed effects i.e. :

                            Estimate Std. Error z value Pr(>|z|)
DV1:Time2 - DV1 == 0        --          --       --       --
DV1:Time3 - DV1:Time2 == 0  --          --       --       --
DV2:Time2 - DV2 == 0        --          --       --       --
DV2:Time3 - DV2:Time2 == 0  --          --       --       --
DV3:Time2 - DV3 == 0        --          --       --       --
DV3:Time3 - DV3:Time2 == 0  --          --       --       --


When I attempt this procedure though multicomp::glht, or other comparable post hoc tests, I get an output that contains every possible combination of the multivariate space e.g.:

                            Estimate Std. Error z value Pr(>|z|)
DV1:Time2 - DV1 == 0        --          --       --       --
DV1:Time3 - DV1:Time2 == 0  --          --       --       --
DV1:Time2 - DV2 == 0        --          --       --       --
DV1:Time3 - DV2:Time2 == 0  --          --       --       --
etc...


Would it be inappropriate to do planned comparison t-tests by hand? Or is there a specific way in which I can running the post-hoc test that is causing all these comparisons to run?

Thanks!

I was able to figure this out via this post.

Essentially, we need to look at the contrast matrix by calling the emmeans contrasts:

emm_model<-emmeans(model, ~ variable + Time)


Then, look at the contrast matrix:

coef(pairs(emm_model))


Decide on which contrasts to use, then form then save your needed contrasts:

(contr_mat <- coef(pairs(emm_model))[, c("c.3", "c.6", "c.11", "c.14", "c.18", "c.21")])


Finally, run the planned comparisons

emmeans(model, ~ variable + Time, contr = contr_mat, adjust = "holm")

• OK, but sometimes it helps to read the documentation. It seems to me that using by = "Time" gets you a long ways there. There are also methods to subset and combine emmGrid objects. – Russ Lenth Jun 3 '20 at 22:09