In strengejacke/ggeffects: Create Tidy Data Frames of Marginal Effects for 'ggplot' from Model Outputs

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  dev = "png",
  fig.width = 7,
  fig.height = 3.5,
  message = FALSE, warning = FALSE
)
options(width = 800)
arrow_color <- "#FF00cc"

pkgs <- c(
  "ggplot2",
  "marginaleffects",
  "emmeans",
  "parameters",
  "htmltools"
)

if (!all(vapply(pkgs, requireNamespace, quietly = TRUE, FUN.VALUE = logical(1L)))) {
  knitr::opts_chunk$set(eval = FALSE)
}

library(htmltools)
callout_tip <- function(header = NULL, ...) {
  div(
    class = "callout-tip",
    tags$h1(
      tags$img(src = "../man/figures/summary.png", width = "20", height = "17", style = "vertical-align:middle"), # nolint
      header
    ),
    ...
  )
}
includeCSS("../man/figures/callout.css")

This vignette is roughly a duplication of the first vignette about Contrasts and Pairwise Comparisons, but demonstrating the different backends for the calculation of pairwise comparisons. The default backend is the marginaleffects package. If desired, engine = "emmeans" can be used to switch to the emmeans package. engine = "ggeffects" is currently experimental and work-in-progress. It is less feature-rich than the other backends, but it also works for predictions of random effects.

callout_tip(
  "Summary of most important points:",
  tags$ul(
    tags$li("The ", tags$em("ggeffects-package"), " relies on the ", tags$em("marginaleffects-package"), " by default to calculate contrasts and pairwise comparisons."), # nolint
    tags$li("Although this covers many different ways to test contrasts and comparisons, sometimes it can be convenient or necessary to calculate specific contrasts, like consecutive, interaction or custom contrasts. In this case (e.g. when ", tags$code("test=\"consecutive\""), ", or a data frame for custom contrasts), ", tags$code("test_predictions()"), " automatically switches the backend to the ", tags$em("emmeans-package"), "."), # nolint
    tags$li("The backend can also be changed explicitly by using the ", tags$code("engine"), " argument. Usually, this is not necessary, unless you want to calculate contrasts of random effects levels. This is currently only possible for ", tags$code("engine=\"ggeffects\""), " (see vignette \"Case Study: Intersectionality Analysis Using The MAIHDA Framework\").") # nolint
  )
)

Within episode, do levels differ?

library(ggeffects)
library(ggplot2)

set.seed(123)
n <- 200
d <- data.frame(
  outcome = rnorm(n),
  grp = as.factor(sample(c("treatment", "control"), n, TRUE)),
  episode = as.factor(sample(1:3, n, TRUE)),
  sex = as.factor(sample(c("female", "male"), n, TRUE, prob = c(0.4, 0.6)))
)
model1 <- lm(outcome ~ grp + episode + grp, data = d)

Predictions

my_pred <- predict_response(model1, "episode", margin = "marginalmeans")
my_pred

Pairwise comparisons

p <- plot(my_pred)
line_data <- as.data.frame(my_pred, terms_to_colnames = FALSE)[1:2, ]
p + geom_segment(
  data = line_data,
  aes(
    x = as.numeric(x[1]) + 0.06, xend = as.numeric(x[2]) - 0.06,
    y = predicted[1], yend = predicted[2], group = NULL, color = NULL
  ),
  color = arrow_color,
  arrow = arrow(length = unit(0.1, "inches"), ends = "both", angle = 40)
) +
ggtitle("Within \"episode\", do levels 1 and 2 differ?")

# comparisons based on estimated marginal means, using "marginaleffects" package
test_predictions(model1, "episode", margin = "marginalmeans")

# comparisons using "emmeans" package
test_predictions(model1, "episode", engine = "emmeans")

# comparisons using "ggeffects" backend. This engine requires the
# ggeffects-object as input
test_predictions(my_pred, engine = "ggeffects")

Does same level of episode differ between groups?

