Hedges’ g for MLMs – exploring different methods

library(lmeInfo)
library(nlme)
library(lme4)
library(eefAnalytics)
library(tictoc)
library(beepr)
library(brms)
library(tidybayes)
library(tidyverse)

all_results <- list()

# vec_res order is estimate, lower, and upper 95% interval
add_result <- function(name, vec_res) {
  stopifnot(is.vector(vec_res))
  stopifnot(length(vec_res) == 3)
  stopifnot(vec_res[1] >= vec_res[2])
  stopifnot(vec_res[1] <= vec_res[3])
  stopifnot(vec_res[2] <= vec_res[3])
  res <- tibble(
    source  = name,
    est     = vec_res[1],
    lower95 = vec_res[2],
    upper95 = vec_res[3]
  )
  
  all_results <<- append(all_results, list(res))
}

clear_results <- function() {
  all_results <<- list()
}

tibble_results <- function() {
  bind_rows(all_results)
}

eefAnalytics

output1 <- crtFREQ(
  Posttest ~ Intervention + Prettest,
  random = "School",
  intervention = "Intervention",
  data = crtData
)

summary(output1)

output1$ES$Intervention1[2,]

output1$Unconditional$ES$Intervention1[2,]

add_result("eefAnalytics total conditional",
           output1$ES$Intervention1[2, ] |> as.numeric())
add_result(
  "eefAnalytics total unconditional",
  output1$Unconditional$ES$Intervention1[2, ] |> as.numeric()
)

tibble_results()

Naive approach using lme4

test_mod <- lmer(Posttest ~ Intervention + Prettest +
                   (1|School), data = crtData)
test_mod |> summary()

test_mod_0 <- lmer(Posttest ~ Intervention + (1|School), data = crtData)
test_mod_0

vars <- VarCorr(test_mod_0) |> as.data.frame()
the_sd <- vars$vcov |> sum() |> sqrt()
the_sd

vars 

add_result("Naively divide CI by SD", 
c(fixef(test_mod)["Intervention"], confint(test_mod)["Intervention",]) |> as.numeric() / the_sd)

tibble_results()

g_mlm in {lmeInfo}

lme_num <- lme(
  fixed = Posttest ~ Intervention + Prettest,
  random = ~ 1 | School,
  data = crtData,
  method = "REML"
)

lme_den <- lme(
  fixed = Posttest ~ Intervention,
  random = ~ 1 | School,
  data = crtData,
  method = "REML"
)

summary(lme_num)

the_g <- g_mlm(
  lme_num,
  p_const = c(0,1,0),
  mod_denom = lme_den,
  r_const = c(1,1),
  infotype = "expected",
  separate_variances = FALSE
)

the_g

Without the df adjustment:

unadj_SE <- sqrt(the_g$SE_g_AB^2 / the_g$J_nu^2)
the_g$delta_AB
c(the_g$delta_AB - 1.96 * unadj_SE,
  the_g$delta_AB + 1.96 * unadj_SE)

“Manual” df adjustment:

the_g$J_nu * the_g$delta_AB

From the model:

the_g$g_AB
c(the_g$g_AB - 1.96 * the_g$SE_g_AB,
  the_g$g_AB + 1.96 * the_g$SE_g_AB)

add_result("lmeInfo Hedges df adjusted", c(the_g$g_AB,
                                           the_g$g_AB - 1.96 * the_g$SE_g_AB,
                                           the_g$g_AB + 1.96 * the_g$SE_g_AB))
add_result("lmeInfo unadjusted", c(the_g$delta_AB,
                                   the_g$delta_AB - 1.96 * unadj_SE,
                                   the_g$delta_AB + 1.96 * unadj_SE))

tibble_results()

Manually, with the help of {merDeriv}

my_g <- fixef(test_mod)["Intervention"] / the_sd
my_g

b_var <- vcov(test_mod)["Intervention", "Intervention"]
b_var

library(merDeriv)
ranef_var <- vcov(test_mod_0,
                  full = TRUE,
                  ranpar = "var",
                  information = "expected")[3, 3]

sum_vars <- vars$vcov |> sum()
the_se <- sqrt((b_var / sum_vars) +
       ((fixef(test_mod)["Intervention"]^2 * ranef_var) / (4 * (sum_vars)^3)))

the_se

