anova_stats returns wrong partial eta squared values for repeated measures anova

[frankpapenmeier]

I noticed that the anova_stats function returns wrong partial eta squared values for repeated measures ANOVAs.

Please see the following reproducible example for details regarding the issue and its source:

library(sjstats)

# generate some data

set.seed(325)

dat <- expand.grid(
  factor1 = c("A","B"),
  factor2 = c("X","Y"),
  participant = 1:10
)
dat$participant <- factor(dat$participant)

dat$value <- rnorm(length(dat$participant))

# add some main effects
dat$value[dat$factor1 == "A"] <- dat$value[dat$factor1 == "A"] + rnorm(dat$value[dat$factor1 == "A"], 5)
dat$value[dat$factor2 == "X"] <- dat$value[dat$factor2 == "X"] + rnorm(dat$value[dat$factor2 == "X"], 2)


# 2x2 repeated measures ANOVA
model <- aov(value ~ factor1*factor2 + Error(participant/(factor1*factor2)), dat)
summary(model)

# Error: participant
# Df Sum Sq Mean Sq F value Pr(>F)
# Residuals  9  17.41   1.935               
# 
# Error: participant:factor1
# Df Sum Sq Mean Sq F value   Pr(>F)    
# factor1    1  317.9   317.9   225.6 1.12e-07 ***
#   Residuals  9   12.7     1.4                     
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Error: participant:factor2
# Df Sum Sq Mean Sq F value  Pr(>F)   
# factor2    1  28.70   28.70   12.16 0.00686 **
#   Residuals  9  21.24    2.36                   
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Error: participant:factor1:factor2
# Df Sum Sq Mean Sq F value Pr(>F)  
# factor1:factor2  1  4.814   4.814   4.581  0.061 .
# Residuals        9  9.457   1.051                 
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Thus, partial eta squared values should be:
#
# main effect factor 1: 317.9 / (317.9 + 12.7) = 0.961585
# main effect factor 2: 28.70 / (28.70 + 21.24) = 0.5746896
# interaction effect factor 1:factor2: 4.814 / (4.814 + 9.457) = 0.3373274
#

# This is the results of anova_stats from the sjstats package
anova_stats(model)

# main effect factor 1: 0.971
# main effect factor 2: 0.752
# interaction effect factor 1:factor2: 0.337

# --> partial eta squared of last effect is correct, but of the other effects is wrong

# the source of this issue lies in the following line 160 of the anova_stats.R script:
#
# "ss.resid <- aov.sum[["sumsq"]][nrow(aov.sum)]"

# --> this line assumes that the last row in the anova structure contains the residuals

# --> thus it calculates the partial eta squared effect size for all effects with the residuals of the last effect in the model

# So here is a reproduction of the false partial eta squared values as calculated by the sjstats package for the above example:

# false ss.resid (9.457) for all effects:
# main effect factor 1: 317.9 / (317.9 + 9.457) = 0.9711111 -> false result reported by anova_stats
# main effect factor 2: 28.70 / (28.70 + 9.457) = 0.7521556 -> false result reported by anova_stats
# interaction effect factor 1:factor2: 4.814 / (4.814 + 9.457) = 0.3373274 -> this one is correct

[strengejacke]

The development of effectsizes from anova models is better supported in the effectsize package, which is also used by the parameters package. Example:

library(parameters)

# generate some data

set.seed(325)

dat <- expand.grid(
  factor1 = c("A","B"),
  factor2 = c("X","Y"),
  participant = 1:10
)
dat$participant <- factor(dat$participant)

dat$value <- rnorm(length(dat$participant))

# add some main effects
dat$value[dat$factor1 == "A"] <- dat$value[dat$factor1 == "A"] + rnorm(dat$value[dat$factor1 == "A"], 5)
dat$value[dat$factor2 == "X"] <- dat$value[dat$factor2 == "X"] + rnorm(dat$value[dat$factor2 == "X"], 2)


# 2x2 repeated measures ANOVA
model <- aov(value ~ factor1*factor2 + Error(participant/(factor1*factor2)), dat)


model_parameters(model, eta_squared = "partial")
#> # participant
#> 
#> Parameter | Sum_Squares | df | Mean_Square
#> ------------------------------------------
#> Residuals |       17.41 |  9 |        1.93
#> 
#> # participant:factor1
#> 
#> Parameter | Sum_Squares | df | Mean_Square |      F |      p | Eta2 (partial)
#> -----------------------------------------------------------------------------
#> factor1   |      317.89 |  1 |      317.89 | 225.57 | < .001 |           0.96
#> Residuals |       12.68 |  9 |        1.41 |        |        |               
#> 
#> # participant:factor2
#> 
#> Parameter | Sum_Squares | df | Mean_Square |     F |     p | Eta2 (partial)
#> ---------------------------------------------------------------------------
#> factor2   |       28.70 |  1 |       28.70 | 12.16 | 0.007 |           0.57
#> Residuals |       21.24 |  9 |        2.36 |       |       |               
#> 
#> # participant:factor1:factor2
#> 
#> Parameter       | Sum_Squares | df | Mean_Square |    F |     p | Eta2 (partial)
#> --------------------------------------------------------------------------------
#> factor1:factor2 |        4.81 |  1 |        4.81 | 4.58 | 0.061 |           0.34
#> Residuals       |        9.46 |  9 |        1.05 |      |       |               
#> 
#> Anova Table (Type 1 tests)

^{Created on 2022-01-21 by the reprex package (v2.0.1)}

I suggest using either effectsize or parameters. The anova-functions in sjstats are going to be deprecated.