I'm doing a mixed ANOVA with one within-subjects factor (time) and one between subjects factor (group). The dependent variable is score. When I find an interaction, I want to do simple main effects, i.e. test for an effect of time for each of the groups separately, but using a common error term, as is customary in two-way between ANOVA. I know how to do the latter with rstatix::anova_test() but not the former.
## Reading data and setting up
library(tidyverse)
library(ggpubr)
library(rstatix)
data("anxiety", package = "datarium")
anxiety <- anxiety %>%
gather(key = "time", value = "score", t1, t2, t3) %>%
convert_as_factor(id, time)
## source: https://www.datanovia.com/en/lessons/mixed-anova-in-r/#two-way-mixed
## Two-way mixed ANOVA test
anxiety %>%
anova_test(dv=score,wid=id,
between=group,within=time) %>%
get_anova_table()
## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 group 2 42 4.352 1.90e-02 * 0.168
## 2 time 2 84 394.909 1.91e-43 * 0.179
## 3 group:time 4 84 110.188 1.38e-32 * 0.108
## Effect of time at each level of exercises group
one.way2 <- anxiety %>%
group_by(group) %>%
anova_test(dv = score, wid = id, within = time) %>%
get_anova_table() %>%
adjust_pvalue(method = "bonferroni")
one.way2
## # A tibble: 3 x 9
## group Effect DFn DFd F p `p<.05` ges p.adj
## <fct> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl> <dbl>
## 1 grp1 time 2 28 14.8 4.05e- 5 * 0.024 1.21e- 4
## 2 grp2 time 2 28 77.5 3.88e-12 * 0.086 1.16e-11
## 3 grp3 time 2 28 490. 1.64e-22 * 0.531 4.92e-22
Here, we can see from the DF that it uses a separate error term for each level of the group variable. I understand that it would be more powerful if it uses the entire dataset to calculate the error term.
For reference, here is how to achieve that in a two-way between ANOVA (not mixed):
data("jobsatisfaction", package = "datarium")
## Compute a two-way between-subjects ANOVA and simple main effects
## using a shared (global) error term:
model <- lm(score ~ gender * education_level, data = jobsatisfaction)
jobsatisfaction %>%
group_by(gender) %>%
anova_test(score ~ education_level, error = model)
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## # A tibble: 2 x 8
## gender Effect DFn DFd F p `p<.05` ges
## * <fct> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 male education_level 2 52 132. 3.92e-21 * 0.836
## 2 female education_level 2 52 62.8 1.35e-14 * 0.707
## which is different from the following, where the error term is specific
## to each group (not global)
jobsatisfaction %>%
group_by(gender) %>%
anova_test(score ~ education_level)
## Coefficient covariances computed by hccm()
## Coefficient covariances computed by hccm()
## # A tibble: 2 x 8
## gender Effect DFn DFd F p `p<.05` ges
## * <fct> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <dbl>
## 1 male education_level 2 25 266. 1.39e-17 * 0.955
## 2 female education_level 2 27 42.8 4.19e- 9 * 0.76
So the above achieves a common error term in a two-way between ANOVA; is there a way to achieve the same thing in a two-way mixed ANOVA?
Thank you!
