# 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!

• You seem to be using as the correlation pattern an exchangeable setup where the serial nature of correlations is very unlikely to be modeled well. In most studies, correlation decreases as time spacing increases within a subject, but you are not capturing that. Consider using a serial correlation structure, not necessarily needing random effects. Nov 24, 2021 at 13:48

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. Jun 3, 2020 at 22:09