Calculate Cohen's d_r and mean differences in measurement scale

cohens_d_r.R

mean_difference_r <- function(x, object) UseMethod("mean_difference_r", object)

mean_difference_r.lm <- function(x, object) {
  x * sigma(object)
}

mean_difference_r.mlm <- function(x, object) {
  x * sqrt(mean(sigma(object)^2))
}

mean_difference_r.afex_aov <- function(x, object) {
  mean_difference_r.mlm(x, object$lm)
}

mean_difference_r.merMod <- function(x, object) {
  x * sigma(object)
}

mean_difference_r.merModLmerTest <- function(x, object) {
  mean_difference_r.merMod(x, object)
}

mean_difference_r.mixed <- function(x, object) {
  mean_difference_r.merMod(x, object$full_model)
}

mean_difference_r.gam <- function(x, object) {
  x * sqrt(object$scale)
}

cohens_d_r <- function(x, object) {
  1 / mean_difference_r(x)
}

These functions can be used to specify an equivalence threshold Δ in units of d_r (see this blog post by Jake Westfall) for two one-sided equivalence tests (Lakens, 2017) in the R-package lsmeans. The following example is adapted from the afex vignette ANOVA and Post-Hoc Contrasts: Reanalysis of Singmann and Klauer (2011).

library("afex")
data("sk2011.1")

source("https://gist.githubusercontent.com/crsh/44a44dee90cba1198288b8c766835809/raw/dc643cc256d58449f403eaaef0b8bd8c91c5867c/cohens_d_r.R")

# ANOVA
anova_model <- aov_ez(
  id = "id"
  , dv = "response"
  , data = sk2011.1
  , between = "instruction"
  , within = c("inference", "plausibility")
)

# Least squares means grid
anova_lsmeans <- lsmeans(anova_model, ~ inference:instruction)

# Define custom contrasts
contrast_matrix <- list(
  v_i.ded = c(0.5, 0.5, -0.5, -0.5, 0, 0, 0, 0),
  v_i.prob = c(0, 0, 0, 0, 0.5, 0.5, -0.5, -0.5)
)

# Calculate contrasts adjusted for multiple comparisons
anova_contrasts <- contrast(anova_lsmeans, contrast_matrix, adjust = "holm")

# Null hypothesis significance tests
summary(
  anova_contrasts
  , infer = TRUE
  , level = 0.95
)

# Two one-sided equivalence tests with delta = 0.5 d_r
summary(
  anova_contrasts
  , infer = TRUE
  , side = "equivalence"
  , delta = mean_difference_r(0.5, anova_model)
  , level = 0.9
)