How can I run post-hoc tests for a significant interaction in Two-Way Mixed ANOVA with the correct (within subjects) error term in R?

Question

I am attempting run a Fisher's LSD post hoc test on a Two-Way Mixed Model ANOVA using the "afex" and "emmeans" packages. The data I am using has one between-subjects factor "group" which has 2 levels, and one within-subjects factor "time" which has 3 levels (i.e. its a 2 x 3 design). The DV is "score". Here is some model data that I created to replicate the error:

library(tidyverse)
library(dplyr)
library(afex)
library(emmeans)
# Create data
set.seed(35)
df1 <- data.frame(id = factor(rep(1:12, each = 3)),
                  time = factor(rep(c(1:3), 12)),
                  group = factor(rep(1:2, each = 18)),
                  score = rnorm(36, 20, 5))

...And here is how I am constructing the ANOVA model:

# Run Two-way Mixed Model ANOVA
model1 <- aov_car(score ~ time*group +
                            Error(id/time), 
                          data = df1)

The output F-table for this two-way ANOVA looks like this, though, when I run it on my actual data there is a significant interaction effect:

Anova Table (Type 3 tests)

Response: score
      Effect          df   MSE      F  ges p.value
1      group       1, 10 27.50 4.11 + .130    .070
2       time 1.72, 17.19 28.02   0.50 .031    .589
3 group:time 1.72, 17.19 28.02   1.09 .065    .350
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

Because there is a significant interaction in my actual data, I would like to conduct follow-up tests. Specifically, I would like to look at the effect of time for each group separately (in a within-subjects manner). I tried to run both a simple main effects test and multiple comparisons(Fisher's LSD) using this code:

# Run whithin-subject simple main effect  
simple.effect <-emmeans(model1, ~time|group, model = 
                "multivariate")
test(pairs(simple.effect), joint=TRUE)

# Run Fisher's LSD multiple comparison 
pairs(simple.effect, adjust='none')

Which returns this:

> # Run whithin-subject simple main effect  
> simple.effect <-emmeans(model1, ~time|group, model = "multivariate")
> test(pairs(simple.effect), joint=TRUE)
 group df1 df2 F.ratio p.value note
 1       2  10   1.009  0.3989  d  
 2       2  10   0.374  0.6969  d  

d: df1 reduced due to linear dependence 
> # Run Fisher's LSD multiple comparison 
> pairs(simple.effect, adjust='none')
group = 1:
 contrast estimate   SE df t.ratio p.value
 X1 - X2    -3.792 2.74 10  -1.386  0.1959
 X1 - X3     0.231 2.35 10   0.098  0.9236
 X2 - X3     4.023 3.33 10   1.209  0.2544

group = 2:
 contrast estimate   SE df t.ratio p.value
 X1 - X2    -0.199 2.74 10  -0.073  0.9435
 X1 - X3    -2.030 2.35 10  -0.863  0.4082
 X2 - X3    -1.831 3.33 10  -0.550  0.5942

For both of these sets of tests the df error reads 10 (which is the error associated with "group" - the between-subjects factor). I believe that the df error should actually be 17.19 (which would correspond to the error term for the within-subjects factor ("time").

It seems to me as though "emmeans" is using the wrong df/MS error to run the follow-up tests, but I don't know how to tell it to use the right one. Any help would be greatly appreciated.

Answer 1

This is what I get:

> simple.effect <- emmeans(model1, ~time|group, model = "univariate")
aov object missing, substituting multivariate/lm model.
to get univariate tests, refit ANOVA with include_aov = TRUE

Note from the message that asking for the univariate model did not produce it. Doing what is recommended in the message, I re-fitted the model:

> model2 <- aov_car(score ~ time*group +
+                       Error(id/time), include_aov = TRUE, 
+                   data = df1)
Contrasts set to contr.sum for the following variables: group

> simple.effect2 <- emmeans(model2, ~time|group, model = "univariate")
> pairs(simple.effect2)
group = 1:
 contrast estimate   SE df t.ratio p.value
 X1 - X2    -3.792 2.83 20  -1.338  0.3912
 X1 - X3     0.231 2.83 20   0.082  0.9963
 X2 - X3     4.023 2.83 20   1.420  0.3499

group = 2:
 contrast estimate   SE df t.ratio p.value
 X1 - X2    -0.199 2.83 20  -0.070  0.9973
 X1 - X3    -2.030 2.83 20  -0.716  0.7568
 X2 - X3    -1.831 2.83 20  -0.646  0.7966

P value adjustment: tukey method for comparing a family of 3 estimates

The pairwise differences each show 20 d.f. Also (not shown), the means themselves have 29.9 d.f. (Note that you have to specify the univariate model because it still defaults to the multivariate one).

By the way, I don't think most users would deem an interaction "significant" when it has a P value of 0.35. They would just obtain the EMMs for time without conditioning on group. The d.f. remain at 10 with the multivariate model and at 20 with the univariate model.