MANOVA effect size (partial eta squared) in R

For ANOVA, one can easily get the partial eta squared (np2) effect size with effectsize::eta_squared:

> model <- aov(mpg ~ factor(cyl), data = mtcars)
> effectsize::eta_squared(model)
For one-way between subjects designs, partial eta squared is equivalent to eta squared.
Returning eta squared.
# Effect Size for ANOVA

Parameter   | Eta2 |       90% CI
---------------------------------
factor(cyl) | 0.73 | [0.57, 0.82]

Furthermore, while conducting MANOVA in SPSS, one can get np2 easily for each DV by simply checking the "Estimates of effect size" box.

However, trying to get the np2 for each DV in R:

> model <- manova(cbind(mpg, hp) ~ factor(cyl), data = mtcars)
> effectsize::eta_squared(model)
# Effect Size for ANOVA (Type I)

Parameter   | Eta2 (partial) |       90% CI
-------------------------------------------
factor(cyl) |           0.46 | [0.28, 0.57]

Or yet:

> heplots::etasq(model)
                eta^2
factor(cyl) 0.4564352

Only outputs a single np2 (for the global model?), rather than a separate effect size for each DV.

Documentation for ?eta_squared notes:

for mlm / maov models, effect sizes are computed for each response separately

But this is obviously not the case here. Yet, when I check class of object, it includes class maov:

> class(model)
[1] "manova" "maov"   "aov"    "mlm"    "lm"    

One could of course rerun an anova for each effect size by defining separate models:

> model <- aov(mpg ~ factor(cyl), data = mtcars)
> effectsize::eta_squared(model)
Parameter   | Eta2 |       90% CI
---------------------------------
factor(cyl) | 0.73 | [0.57, 0.82]

> model <- aov(hp ~ factor(cyl), data = mtcars)
> effectsize::eta_squared(model)
Parameter   | Eta2 |       90% CI
---------------------------------
factor(cyl) | 0.71 | [0.55, 0.80]

Or write a lapply call for each model:

aov.list <- list(model1 <- aov(mpg ~ factor(cyl), data = mtcars),
                 model2 <- aov(hp ~ factor(cyl), data = mtcars))

lapply(aov.list, effectsize::eta_squared)

[[1]]
Parameter   | Eta2 |       90% CI
---------------------------------
factor(cyl) | 0.73 | [0.57, 0.82]

[[2]]    
Parameter   | Eta2 |       90% CI
---------------------------------
factor(cyl) | 0.71 | [0.55, 0.80]

But it doesn't seem to be the most efficient. Surely, there is a way to do this from effectsize::eta_squared?

Question

How to get the partial eta squared for MANOVA in R for each DV the correct way?

You can fit the multivariate analysis of variance model using lm instead. You then have the option to access partial effect sizes for individual variables or globally.

options(contrasts = c("contr.sum", "contr.poly"))
model <- lm(cbind(mpg, hp) ~ factor(cyl), 
            data = mtcars)
# Compute effect size variable by variable
effectsize::eta_squared(model, partial = FALSE)
#> # Effect Size for ANOVA (Type I)
#> 
#> Response |   Parameter | Eta2 |       95% CI
#> --------------------------------------------
#> mpg      | factor(cyl) | 0.73 | [0.57, 1.00]
#> hp       | factor(cyl) | 0.71 | [0.55, 1.00]

# Compute MANOVA statistic and global effect size
mtest <- car::Anova(model)
effectsize::eta_squared(mtest)
#> # Effect Size for ANOVA (Type II)
#> 
#> Parameter   | Eta2 (partial) |       95% CI
#> -------------------------------------------
#> factor(cyl) |           0.46 | [0.28, 1.00]
#> 
#> - One-sided CIs: upper bound fixed at [1.00].

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