Appropriate to use t_to_d() on estimated marginal mean contrasts from ANOVA?

[DallasNovakowski]

Question and context

I'm running some factorial anova analyses, and would like to use t_to_d() with emmeans to retrieve cohen's d effect sizes for my pairwise contrasts. I'm just wondering whether this is statistically sound? I see the degrees of freedom are the same all across the contrasts, and am not sure whether these large df values might lead to some inflated and inaccurate estimates of cohen's d.

Reprex:

library(palmerpenguins) # dataset
#> Warning: package 'palmerpenguins' was built under R version 4.1.3
library(car) # anova functions
#> Warning: package 'car' was built under R version 4.1.3
#> Loading required package: carData
#> Warning: package 'carData' was built under R version 4.1.3
library(effectsize) # effect sizes
library(dplyr)
#> Warning: package 'dplyr' was built under R version 4.1.3
#> 
#> Attaching package: 'dplyr'
#> The following object is masked from 'package:car':
#> 
#>     recode
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

# assigns data to a dataframe we call "df"
df <- palmerpenguins::penguins

# drop rows with missing values
df <- df[complete.cases(df)==TRUE, ]


df %>%
    group_by(species, sex) %>% # Group by the specified variables
    dplyr::summarise(n()) %>%
  knitr::kable()
#> `summarise()` has grouped output by 'species'. You can override using the
#> `.groups` argument.

species	sex	n()
Adelie	female	73
Adelie	male	73
Chinstrap	female	34
Chinstrap	male	34
Gentoo	female	58
Gentoo	male	61

# Fit data
flipper_fit <- stats::aov(flipper_length_mm ~ species, data = df)
# Run anova
flipper_anova <- car::Anova(flipper_fit)


# Extract estimated marginal means
flipper_emmeans <- emmeans::emmeans(flipper_fit, specs = pairwise ~ species)

# convert estimated marginal mean contrasts to dataframe
flipper_emmeans_contrasts <- data.frame(flipper_emmeans$contrasts)

knitr::kable(flipper_emmeans_contrasts)

contrast	estimate	SE	df	t.ratio	p.value
Adelie - Chinstrap	-5.72079	0.9796493	330	-5.83963	0
Adelie - Gentoo	-27.13255	0.8240767	330	-32.92479	0
Chinstrap - Gentoo	-21.41176	1.0143492	330	-21.10887	0

effectsize::t_to_d(t = flipper_emmeans_contrasts[flipper_emmeans_contrasts$contrast == "Adelie - Chinstrap", "t.ratio"],
                   df_error = flipper_emmeans_contrasts[flipper_emmeans_contrasts$contrast == "Adelie - Chinstrap", "df"]
                   )
#> d     |         95% CI
#> ----------------------
#> -0.64 | [-0.86, -0.42]

^{Created on 2023-11-05 with reprex v2.0.2}

[mattansb]

Hey, this will give approximate d effect sizes (see discussion here #212). But if you are using {emmeans} you can use emmeans::eff_size() to get d effect sizes directly:

library(palmerpenguins) # dataset
library(car) # anova functions
library(effectsize) # effect sizes
library(dplyr)

# assigns data to a dataframe we call "df"
df <- palmerpenguins::penguins

# drop rows with missing values
df <- df[complete.cases(df)==TRUE, ]

flipper_fit <- stats::aov(flipper_length_mm ~ species, data = df)

# Extract estimated marginal means
flipper_emmeans <- emmeans::emmeans(flipper_fit, specs = ~ species)

emmeans::eff_size(
  flipper_emmeans,
  method = "pairwise",
  sigma = sigma(flipper_fit),
  edf = df.residual(flipper_fit)
)
#>  contrast             effect.size    SE  df lower.CL upper.CL
#>  (Adelie - Chinstrap)      -0.857 0.151 330    -1.15   -0.561
#>  (Adelie - Gentoo)         -4.066 0.201 330    -4.46   -3.671
#>  (Chinstrap - Gentoo)      -3.209 0.197 330    -3.60   -2.822
#> 
#> sigma used for effect sizes: 6.673 
#> Confidence level used: 0.95

(This suggestion also appears in t_to_d() documentation:

The resulting d effect size is an approximation to Cohen's d, and assumes two equal group sizes. When possible, it is advised to directly estimate Cohen's d, with cohens_d(), emmeans::eff_size(), or similar functions.

I suggest reading emmeans::eff_size() docs about how to select sigma in more complex models.

Hope this helps!

[DallasNovakowski]

This is terrific. Thank you!

If it is of interest to anyone who stumbles here, I wrote a small function for calculating emmeans::eff_size() and merging it into the same dataframe as the raw contrasts:

calculate_and_merge_effect_sizes <- function(emmeans, model) {
    
  contrasts <- data.frame(emmeans$contrasts)

  emmean_d <- data.frame(emmeans::eff_size(
    emmeans,
    method = "pairwise",
    sigma = sigma(model),
    edf = df.residual(model)))
  
 combined_dataframe <- data.frame(contrasts, emmean_d)
  
    # Rename some columns for clarity
    combined_dataframe <- combined_dataframe %>%
  select(-contrast.1, -df.1)%>%
      rename(d = effect.size,
            d_ci_low = lower.CL,
            d_ci_high = upper.CL,
            d_se = SE.1,
            df_error = df,
            p = p.value)

    # Return the combined dataframe with effect size results
    return(combined_dataframe)
}