-
-
Notifications
You must be signed in to change notification settings - Fork 25
Appropriate to use t_to_d() on estimated marginal mean contrasts from ANOVA? #617
|
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.
# 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)
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 |
All reactions
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 =
Replies: 1 comment 1 reply
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!
All reactions
-
🎉 1
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)
}