Unbiased estimator of Cohen's f

[Daiki-Nakamura-git]

Hello, thank you for your continued development of effectsize package. I have a suggestion on how to calculate Cohen's f.

In the current code, Cohen's f appears to be calculated by transforming $\eta^2$. However, $\eta^2$ is a biased estimator and the Cohen's f derived from it is also biased.

The unbiased estimator of Cohen's f has been found to be derived by the following equation (Grissom & Kim, 2012, p.181, eq6.7).

How about modifying the code to calculate the above unbiased estimator? The following simulation shows that the above estimator is unbiased.

muvec <- c(3,4,5) # mu vector
sigma <- 1　# common population variance across groups
k <- 3      # number of groups
n <- 50     # sample size for each group
N <- n*k    # total sample size
R <- 10000  # number of repetitions

sigma_b_square <- sum((muvec-mean(muvec))^2)/k
f <- sqrt(sigma_b_square)/sigma # true Cohen's f

f_biased <- NULL
f_unbiased <- NULL
eta2_biased <- NULL

library(parameters)
set.seed(123)
for(i in 1:R){
  A <- rnorm(n, muvec[1], sigma)
  B <- rnorm(n, muvec[2], sigma)
  C <- rnorm(n, muvec[3], sigma)
  group <- rep(c("A","B","C"), each=n)
  dat <- data.frame(group=group, outcome=c(A, B, C))
  res <- model_parameters(anova(lm(outcome ~ group, data = dat)))
  F <- res$Mean_Square[1]/res$Mean_Square[2]
  f_unbiased[i] <- sqrt((k-1)*(F-1)/N)
  eta2_biased[i] <- res$Sum_Squares[1]/(res$Sum_Squares[1]+res$Sum_Squares[2])
  f_biased[i] <- sqrt(eta2_biased[i]/(1-eta2_biased[i]))
}

# bias
mean(f_biased)-f
mean(f_unbiased)-f

Reference Grissom, R. J., & Kim, J. J. (2012). Effect sizes for research: Univariate and multivariate applications. Routledge.

[mattansb]

Hi @Daiki-Nakamura-git - thanks for the suggestion! This is also the same as basing Cohen's f on Omega squared:

A <- aov(mpg ~ factor(cyl), data = mtcars)

library(parameters)
library(effectsize)
res <- model_parameters(A)

F <- res$Mean_Square[1]/res$Mean_Square[2]
(f_unbiased <- sqrt((k - 1) * (F - 1) / N))
#> Error in eval(expr, envir, enclos): object 'k' not found

# Same as:
omega <- omega_squared(A)[[2]]
#> For one-way between subjects designs, partial omega squared is
#>   equivalent to omega squared. Returning omega squared.
sqrt(omega / (1 - omega))
#> [1] 1.555183

Now added as a "method" argument added in https://github.com/easystats/effectsize/pull/552/

cohens_f(A)
#> For one-way between subjects designs, partial eta squared is equivalent
#>   to eta squared. Returning eta squared.
#> # Effect Size for ANOVA
#> 
#> Parameter   | Cohen's f |      95% CI
#> -------------------------------------
#> factor(cyl) |      1.65 | [1.16, Inf]
#> 
#> - One-sided CIs: upper bound fixed at [Inf].
cohens_f(A, method = "omega")
#> For one-way between subjects designs, partial omega squared is
#>   equivalent to omega squared. Returning omega squared.
#> # Effect Size for ANOVA
#> 
#> Parameter   | Cohen's f |      95% CI
#> -------------------------------------
#> factor(cyl) |      1.56 | [1.08, Inf]
#> 
#> - Based on Omega squared.
#> - One-sided CIs: upper bound fixed at [Inf].

^{Created on 2023-07-25 with reprex v2.0.2}