Is there an equivalent to the del.eff method in emmeans?

[walrossker]

Hi @vincentarelbundock, for a categorical predictor with 3+ levels, I'm hoping to compare each group's mean to the mean of everyone else in the sample that's not in the group. Below is an example of how this can be done in emmeans (i.e., weigh each group by its size since the default reference group is an average of the other group averages).

Is there a way to accomplish this using the marginaleffects package? Similar question for the eff method in emmeans which compares each group to the overall sample mean.

library(tidyverse)

# Create dataset
dat <- mtcars %>% 
  mutate(across(cyl, factor))

# Run model
mod <- lm(mpg ~ cyl, data = dat)

# Run weighted contrast to compare mean of each group 
# to the mean of everyone else
emm <- emmeans::emmeans(mod, "cyl")
emmeans::contrast(emm, method = "del.eff", wts = emm@grid$.wgt.)
#>  contrast    estimate   SE df t.ratio p.value
#>  cyl4 effect   10.016 1.20 29   8.349  <.0001
#>  cyl6 effect   -0.445 1.38 29  -0.323  0.7490
#>  cyl8 effect   -8.872 1.15 29  -7.725  <.0001
#> 
#> P value adjustment: fdr method for 3 tests

# Manually replicate mean differences (but not the std errors or p-values)
dat$cyl %>% 
  unique() %>% 
  sort() %>% 
  set_names() %>% 
  map_dfr(\(.x){
    
    tmp_dat <- dat %>% mutate(group = cyl == .x)
    mod <- lm(mpg ~ group, data = tmp_dat)
    broom::tidy(mod) %>% 
      filter(term == "groupTRUE") %>% 
      select(-term)
    
  }, .id = "cyl")
#> # A tibble: 3 × 5
#>   cyl   estimate std.error statistic      p.value
#>   <chr>    <dbl>     <dbl>     <dbl>        <dbl>
#> 1 4       10.0        1.36     7.35  0.0000000344
#> 2 6       -0.445      2.62    -0.170 0.866       
#> 3 8       -8.87       1.46    -6.06  0.00000117

[vincentarelbundock]

Hi @walrossker,

Unfortunately, there currently there isn't an easy shortcut for this. I opened a issue for the feature request here: #1092

In that issue, you'll see that I included some code to replicate your emmeans results. It works, but it is not exactly user-friendly. I included it there not because I expect people to use that, but because it can give us some ideas for internal implementation of a shortcut.

There are a couple tricky challenges for implementation, but those can probably be overcome.

[vincentarelbundock]

Also note that the only reason this is anything more than trivial is that you are using the wts argument in emmeans. Without it, we fan do things like:

library(emmeans)
library(marginaleffects)
dat <- transform(mtcars, cyl = factor(cyl))
mod <- lm(mpg ~ cyl, data = dat)

emm <- emmeans::emmeans(mod, "cyl")
emmeans::contrast(emm, method = "eff")

avg_predictions(mod, by = "cyl", transform = \(z) z - mean(z))

[walrossker]

Thanks for looking into this. Using your solution in #1092, I noticed that results between packages don't match when there is more than one variable in the model. Maybe the solution wasn't designed for this scenario, or maybe this difference is expected/desirable given a underlying difference in how emmeans and marginaleffects operate (i.e., holding covariates at mean/typical values vs. evaluating all original rows then aggregating)?

# Create dataset
dat <- transform(mtcars, cyl = factor(cyl))

# Run model
mod <- lm(mpg ~ cyl + disp, data = dat)

# Run weighted contrast to compare mean of each group 
# to the mean of everyone else
emm <- emmeans::emmeans(mod, "cyl")
emmeans::contrast(emm, method = "del.eff", wts = emm@grid$.wgt.)
#>  contrast    estimate   SE df t.ratio p.value
#>  cyl4 effect     4.79 2.31 28   2.075  0.1418
#>  cyl6 effect    -2.10 1.42 28  -1.485  0.2233
#>  cyl8 effect    -2.93 2.54 28  -1.156  0.2576
#> 
#> P value adjustment: fdr method for 3 tests

# Run the `marginaleffects` solution in issue #1092 
assert_nrow_equal <- function(x, newdata) {
  if (nrow(x) != nrow(newdata)) {
    msg <- c(
      "`newdata` must have the same number of rows as the output.",
      "Do not use `avg_` or `by` with this hypothesis function."
    )
    insight::format_error(msg)
  }
  return(invisible(TRUE))
}
hyp <- function(x, newdata) {
  assert_nrow_equal(x, newdata)
  out <- cbind(x, newdata)
  s <- sum(out$estimate)
  f <- function(z) mean(z - (s - sum(z)) / (nrow(x) - length(z)))
  out <- aggregate(estimate ~ cyl, FUN = f, data = out)
  out$term <- out$cyl
  out
}
marginaleffects::predictions(mod, hypothesis = hyp)
#> 
#>  Term Estimate Std. Error      z Pr(>|z|)    S  2.5 % 97.5 %
#>     4   10.016       1.10  9.123   <0.001 63.6   7.86  12.17
#>     6   -0.445       1.26 -0.353    0.724  0.5  -2.92   2.03
#>     8   -8.872       1.05 -8.441   <0.001 54.8 -10.93  -6.81
#> 
#> Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, cyl 
#> Type:  response

[vincentarelbundock]

I don’t know what emmeans does with the wts argument in this context, but I can show you a simple base R replication of marginaleffects:

dat <- transform(mtcars, cyl = factor(cyl))
mod <- lm(mpg ~ cyl + disp, data = dat)

p <- predict(mod)
del <- mean(p[mtcars$cyl != "6"])
mean(p[mtcars$cyl == "6"] - del)

    [1] -0.4451429

[walrossker]

That makes sense. Below is a close replication of the emmeans output using your solution (only p-values differ). So it looks like the difference is mostly due to emmeans setting disp to its mean rather than evaluating all original rows as marginaleffects does:

dat <- transform(mtcars, cyl = factor(cyl))

mod <- lm(mpg ~ cyl + disp, data = dat)

assert_nrow_equal <- function(x, newdata) {
  if (nrow(x) != nrow(newdata)) {
    msg <- c(
      "`newdata` must have the same number of rows as the output.",
      "Do not use `avg_` or `by` with this hypothesis function."
    )
    insight::format_error(msg)
  }
  return(invisible(TRUE))
}
hyp <- function(x, newdata) {
  assert_nrow_equal(x, newdata)
  out <- cbind(x, newdata)
  s <- sum(out$estimate)
  f <- function(z) mean(z - (s - sum(z)) / (nrow(x) - length(z)))
  out <- aggregate(estimate ~ cyl, FUN = f, data = out)
  out$term <- out$cyl
  out
}

marginaleffects::predictions(
  mod, 
  hypothesis = hyp, 
  newdata = marginaleffects::datagrid(disp = mean,
                                      grid_type = "counterfactual"))
#> 
#>  Term cyl Estimate Std. Error     z Pr(>|z|)   S  2.5 % 97.5 %
#>     4   4     4.79       2.31  2.08    0.038 4.7  0.266  9.314
#>     6   6    -2.10       1.42 -1.48    0.138 2.9 -4.878  0.673
#>     8   8    -2.93       2.54 -1.16    0.248 2.0 -7.902  2.040
#> 
#> Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, cyl 
#> Type:  response

[vincentarelbundock]

Ah yeah, that makes sense. Thanks for the investigation!