Comparing Total MEs across models
Example 5.2.e of Mize and Han (2025), in R: mediation/attenuation
Does the total effect of a college degree on self-rated health shrink once family income is in the model? Fit the model without income and the model with it, and combine them with suest to test the difference in the Total ME. The same example in Stata is on the totalme page.
Load and prepare the data
library(haven) # Read Stata data
library(marginaleffects) # Marginal effects and hypotheses
library(nnet) # Multinomial logit
library(suest) # Combine models
source("https://raw.githubusercontent.com/tdmize/Rfunctions/main/ME_helper_functions.R")
gss <- read_dta("https://tdmize.github.io/data/data/cda_gss.dta")
gss <- gss[gss$year >= 2000 & gss$year <= 2021, ]
vars <- c("healthR", "college", "race4", "age", "woman", "parent", "married", "faminc")
gss <- gss[complete.cases(gss[vars]), vars]
fvars <- c("healthR", "college", "race4", "woman", "parent", "married")
gss[fvars] <- lapply(gss[fvars], as_factor)
nrow(gss)[1] 19292Fit and combine the models
basemod <- multinom(healthR ~ college + race4 + age + woman + parent +
married, data = gss, trace = FALSE)
medmod <- multinom(healthR ~ college + race4 + age + woman + parent +
married + faminc, data = gss, trace = FALSE)
fit <- suest(basemod, medmod, model_names = c("Base", "Income"))Total ME in each model
totalme() returns the Total ME in each model:
tme <- avg_comparisons(fit,
variables = list(college = "reference"),
newdata = gss,
hypothesis = totalme(fit))
tme Group Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
Base 0.159 0.00605 26.3 <0.001 505.4 0.147 0.171
Income 0.118 0.00673 17.5 <0.001 225.7 0.105 0.131
Term: college
Type: responseTest the difference
hypotheses(tme, hypothesis = difference ~ revpairwise) Hypothesis Estimate Std. Error z Pr(>|z|) S 2.5 % 97.5 %
(Base) - (Income) 0.0414 0.00276 15 <0.001 166.6 0.036 0.0468The difference is the part of the total college effect that family income accounts for, with its standard error.