bootstap_me_inequality.R

# Example of how to estimate ME inequality measure analytically and via 
# bootstrap based on:
# https://sociologicalscience.com/articles-v12-7-115/

# load libraries ----------------------------------------------------------

library(palmerpenguins)
library(marginaleffects)
library(tidyverse)

# prep the model ----------------------------------------------------------

# get rid of missing values
penguins <- penguins |> drop_na()

# make a model
model <- lm(body_mass_g ~ species + island + sex, data = penguins)

# actual comparison -------------------------------------------------------

p <- prop.table(table(penguins$species))
w_ac <- (p[1] + p[2]) / 2
w_ag <- (p[1] + p[3]) / 2
w_cg <- (p[2] + p[3]) / 2

# set vcov = FALSE if you don't want SEs and CIs here
estimate_analytical <- avg_comparisons(
  model, 
  variables = list(species = "pairwise"),
  vcov = TRUE, 
  hypothesis = "`w_ac`*abs(b1) + `w_ag`*abs(b2) + `w_cg`*abs(b3) = 0"
)

# bootstrap comparison -----------------------------------------------------

B <- 1000
me_ineq <- rep(NA, B)
for(i in 1:B) {
  # create the replicate sample
  penguins_rep <- slice_sample(penguins, n = nrow(penguins), replace = TRUE)
  
  # calculate model
  model <- lm(body_mass_g ~ species + island + sex, data = penguins_rep)
  
  # get estimate
  p <- prop.table(table(penguins_rep$species))
  w_ac <- (p[1] + p[2]) / 2
  w_ag <- (p[1] + p[3]) / 2
  w_cg <- (p[2] + p[3]) / 2
  estimate <- avg_comparisons(
    model, 
    variables = list(species = "pairwise"),
    vcov = FALSE, 
    hypothesis = "`w_ac`*abs(b1) + `w_ag`*abs(b2) + `w_cg`*abs(b3) = 0"
  )
  
  # add it to the vector
  me_ineq[i] <- estimate$estimate
}

# get SE
sd(me_ineq)
# get 95% CI
quantile(me_ineq, probs = c(0.025, 0.975))