model2 <- lm(outcome ~ grp * episode + grp, data = d)

Predictions

my_pred <- predict_response(model2, c("episode", "grp"), margin = "marginalmeans")
my_pred

Pairwise comparisons

p <- plot(my_pred)
line_data <- as.data.frame(my_pred, terms_to_colnames = FALSE)[3:4, 1:2]
line_data$group_col <- "control"
p + geom_segment(
  data = line_data,
  aes(
    x = as.numeric(x[1]) - 0.06, xend = as.numeric(x[2]) + 0.06,
    y = predicted[1], yend = predicted[2], group = NULL, color = NULL
  ),
  color = arrow_color,
  arrow = arrow(length = unit(0.1, "inches"), ends = "both", angle = 40)
) +
ggtitle("Within level 2 of \"episode\", do treatment and control group differ?")

# we want "episode = 2-2" and "grp = control-treatment"

# comparisons based on estimated marginal means, using "marginaleffects" package
test_predictions(model2, c("episode [2]", "grp"), margin = "marginalmeans")

# comparisons based using "emmeans" package
test_predictions(model2, c("episode [2]", "grp"), engine = "emmeans")

# comparisons using "ggeffects" backend
my_pred <- predict_response(model2, c("episode [2]", "grp"), margin = "marginalmeans")
test_predictions(my_pred, engine = "ggeffects")

Does difference between two levels of episode in the control group differ from difference of same two levels in the treatment group?

The test argument also allows us to compare difference-in-differences. When engine = "emmeans" or "ggeffects", we need to set test = "interaction" to get interaction contrasts, i.e. differences-in-differences.

my_pred <- predict_response(model2, c("grp", "episode"))
p <- plot(my_pred)
line_data <- as.data.frame(my_pred, terms_to_colnames = FALSE)[, 1:2, ]
line_data$group_col <- "1"
p + geom_segment(
  data = line_data,
  aes(
    x = as.numeric(x[1]) - 0.05, xend = as.numeric(x[1]) - 0.05,
    y = predicted[1], yend = predicted[2], group = NULL, color = NULL
  ),
  color = "orange",
  arrow = arrow(length = unit(0.1, "inches"), ends = "both", angle = 40, type = "closed")
) + geom_segment(
  data = line_data,
  aes(
    x = as.numeric(x[4]) - 0.05, xend = as.numeric(x[4]) - 0.05,
    y = predicted[4], yend = predicted[5], group = NULL, color = NULL
  ),
  color = "orange",
  arrow = arrow(length = unit(0.1, "inches"), ends = "both", angle = 40, type = "closed")
) + geom_segment(
  data = line_data,
  aes(
    x = as.numeric(x[1]) - 0.05, xend = as.numeric(x[4]) - 0.05,
    y = (predicted[1] + predicted[2]) / 2,
    yend = (predicted[4] + predicted[5]) / 2, group = NULL, color = NULL
  ),
  color = arrow_color,
  arrow = arrow(length = unit(0.1, "inches"), ends = "both", angle = 40)
) +
ggtitle("Differnce-in-differences")

# specifying the difference-in-difference when using "marginaleffects"
test_predictions(model2, c("episode", "grp"), test = "(b1 - b3) = (b2 
- b4)", margin = "marginalmeans")

# "emmeans" provides similar comparisons when we set test = "interaction".
# This displays *all* possible differences-in-differences. The first row in
# this output is identical to the above result from "marginaleffects". The
# "emmeans" package is used automatically, when test = "interaction".
test_predictions(model2, c("episode", "grp"), test = "interaction")

# using "ggeffects", we also need to set test = "interaction" to get the same
# results. However, since by default "emmeans" us used, we also need to specify
# the "engine" argument
my_pred <- predict_response(model2, c("episode", "grp"), margin = "marginalmeans")
test_predictions(my_pred, test = "interaction", engine = "ggeffects")

strengejacke/ggeffects documentation built on Dec. 24, 2024, 3:27 a.m.