Use {marginaleffects} as default backend

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: modelbased
 Title: Estimation of Model-Based Predictions, Contrasts and Means
-Version: 0.8.9.31
+Version: 0.8.9.32
 Authors@R:
     c(person(given = "Dominique",
              family = "Makowski",
@@ -89,4 +89,4 @@ Config/testthat/parallel: true
 Roxygen: list(markdown = TRUE)
 Config/Needs/check: stan-dev/cmdstanr
 Config/Needs/website: easystats/easystatstemplate
-Remotes: easystats/bayestestR
+Remotes: easystats/bayestestR, vincentarelbundock/marginaleffects

R/estimate_means.R

@@ -6,10 +6,12 @@
 #' [estimate_slopes()].
 #'
 #' @param model A statistical model.
-#' @param by The predictor variable(s) at which to evaluate the desired effect
-#' / mean / contrasts. Other predictors of the model that are not included
-#' here will be collapsed and "averaged" over (the effect will be estimated
-#' across them).
+#' @param by The (focal) predictor variable(s) at which to evaluate the desired
+#' effect / mean / contrasts. Other predictors of the model that are not
+#' included here will be collapsed and "averaged" over (the effect will be
+#' estimated across them). `by` can a character (vector) naming the focal
+#' predictors (and optionally, representative values or levels), or a list of named
+#' elements. See details in [`insight::get_datagrid()`].
 #' @param predict Is passed to the `type` argument in `emmeans::emmeans()` (when
 #' `backend = "emmeans"`) or in `marginaleffects::avg_predictions()` (when
 #' `backend = "marginaleffects"`). For emmeans, see also
@@ -77,8 +79,20 @@
 #' as default backend.
 #' @param transform Deprecated, please use `predict` instead.
 #' @param verbose Use `FALSE` to silence messages and warnings.
-#' @param ... Other arguments passed, for instance, to [insight::get_datagrid()]
-#' or functions from the **emmeans** or **marginaleffects** package.
+#' @param ... Other arguments passed, for instance, to [insight::get_datagrid()],
+#' to functions from the **emmeans** or **marginaleffects** package, or to process
+#' Bayesian models via [bayestestR::describe_posterior()]. Examples:
+#' - `insight::get_datagrid()`: Argument such as `length` or `range` can be used
+#'   to control the (number of) representative values.
+#' - **marginaleffects**: Internally used functions are `avg_predictions()` for
+#'   means and contrasts, and `avg_slope()` for slopes. Therefore, arguments
+#'   for instance like `vcov`, `transform`, `equivalence` or `slope` can be
+#'   passed to those functions.
+#' - **emmeans**: Internally used functions are `emmeans()` and `emtrends()`.
+#'   Additional arguments can be passed to these functions.
+#' - Bayesian models: For Bayesian models, parameters are cleaned using
+#'   `describe_posterior()`, thus, arguments like, for example, `centrality`,
+#'   `rope_range`, or `test` are passed to that function.
 #'
 #' @inheritParams parameters::model_parameters.default
 #' @inheritParams estimate_expectation

R/estimate_slopes.R

@@ -129,8 +129,12 @@ estimate_slopes <- function(model,
   info <- attributes(estimated)
 
   # Table formatting
+  table_footer <- paste("\nMarginal effects estimated for", info$trend)
+  if (!is.null(attributes(trends)$slope)) {
+    table_footer <- paste0(table_footer, "\nType of slope was ", attributes(trends)$slope)
+  }
   attr(trends, "table_title") <- c("Estimated Marginal Effects", "blue")
-  attr(trends, "table_footer") <- c(paste("Marginal effects estimated for", info$trend), "blue")
+  attr(trends, "table_footer") <- c(table_footer, "blue")
 
   # Add attributes
   attr(trends, "model") <- model

R/format.R

@@ -10,7 +10,7 @@
   }
 
   # arrange columns (not for contrast now)
   by <- rev(attr(x, "focal_terms", exact = TRUE))
   if (!is.null(by) && all(by %in% colnames(x))) {
     # arrange predictions
     x <- datawizard::data_arrange(x, select = by)
@@ -114,6 +114,13 @@
   if ("term" %in% colnames(x) && insight::n_unique(x$term) == 1) {
     remove_columns <- c("Parameter", remove_columns)
   }
+  # there are some exceptions for `estimate_slope()`, when the `Comparison`
+  # column contains information about the type of slope (dx/dy etc.). we want
+  # to remove this here, but add information as attribute.
+  if ("contrast" %in% colnames(x) && all(x$contrast %in% .marginaleffects_slopes())) {
+    remove_columns <- c("Comparison", "contrast", remove_columns)
+    attr(x, "slope") <- unique(x$contrast)
+  }
   # reshape and format columns
   params <- .standardize_marginaleffects_columns(
     x,
@@ -132,10 +139,11 @@
 
 #' @export
 format.marginaleffects_contrasts <- function(x, model, p_adjust, comparison, ...) {
   predict <- attributes(x)$predict
   by <- attributes(x)$by
   contrast <- attributes(x)$contrast
   focal_terms <- attributes(x)$focal_terms
+  dgrid <- attributes(x)$datagrid
 
   # clean "by" and contrast variable names, for the special cases. for example,
   # if we have `by = "name [fivenum]"`, we just want "name"
@@ -163,7 +171,7 @@
   if (!is.null(comparison) && is.character(comparison) && comparison %in% valid_options) {
     #  the goal here is to create tidy columns with the comparisons.
     # marginaleffects returns a single column that contains all levels that
-    # are contrastet. We want to have the contrasted levels per predictor in
+    # are contrasted. We want to have the contrasted levels per predictor in
     # a separate column. This is what we do here...
 
     # split parameter column into comparison groups.
@@ -172,6 +180,38 @@
       lapply(x$Parameter, .split_at_minus_outside_parentheses)
     ))
 
+    # When we filter contrasts, e.g. `contrast = c("vs", "am='1'")` or
+    # `contrast = c("vs", "am"), by = "gear='5'"`, we get no contrasts if one
+    # of the focal terms only has one unique value in the data grid. Thus,
+    # we need to exclude all those focal terms that only have one unique value
+    # in the data grid now. Fingers crossed that it works...
+    focal_terms <- focal_terms[lengths(lapply(dgrid[focal_terms], unique)) > 1]
+
+    # in the second example, `contrast = c("vs", "am"), by = "gear='5'"`, the
+    # `by` column is the one with one unique value only, we thus have to update
+    # `by` as well, and also `contrast` (the latter not(!) for numerics)...
+    by <- by[lengths(lapply(dgrid[by], unique)) > 1]
+
+    # for contrasts, we also filter variables with one unique value, but we
+    # keep numeric variables. When these are hold constant in the data grid,
+    # they are set to their mean value - meaning, they only have one unique
+    # value in the data grid, anyway. so we need to keep them
+    keep_contrasts <- lengths(lapply(dgrid[contrast], unique)) > 1 | vapply(dgrid[contrast], is.numeric, logical(1)) # nolint
+    contrast <- contrast[keep_contrasts]
+
+    # set to NULL, if all by-values have been removed here
+    if (!length(by)) by <- NULL
+
+    # if we have no contrasts left, e.g. due to `contrast = "time = factor(2)"`,
+    # we error here - we have no contrasts to show
+    if (!length(contrast)) {
+      insight::format_error("No contrasts to show. Please adjust `contrast`.")
+    }
+
+    # contrasts can't be longer than focal terms - make sure we have not
+    # removed too much (and that we now have captured all exceptions...)
+    if (length(contrast) > length(focal_terms)) focal_terms <- contrast
+
     # for more than one term, we have comma-separated levels.
     if (length(focal_terms) > 1) {
       # we now have a data frame with each comparison-pairs as single column.
@@ -201,7 +241,7 @@
       # unite back columns with focal contrasts - only needed when not slopes
       if (inherits(x, "estimate_slopes")) {
         contrast <- by
         by <- NULL
       }
 
       for (i in seq_along(contrast)) {
@@ -315,7 +355,7 @@
   params <- params[c(setdiff(colnames(params), relocate_columns), relocate_columns)]
 
   # relocate focal terms to the beginning
   by <- attr(x, "focal_terms", exact = TRUE)
   if (!is.null(by) && all(by %in% colnames(params))) {
     params <- datawizard::data_reorder(params, by, verbose = FALSE)
   }
@@ -431,7 +471,7 @@
     if (substring(input_string, match_positions[i], match_positions[i]) == "-") {
       inside_parentheses <- FALSE
       for (j in seq_along(match_positions)) {
         if (i != j && match_positions[i] > match_positions[j] && match_positions[i] < (match_positions[j] + match_lengths[j])) {
           inside_parentheses <- TRUE
           break
         }

R/get_marginalmeans.R

@@ -184,7 +184,8 @@ get_marginalmeans <- function(model,
 .info_elements <- function() {
   c(
     "at", "by", "focal_terms", "adjusted_for", "predict", "trend", "comparison",
-    "contrast", "marginalize", "p_adjust", "datagrid", "preserve_range", "coef_name"
+    "contrast", "marginalize", "p_adjust", "datagrid", "preserve_range",
+    "coef_name", "slope"
   )
 }
 

R/get_marginaltrends.R

@@ -122,3 +122,9 @@ get_marginaltrends <- function(model,
 
   trend
 }
+
+
+#' @keywords internal
+.marginaleffects_slopes <- function() {
+  c("dY/dX", "eY/eX", "eY/dX", "dY/eX")
+}

R/visualisation_recipe_internal.R

@@ -3,40 +3,49 @@
 
 #' @keywords internal
 .find_aes <- function(x) {
   data <- as.data.frame(x)
+  data$.group <- 1
+
   att <- attributes(x)
   aes <- list(
     y = "Predicted",
-    group = 1
+    group = ".group"
   )
 
   # extract information for labels
   model_data <- .safe(insight::get_data(attributes(x)$model, verbose = FALSE))
   model_response <- attributes(x)$response
 
+  # Find predictors
+  by <- att$focal_terms
+
   # Main geom
   if ("estimate_contrasts" %in% att$class) {
     insight::format_error("Automated plotting is not yet implemented for this class.")
   } else if ("estimate_means" %in% att$class) {
     aes$y <- att$coef_name
   } else if ("estimate_slopes" %in% att$class) {
     aes$y <- "Slope"
+    if ("Comparison" %in% names(data)) {
+      # Insert "Comparison" column as the 2nd by so that it gets plotted as color
+      if (length(by) > 1) by[3:(length(by) + 1)] <- by[2:length(by)]
+      by[2] <- "Comparison"
+    }
   } else if ("estimate_grouplevel" %in% att$class) {
     aes$x <- "Level"
     aes$y <- "Coefficient"
     aes$type <- "grouplevel"
     if (length(unique(data$Parameter)) > 1) {
       aes$color <- "Parameter"
-      aes$group <- "Parameter"
+      data$.group <- paste(data$.group, data$Parameter)
     }
     if (length(unique(data$Group)) > 1) aes$facet <- "Group"
     aes <- .find_aes_ci(aes, data)
     return(list(aes = aes, data = data))
   }
 
-  # Find predictors
-  by <- att$focal_terms
-  # 2nd try
+
+  # Assign predictors to aes
   if (is.null(by)) {
     by <- att$by
   }
@@ -53,16 +62,15 @@
   }
   if (length(by) > 1) {
     aes$color <- by[2]
-    aes$group <- by[2]
+    data$.group <- paste(data$.group, data[[by[2]]])
   }
   if (length(by) > 2) {
     if (is.numeric(data[[by[3]]])) {
       aes$alpha <- by[3]
     } else {
       aes$facet <- stats::as.formula(paste("~", paste(utils::tail(by, -2), collapse = " * ")))
     }
-    data$.group <- paste(data[[by[2]]], "_", data[[by[3]]])
-    aes$group <- ".group"
+    data$.group <- paste(data$.group, data[[by[3]]])
   }
   if (length(by) > 3) {
     aes$facet <- NULL
@@ -305,7 +313,17 @@
     stroke = stroke
   )
 
-  # set default alpha, it not mapped by aes
+  # check if we have matching columns in the raw data - some functions,
+  # likes slopes, have mapped these aes to other columns that are not part
+  # of the raw data - we set them to NULL
+  if (!is.null(aes$color) && !aes$color %in% colnames(rawdata)) {
+    out$aes$color <- NULL
+  }
+  if (!is.null(aes$alpha) && !aes$alpha %in% colnames(rawdata)) {
+    out$aes$alpha <- NULL
+  }
+
+  # set default alpha, if not mapped by aes
   if (is.null(aes$alpha)) {
     out$alpha <- 1 / 3
   } else {