Displaying margin effects, average margin effects or conditional effects

DESCRIPTION

@@ -36,13 +36,15 @@ Suggests:
     cmprsk,
     covr,
     datasets,
+    effects,
     emmeans,
     fixest (>= 0.10.0),
     forcats,
     gam,
     gee,
     geepack,
     ggplot2,
+    ggeffects,
     glmmTMB,
     glue,
     gt,
@@ -52,6 +54,7 @@ Suggests:
     lfe,
     lme4 (>= 1.1.28),
     logitr (>= 0.8.0),
+    margins,
     MASS,
     mgcv,
     mice,

NAMESPACE

@@ -143,12 +143,15 @@ export(tidy_add_pairwise_contrasts)
 export(tidy_add_reference_rows)
 export(tidy_add_term_labels)
 export(tidy_add_variable_labels)
+export(tidy_all_effects)
 export(tidy_and_attach)
 export(tidy_attach_model)
 export(tidy_detach_model)
 export(tidy_disambiguate_terms)
 export(tidy_get_model)
+export(tidy_ggpredict)
 export(tidy_identify_variables)
+export(tidy_margins)
 export(tidy_multgee)
 export(tidy_parameters)
 export(tidy_plus_plus)

NEWS.md

@@ -8,6 +8,9 @@
   to compute pairwise contrasts of categorical variables with `emmeans`,
   and corresponding new arguments in `tidy_plus_plus()` (#192)
 - `tidy_select_variables()` now sorts the variables according to `include` (#183)
+- New tidier `tidy_margins()` to display Average Marginal Effects (#195)
+- New tidier `tidy_all_effects()` to display Marginal Effects (#195)
+- New tidier `tidy_ggpredict()` to display Conditional Effects (#195)
 
 **New supported models**
 

R/tidiers.R

@@ -170,3 +170,142 @@ tidy_multgee <- function(x, conf.int = TRUE, conf.level = .95, ...) {
     return(res)
   }
 }
+
+#' Average Marginal Effects Estimation
+#'
+#' Use `margins::margins()` to estimate "average marginal effects" and return a
+#' tibble tidied in a way that could be used by `broom.helpers`functions. See
+#' `margins::margins()` for a list of supported models.
+#' @note When applying `margins::margins()`, custom contrasts are ignored.
+#' Treatment contrasts (`stats::contr.treatment()`) are applied to all
+#' categorical variables. Interactions are also ignored.
+#' @param x a model
+#' @param conf.int logical indicating whether or not to include a confidence
+#' interval in the tidied output
+#' @param conf.level the confidence level to use for the confidence interval
+#' @param ... additional parameters passed to `margins::margins()`
+#' @family custom_tieders
+#' @seealso `margins::margins()`
+#' @export
+#' @examples
+#' df <- Titanic %>%
+#'   dplyr::as_tibble() %>%
+#'   dplyr::mutate(Survived = factor(Survived, c("No", "Yes")))
+#' mod <- glm(
+#'   Survived ~ Class + Age + Sex,
+#'   data = df, weights = df$n, family = binomial
+#' )
+#' tidy_margins(mod)
+#' tidy_plus_plus(mod, tidy_fun = tidy_margins)
+tidy_margins <- function(x, conf.int = TRUE, conf.level = 0.95, ...) {
+  .assert_package("margins")
+
+  dots <- rlang::dots_list(...)
+  if (isTRUE(dots$exponentiate))
+    cli::cli_abort("{.arg exponentiate = TRUE} is not relevant for {.fun broom.helpers::tidy_margins}.") # nolint
+
+  res <- broom::tidy(
+    margins::margins(x, ...),
+    conf.int = conf.int,
+    conf.level = conf.level
+  )
+  attr(res, "coefficients_type") <- "average_marginal_effects"
+  res
+}
+
+#' Marginal Effects Estimation
+#'
+#' Use `effects::allEffects()` to estimate "marginal effects" and return a
+#' tibble tidied in a way that could be used by `broom.helpers`functions.
+#' See `vignette("functions-supported-by-effects", package = "effects")` for
+#' a list of supported models.
+#' @param x a model
+#' @param conf.int logical indicating whether or not to include a confidence
+#' interval in the tidied output
+#' @param conf.level the confidence level to use for the confidence interval
+#' @param ... additional parameters passed to `effects::allEffects()`
+#' @family custom_tieders
+#' @seealso `effects::allEffects()`
+#' @export
+#' @examples
+#' df <- Titanic %>%
+#'   dplyr::as_tibble() %>%
+#'   dplyr::mutate(Survived = factor(Survived, c("No", "Yes")))
+#' mod <- glm(
+#'   Survived ~ Class + Age + Sex,
+#'   data = df, weights = df$n, family = binomial
+#' )
+#' tidy_all_effects(mod)
+#' tidy_plus_plus(mod, tidy_fun = tidy_all_effects)
+tidy_all_effects <- function(x, conf.int = TRUE, conf.level = .95, ...) {
+  .assert_package("effects")
+
+  dots <- rlang::dots_list(...)
+  if (isTRUE(dots$exponentiate))
+    cli::cli_abort("{.arg exponentiate = TRUE} is not relevant for {.fun broom.helpers::tidy_all_effects}.") # nolint
+
+  .clean <- function(x) {
+    # merge first columns if interaction
+    x <- tidyr::unite(x, "term", 1:(ncol(x) - 4), sep = ":")
+    names(x) <- c("term", "estimate", "std.error", "conf.low", "conf.high")
+    x$term <- as.character(x$term)
+    rownames(x) <- NULL
+    x
+  }
+  res <- x %>%
+    effects::allEffects(se = conf.int, level = conf.level, ...) %>%
+    as.data.frame() %>%
+    purrr::map(.clean) %>%
+    dplyr::bind_rows(.id = "variable")
+  attr(res, "coefficients_type") <- "marginal_effects"
+  res
+}
+
+#' Conditional Effects Estimation
+#'
+#' Use `ggeffects::ggpredict()` to estimate "conditional effects" and return a
+#' tibble tidied in a way that could be used by `broom.helpers`functions.
+#' See <https://strengejacke.github.io/ggeffects/> for a list of supported
+#' models.
+#' @param x a model
+#' @param conf.int logical indicating whether or not to include a confidence
+#' interval in the tidied output
+#' @param conf.level the confidence level to use for the confidence interval
+#' @param ... additional parameters passed to `ggeffects::ggpredict()`
+#' @family custom_tieders
+#' @seealso `ggeffects::ggpredict()`
+#' @export
+#' @examples
+#' df <- Titanic %>%
+#'   dplyr::as_tibble() %>%
+#'   dplyr::mutate(Survived = factor(Survived, c("No", "Yes")))
+#' mod <- glm(
+#'   Survived ~ Class + Age + Sex,
+#'   data = df, weights = df$n, family = binomial
+#' )
+#' tidy_ggpredict(mod)
+#' tidy_plus_plus(mod, tidy_fun = tidy_ggpredict)
+tidy_ggpredict <- function(x, conf.int = TRUE, conf.level = .95, ...) {
+  .assert_package("ggeffects")
+
+  dots <- rlang::dots_list(...)
+  if (isTRUE(dots$exponentiate))
+    cli::cli_abort("{.arg exponentiate = TRUE} is not relevant for {.fun broom.helpers::tidy_ggpredict}.") # nolint
+
+  if (isFALSE(conf.int)) conf.level <- NA
+  res <- x %>%
+    ggeffects::ggpredict(ci.lvl = conf.level, ...) %>%
+    purrr::map(
+      ~ .x %>%
+        dplyr::as_tibble() %>%
+        dplyr::mutate(x = as.character(.data$x))
+    ) %>%
+    dplyr::bind_rows() %>%
+    dplyr::rename(
+      variable = "group",
+      term = "x",
+      estimate = "predicted"
+    )
+  attr(res, "coefficients_type") <- "conditional_effects"
+  res
+}

R/tidy_add_coefficients_type.R

@@ -41,23 +41,36 @@ tidy_add_coefficients_type <- function(
   .attributes <- .save_attributes(x)
   .attributes$exponentiate <- exponentiate
 
-  coefficients_type <- model_get_coefficients_type(model)
-  if (exponentiate)
-    coefficients_label <- dplyr::case_when(
-      coefficients_type == "logistic" ~ "OR",
-      coefficients_type == "poisson" ~ "IRR",
-      coefficients_type == "relative_risk" ~ "RR",
-      coefficients_type == "prop_hazard" ~ "HR",
-      TRUE ~ "exp(Beta)"
-    )
-  else
-    coefficients_label <- dplyr::case_when(
-      coefficients_type == "logistic" ~ "log(OR)",
-      coefficients_type == "poisson" ~ "log(IRR)",
-      coefficients_type == "relative_risk" ~ "log(RR)",
-      coefficients_type == "prop_hazard" ~ "log(HR)",
-      TRUE ~ "Beta"
-    )
+  # specific case for marginal models where coefficients_type
+  # is already define by the tidier
+  if (isTRUE(.attributes$coefficients_type == "average_marginal_effects")) {
+    coefficients_type <- "average_marginal_effects"
+    coefficients_label <- "Average Marginal Effects"
+  } else if (isTRUE(.attributes$coefficients_type == "marginal_effects")) {
+    coefficients_type <- "marginal_effects"
+    coefficients_label <- "Marginal Effects"
+  }  else if (isTRUE(.attributes$coefficients_type == "conditional_effects")) {
+    coefficients_type <- "conditional_effects"
+    coefficients_label <- "Conditional Effects"
+  } else {
+    coefficients_type <- model_get_coefficients_type(model)
+    if (exponentiate)
+      coefficients_label <- dplyr::case_when(
+        coefficients_type == "logistic" ~ "OR",
+        coefficients_type == "poisson" ~ "IRR",
+        coefficients_type == "relative_risk" ~ "RR",
+        coefficients_type == "prop_hazard" ~ "HR",
+        TRUE ~ "exp(Beta)"
+      )
+    else
+      coefficients_label <- dplyr::case_when(
+        coefficients_type == "logistic" ~ "log(OR)",
+        coefficients_type == "poisson" ~ "log(IRR)",
+        coefficients_type == "relative_risk" ~ "log(RR)",
+        coefficients_type == "prop_hazard" ~ "log(HR)",
+        TRUE ~ "Beta"
+      )
+  }
 
   attr(x, "coefficients_type") <- coefficients_type
   attr(x, "coefficients_label") <- coefficients_label

R/tidy_add_pairwise_contrasts.R

@@ -68,6 +68,13 @@ tidy_add_pairwise_contrasts <- function(
 
   .attributes <- .save_attributes(x)
 
+  if (isTRUE(.attributes$coefficients_type == "average_marginal_effects"))
+    cli::cli_abort("Pairwise contrasts are not compatible with Average Marginal Effects.") # nolint
+  if (isTRUE(.attributes$coefficients_type == "marginal_effects"))
+    cli::cli_abort("Pairwise contrasts are not compatible with Marginal Effects.") # nolint
+  if (isTRUE(.attributes$coefficients_type == "conditional_effects"))
+    cli::cli_abort("Pairwise contrasts are not compatible with Conditional Effects.") # nolint
+
   if (is.null(conf.level))
     conf.level <- .attributes$conf.level
 

R/tidy_add_reference_rows.R

@@ -66,10 +66,20 @@ tidy_add_reference_rows <- function(
     stop("'model' is not provided. You need to pass it or to use 'tidy_and_attach()'.")
   }
 
+  .attributes <- .save_attributes(x)
+
   # adding reference rows is not meaningful for stats::aov
   if (inherits(model, "aov")) {
     return(x %>% dplyr::mutate(reference_row = NA))
   }
+  # adding reference rows is not meaningful for marginal effects
+  if (isTRUE(.attributes$coefficients_type == "marginal_effects")) {
+    return(x %>% dplyr::mutate(reference_row = NA))
+  }
+  # adding reference rows is not meaningful for conditional effects
+  if (isTRUE(.attributes$coefficients_type == "conditional_effects")) {
+    return(x %>% dplyr::mutate(reference_row = NA))
+  }
 
   if ("header_row" %in% names(x)) {
     stop("`tidy_add_reference_rows()` cannot be applied after `tidy_add_header_rows().`")
@@ -84,8 +94,6 @@ tidy_add_reference_rows <- function(
     return(x)
   }
 
-  .attributes <- .save_attributes(x)
-
   if ("label" %in% names(x)) {
     if (!quiet)
       cli_alert_info(paste0(

R/tidy_add_term_labels.R

@@ -215,7 +215,27 @@ tidy_add_term_labels <- function(x,
   names(interaction_terms) <- interaction_terms
   interaction_terms <-
     interaction_terms %>%
-    strsplit(":") %>%
+    strsplit(":")
+
+  # in the case of marginal/conditional effects
+  # interaction terms are not prefixed by variable names
+  # => need to identify them from interaction_terms directly
+  if (
+    isTRUE(.attributes$coefficients_type == "marginal_effects") ||
+    isTRUE(.attributes$coefficients_type == "conditional_effects")
+  ) {
+    missing_terms <- setdiff(
+      unique(unname(interaction_terms[interaction_terms != ""])),
+      names(term_labels)
+    )
+    if (length(missing_terms) > 0) {
+      names(missing_terms) <- missing_terms
+      term_labels <- term_labels %>%
+        .update_vector(missing_terms)
+    }
+  }
+
+  interaction_terms <- interaction_terms %>%
     lapply(function(x) {
       paste(term_labels[x], collapse = interaction_sep)
     }) %>%

R/tidy_and_attach.R

@@ -31,7 +31,15 @@ tidy_attach_model <- function(x, model, .attributes = NULL) {
     dplyr::as_tibble() %>%
     .order_tidy_columns()
   class(x) <- c("broom.helpers", class(x))
-  attr(x, "model") <- model_get_model(model)
+  model <- model_get_model(model)
+
+  # if average_marginal_effects, force contr.treatment contrasts
+  if (isTRUE(attr(x, "coefficients_type") == "average_marginal_effects")) {
+    for (v in names(model$contrasts))
+      model$contrasts[[v]] <- "contr.treatment"
+  }
+
+  attr(x, "model") <- model
   for (a in names(.attributes)) {
     if (!is.null(.attributes[[a]]))
       attr(x, a) <- .attributes[[a]]

R/tidy_identify_variables.R

@@ -52,6 +52,21 @@ tidy_identify_variables <- function(x, model = tidy_get_model(x),
 
   .attributes <- .save_attributes(x)
 
+  if (
+    isTRUE(.attributes$coefficients_type == "marginal_effects") ||
+    isTRUE(.attributes$coefficients_type == "conditional_effects")
+  ) {
+    variables_list <- model_identify_variables(model) %>%
+      dplyr::select(-dplyr::all_of("term")) %>%
+      dplyr::distinct(.data$variable, .keep_all = TRUE) %>%
+      dplyr::filter(!is.na(.data$variable))
+
+    x <- x %>%
+      dplyr::left_join(variables_list, by = "variable") %>%
+      tidy_attach_model(model = model, .attributes = .attributes)
+    return(x)
+  }
+
   if ("variable" %in% names(x)) {
     x <- x %>% dplyr::select(
       -any_of(c("variable", "var_class", "var_type", "var_nlevels"))

_pkgdown.yml

@@ -50,6 +50,9 @@ reference:
       - tidy_parameters
       - tidy_with_broom_or_parameters
       - tidy_multgee
+      - tidy_margins
+      - tidy_all_effects
+      - tidy_ggpredict
 
   - title: Get information from model objects
   - contents: