Confidence interval for eta^2 #611

Activity

1. mattansb commented on Sep 29, 2023

   

   mattansb

   on Sep 29, 2023

   Member

   More actions

   Here is a reprex for your example, which compares effectsize with MBESS (both using non-central F distributions):

   data(obk.long, package = "afex")
   # modify the data slightly for the demonstration:
   obk.long <- obk.long[1:240 %% 3 == 0, ]
   obk.long$id <- seq_len(nrow(obk.long))
   
   model <- aov(value ~ treatment, data = obk.long)
   
   effectsize::eta_squared(model, alternative = "two.sided", ci = 0.95) |> 
     print(digits = 8)
   #> For one-way between subjects designs, partial eta squared is equivalent
   #>   to eta squared. Returning eta squared.
   #> # Effect Size for ANOVA
   #> 
   #> Parameter |       Eta2 |                   95% CI
   #> -------------------------------------------------
   #> treatment | 0.22348040 | [0.07209814, 0.36848290]
   
   A <- anova(model)
   MBESS::ci.omega2(F.value = A[1,"F value"],
                    df.1 = A[1,"Df"],
                    df.2 = A[2,"Df"],
                    N = nrow(obk.long))
   #> $lower_limit_omega2
   #> [1] 0.0695825
   #> 
   #> $upper_limit_omega2
   #> [1] 0.3596343

   ^{Created on 2023-09-29 with reprex v2.0.2}

   These values are very close, I would say the differences are negligible (2-3% difference in magnitude), and can be attributed to differences in tolerances of the function the finds the non-central F values, as well as MBESS using $\eta^2_p=\frac{\lambda}{\lambda + N}$ (as presented in Steiger, 2004) where as effectsize uses $\eta^2_p=\frac{\lambda}{\lambda + df_2}$, which we believe is more appropriate since it follows the actual conversion of the F statistic to $\eta^2_p$.

   The examples in Steiger, 2004 are hard to generalize to complex within/mixed/non-balanced designs.
2. kazz530 commented on Sep 29, 2023

   

   kazz530

   on Sep 29, 2023

   Author

   More actions

   Thank you for your quick reply.
   I have additional questions.
   How can we confirm the complete document that contains all estimation equations used in your effectsize package?
   Do you have any theoretical justification paper to use DF2 instead of N?
3. mattansb commented on Sep 29, 2023

   

   mattansb

   on Sep 29, 2023

   Member

   More actions

   The justification is based on the fact that the (exact) F to $\eta^2_p$ conversion is F *d1 / (F *d1 + df2), and so, substituting ncp for F * df1, it would yield ncp / (ncp + df2).
   Additionally, methods for obtaining CIs for partial correlations (e.g., RVAideMemoire::pcor.test()) do not use N, only df2. Since $\eta^2_p$ is a squared partial correlation, I figure it would be analogous.

   But here is a little simulation I ran that seems to suggest that even is multi-factor model in small samples the difference is negligible.

   # function to extract CIs from both packages:
   MBESS_eta2p_CI <- function(model) {
     A <- anova(model)
     ci <- MBESS::ci.omega2(F.value = A[1,"F value"],
                            df.1 = A[1,"Df"],
                            df.2 = A[nrow(A),"Df"],
                            N = insight::n_obs(model), 
                            conf.level = 0.9)
     unlist(ci)
   }
   
   effectsize_eta2p_CI <- function(model) {
     es <- effectsize::eta_squared(model, alternative = "two", ci = 0.9, verbose = FALSE)
     ci <- es[1,c("CI_low", "CI_high")]
     unlist(ci)
   }
   
   # Test if a value is within a range
   is_between <- function(x, range) {
     range[1] < x && x < range[2]
   }
   
   data(obk.long, package = "afex")
   obk.long <- obk.long[1:240 %% 3 == 0, ]
   
   model <- aov(value ~ treatment + phase + gender + hour, data = obk.long)
   
   # We will treat the values obtained from this sample as the true population 
   # so the TRUE effect size is this:
   true_es <- effectsize::eta_squared(model, ci = NULL)[1, 2]
   
   # We will simulate data from the model
   sim_model <- function(model) {
     dat <- insight::get_data(model)
     dat[insight::find_response(model)] <- simulate(model, nsim = 1)
     
     update(model, data = dat)
   }
   
   # Results:
   set.seed(1234)
   res <- replicate(1000, {
     mod <- sim_model(model)
     
     c(MBESS = is_between(true_es, MBESS_eta2p_CI(mod)),
       effectsize = is_between(true_es, effectsize_eta2p_CI(mod)))
   })
   
   mean(res["MBESS",])
   #> [1] 0.901
   mean(res["effectsize",])
   #> [1] 0.902

   Both methods obtain the same coverage rates.
4. bwiernik commented on Sep 30, 2023

   

   bwiernik

   on Sep 30, 2023

   Contributor

   More actions

   There is no theoretical justification for N. The distribution of eta^2 follows the F distribution with the appropriate df.

   Steiger used a maximum likelihood asymptotic approximation. Using the df is less biased and more accurate
5. mattansb commented on Oct 1, 2023

   

   mattansb

   on Oct 1, 2023

   Member

   More actions

   Thanks @bwiernik - this makes a lot of sense!
6. easystats

   locked and limited conversation to collaborators on Oct 1, 2023
7. mattansb

   converted this issue into a discussion #612 on Oct 1, 2023