add_result("Manual using merDeriv (no df adjust)", c(my_g,
                                      my_g - 1.96 * the_se,
                                      my_g + 1.96 * the_se))

tibble_results()

Try some bootstrapping with {lmersampler}

library(lmeresampler)

to_boot <- lmer(Posttest ~ Intervention + Prettest +
                  (1 | School), data = crtData)

each_boot <- function(mod) {
  denom_mod <- lmer(Posttest ~ Intervention + (1 | School),
                    data = model.frame(mod))
  
  vars   <- VarCorr(denom_mod) |> as.data.frame()
  the_sd <- vars$vcov |> sum() |> sqrt()
  
  the_b  <- fixef(mod)["Intervention"]
  the_es <- the_b / the_sd
  
  c(b = the_b, sd = the_sd, g = the_es)
}

each_boot(to_boot)

how_many_boots <- 5000

tic()
boo_para <- bootstrap(
  model = to_boot,
  .f    = each_boot,
  type  = "parametric",
  B = how_many_boots 
)
toc()
beep(3)

tic()
boo_resid <- bootstrap(
  model = to_boot,
  .f    = each_boot,
  type  = "residual",
  B = how_many_boots 
)
toc()
beep(2)

tic()
boo_case <- bootstrap(
  model = to_boot,
  .f    = each_boot,
  type  = "case",
  B = how_many_boots ,
  resample = c(TRUE, TRUE)
)
toc()
beep(1)

confint(boo_para, type = "perc")

confint(boo_resid, type = "perc")

confint(boo_case, type = "perc")

get_bootstrap_ests <- function(obj) {
  mean_est <- obj$replicates$g.Intervention  |> mean()
  intervals <- confint(obj, type = "perc") |>
  filter(term == "g.Intervention")
  
  c(mean_est, intervals$lower, intervals$upper) |> as.numeric()
}

add_result("bootstrap (parametric)", get_bootstrap_ests(boo_para))
add_result("bootstrap (resid)",      get_bootstrap_ests(boo_resid))
add_result("bootstrap (case)",       get_bootstrap_ests(boo_case))

tibble_results()

Go Bayesian

Model for the numerator:

bayes_mod_1 <- brm(Posttest ~ Intervention + Prettest +
                     (1 | School), data = crtData)

Model for the denominator:

bayes_mod_0 <- brm(Posttest ~ Intervention +
                     (1 | School), data = crtData)

bayes_mod_0 |> get_variables()

b_draws <- bayes_mod_1 |>
  spread_draws(b_Intervention)
sd_draws <- bayes_mod_0 |>
  spread_draws(sd_School__Intercept, sigma) |>
  mutate(total_sd = sqrt(sd_School__Intercept^2 + sigma^2))

head(b_draws)

head(sd_draws)

combined_draws <- bind_cols(b_draws |> select(b_Intervention),
                      sd_draws |> select(total_sd)) |>
  mutate(g = b_Intervention / total_sd)
head(combined_draws)

combined_draws_Bayes_mean <- combined_draws |>
  mean_hdci(g)
combined_draws_Bayes_median <- combined_draws |>
  median_hdci(g)
combined_draws_Bayes_mode <- combined_draws |>
  mode_hdci(g)

add_result("Bayesian HDCI (mean)",   combined_draws_Bayes_mean[,1:3] |> as.numeric())
add_result("Bayesian HDCI (median)", combined_draws_Bayes_median[,1:3] |> as.numeric())
add_result("Bayesian HDCI (mode)",   combined_draws_Bayes_mode[,1:3] |> as.numeric())

Try {eefAnalytics}’ Bayes model

eefBayes <- crtBayes(
  Posttest ~ Intervention + Prettest,
  random = "School",
  intervention = "Intervention",
  nsim = 4000,
  data = crtData
)

eefBayes

add_result("eefAnalytics Bayes total conditional",
           eefBayes$ES$Intervention1[2, ] |> as.numeric())
add_result(
  "eefAnalytics Bayes total unconditional",
  eefBayes$Unconditional$ES$Intervention1[2, ] |> as.numeric()
)

Summary of all approaches

tibble_results() |>
  mutate(CIwidth = upper95 - lower95) |>
  mutate(across(where(is.numeric), ~ round(.x, 2))) |>
  arrange(CIwidth, est)