Post hoc test in R for repeated measures ANOVA with 2 within-variables

One option may be to use the afex package to fit the ANOVA and the marginaleffects package to compute contrasts. (Disclaimer: I am the author of marginaleffects.)

First, fit the model:

library(marginaleffects)
library(afex)

df = matrix(nrow = 120, ncol = 0) 
df = as.data.frame(df)
df$ppt_id = c(rep(c(1:30),4))
df$color = c(rep("red", 60), rep("blue", 60))
df$side = c(rep(c(rep("left", 30), rep("right", 30)),2))
df$RT = c(rnorm(30, 600, 50),
          rnorm(30, 650, 50),
          rnorm(30, 700, 50),
          rnorm(30, 600, 50))

mod <- aov_ez(
    id = "ppt_id",
    dv = "RT",
    within = c("side", "color"),
    data = df)

Then, do post hoc processing to get the contrasts. This will tell you how the values of the outcome predicted by the model change when we manipulate the explanators (and their pairwise combinations):

cmp <- comparisons(
    # the model
    mod, 
    # hold other regressors at their means if there are any
    newdata = "mean",
    # vector of variables we want to "manipulate" in the contrast
    variables = c("side", "color"),
    # we care about interactions
    interactions = TRUE)

summary(cmp)
#> Average contrasts 
#>           side      color  Effect Std. Error z value   Pr(>|z|)
#> 1  left - left blue - red 110.489      10.71 10.3148 < 2.22e-16
#> 2 right - left  red - red  51.055      13.84  3.6896 0.00022457
#> 3 right - left blue - red   6.274      12.19  0.5147 0.60674367
#>    2.5 % 97.5 %
#> 1  89.49 131.48
#> 2  23.93  78.18
#> 3 -17.62  30.16
#> 
#> Model type:  afex_aov 
#> Prediction type:  response

There is a detailed vignette on contrasts here: https://vincentarelbundock.github.io/marginaleffects/articles/contrasts.html

Edit: A few notes on interpretation

Start by using the afex package’s default predict() method to compute predicted outcome for an individual with color blue and either values of the side variable:

nd <- data.frame(side = c("left", "right"), color = "blue", ppt_id = 1)
nd
#>    side color ppt_id
#> 1  left  blue      1
#> 2 right  blue      1

predict(mod, newdata = nd)
#>        1        2 
#> 678.7812 609.2096

The difference between left+blue and right+blue is:

diff(predict(mod, newdata = nd))
#>         2 
#> -69.57154

We can obtain the same result with standard errors using the predictions() and comparisons() functions from marginaleffects():

predictions(mod, newdata = nd)
#>   rowid     type predicted std.error statistic p.value conf.low conf.high  side color ppt_id
#> 1     1 response  678.7812  11.02046  61.59284       0 657.1815  700.3808  left  blue      1
#> 2     2 response  609.2096   9.27690  65.66953       0 591.0272  627.3920 right  blue      1

comparisons(
    mod,
    variables = "side",
    newdata = datagrid(color = "blue", ppt_id = 1))
#>   rowid     type        term contrast_side comparison std.error statistic      p.value  conf.low conf.high side color ppt_id
#> 1     1 response interaction  right - left  -69.57154  14.12988 -4.923717 8.491553e-07 -97.26559 -41.87748 left  blue      1

This is a simple “contrast”: What happens to the predicted outcome when side changes from left to right, and color is equal to blue for individual 1?

Now we can compute the same contrast between “left” and “right”, but for individuals in either color conditions:

comparisons(
    mod,
    variables = "side",
    newdata = datagrid(color = c("red", "blue"), ppt_id = 1))
#>   rowid     type        term contrast_side comparison std.error statistic      p.value  conf.low conf.high side color ppt_id
#> 1     1 response interaction  right - left   52.41222  13.74332  3.813650 1.369293e-04  25.47581  79.34864 left   red      1
#> 2     2 response interaction  right - left  -69.57154  14.12988 -4.923717 8.491553e-07 -97.26559 -41.87748 left  blue      1

In the first example I initially gave you, you see “interacted” contrasts: What happens when both the side and color variables change simultaneously?

I’m not an expert in the afex package, so you will have to refer to their documentation to answer the “within” question.

Perhaps I'm not able to understand the answer by @Vincent, or maybe I did not explain clearly what I needed.

Ideally, I would like to get a table with all possible comparisons for this dataset, i.e., there should be six comparisons (because there are 4 possible conditions), something like what this function does:

## combine both independent variables into one
df$condition = paste0(df$side, "_", df$color)
df$condition = as.factor(df$condition)

## calculate pairwise tests
mult_comp = df %>%
  pairwise_t_test(
    RT ~ condition, paired = TRUE, 
    p.adjust.method = "holm")%>%
  select(-statistic, -df)

Now, the problem is that this function can only work with one independent variable (the same is true for a standard function pairwise.t.test), while I actually have two independent variables.

One possible solution is to calculate ANOVA by using the function aov and then use the function TukeyHSD for calculating pairwise comparisons:

anova_df = aov(RT ~ side*color, data = df)

TukeyHSD(anova_df)

The downside is that the calculation is then limited to the Tukey method, which might not always be appropriate. E.g., in this case it would be not statistically correct to use it, since I have a repeated-measures design - and it is not possible to use other corrections with this function.