estimate_contrasts() gets a post_process argument

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: modelbased
 Title: Estimation of Model-Based Predictions, Contrasts and Means
-Version: 0.15.0.2
+Version: 0.15.0.3
 Authors@R:
     c(person(given = "Dominique",
              family = "Makowski",
@@ -114,10 +114,10 @@ VignetteBuilder:
     knitr
 Encoding: UTF-8
 Language: en-US
-RoxygenNote: 7.3.3
 Config/testthat/edition: 3
 Config/testthat/parallel: true
 Roxygen: list(markdown = TRUE)
 Config/Needs/check: stan-dev/cmdstanr
 Config/Needs/website: easystats/easystatstemplate
 LazyData: true
+Config/roxygen2/version: 8.0.0

NEWS.md

@@ -5,6 +5,9 @@
 * Informative error message when the `trend` variable in `estimate_slopes()` is
   not numeric and `backend = "emmeans"`.
 
+* `estimate_contrasts()` gets a `post_process` argument, to process subsequent
+  comparisons. It allows for complex, multi-step comparisons.
+
 # modelbased 0.15.0
 
 ## Breaking Changes

R/estimate_contrasts.R

@@ -77,6 +77,11 @@
 #'     `difference` or `ratio`, or skipped).
 #'   * If contrasts should be calculated (or grouped by) factors, `comparison`
 #'     can also be a matrix that specifies factor contrasts (see 'Examples').
+#' @param post_process Optional formula, character string or function (see
+#' `comparison`), or a list of formulas, string or functions, to process
+#' subsequent, multi-step comparisons. After the initial comparison in
+#' `comparison` is completed, the results are then post-processed using the
+#' specified post-process tests. See 'Exmaples'.
 #' @param effectsize Desired measure of standardized effect size, one of
 #' `"emmeans"`, `"marginal"`, or `"boot"`. Default is `NULL`, i.e. no effect
 #' size will be computed.
@@ -368,6 +373,34 @@
 #' # are differences in time trends of context effects statistically significant
 #' # between education levels?
 #' estimate_contrasts(model, c("phq4_within", "phq4_between", "education"))
+#'
+#' # Post-processing of multiple comparisons ---------------------
+#' # Caution! Don't expect this example to be meaningful! # It is
+#' # just to demonstrate the usage of the `post_process` argument.
+#' # -------------------------------------------------------------
+#' data("qol_cancer", package = "parameters")
+#' model <- lme4::lmer(QoL ~ time * education + (1 + time | ID), data = qol_cancer)
+#'
+#' # contrasts (pairwise comparisons) by timepoints - the default
+#' estimate_contrasts(model, "education=c('low', 'mid')", by = "time")
+#'
+#' # contrasts (pairwise comparisons) by timepoints, the default for `comparison`
+#' # additionally, we compare the differences of these contrasts across timepoints
+#' # against the reference time point (contrasts at times 2 and 3 against
+#' # contrasts at time 1)
+#' estimate_contrasts(model,
+#'   contrasts = "education=c('low', 'mid')",
+#'   by = "time",
+#'   post_process = ~reference
+#' )
+#'
+#' # multiple post-processing steps - same as before, but calculates
+#' # poly-contrasts in addition to reference contrasts
+#' estimate_contrasts(model,
+#'   contrasts = "education=c('low', 'mid')",
+#'   by = "time",
+#'   post_process = list(~reference, ~poly)
+#' )
 #' }
 #'
 #' @return A data frame of estimated contrasts.
@@ -389,6 +422,7 @@ estimate_contrasts.default <- function(
   estimate = NULL,
   p_adjust = "none",
   transform = NULL,
+  post_process = NULL,
   keep_iterations = FALSE,
   effectsize = NULL,
   iterations = 200,
@@ -441,6 +475,7 @@ estimate_contrasts.default <- function(
       ci = ci,
       estimate = estimate,
       transform = transform,
+      post_process = post_process,
       keep_iterations = keep_iterations,
       verbose = verbose,
       ...

R/get_marginalcontrasts.R

@@ -9,6 +9,7 @@ get_marginalcontrasts <- function(
   comparison = "pairwise",
   estimate = NULL,
   transform = NULL,
+  post_process = NULL,
   p_adjust = "none",
   keep_iterations = FALSE,
   verbose = TRUE,
@@ -123,6 +124,7 @@ get_marginalcontrasts <- function(
       hypothesis = my_args$comparison_slopes,
       backend = "marginaleffects",
       transform = transform,
+      post_process = post_process,
       keep_iterations = keep_iterations,
       verbose = verbose,
       ...
@@ -139,6 +141,7 @@ get_marginalcontrasts <- function(
       backend = "marginaleffects",
       estimate = estimate,
       transform = transform,
+      post_process = post_process,
       keep_iterations = keep_iterations,
       verbose = verbose,
       .joint_test = my_args$joint_test,

R/get_marginalmeans.R

@@ -40,6 +40,7 @@ get_marginalmeans <- function(
 
   dots <- list(...)
   comparison <- dots$hypothesis
+  post_process <- dots$post_process
   joint_test <- isTRUE(dots$.joint_test)
 
   # validate input
@@ -86,7 +87,7 @@ get_marginalmeans <- function(
   # fmt: skip
   dots[c(
     "by", "conf_level", "type", "digits", "bias_correction", "sigma",
-    "offset", ".joint_test"
+    "offset", ".joint_test", "post_process"
   )] <- NULL
 
   # model df - can be passed via `...`
@@ -184,7 +185,7 @@ get_marginalmeans <- function(
 
   # we can use this function for contrasts as well,
   # just need to add "hypothesis" argument
-  means <- .call_marginaleffects(fun_args)
+  means <- .call_marginaleffects(fun_args, post_process, verbose)
   vcov_means <- .safe(stats::vcov(means))
 
   # intermediate step: joint tests --------------------------------------------
@@ -253,9 +254,14 @@ get_marginalmeans <- function(
 
 # call marginaleffects and process potential errors ---------------------------
 
-.call_marginaleffects <- function(fun_args, type = "means") {
+.call_marginaleffects <- function(fun_args, post_process = NULL, verbose = TRUE) {
   out <- tryCatch(
-    suppressWarnings(do.call(marginaleffects::avg_predictions, fun_args)),
+    {
+      # call marginaleffects
+      result <- suppressWarnings(do.call(marginaleffects::avg_predictions, fun_args))
+      # process subsequential comparisons, if any
+      .post_process_comparisons(result, post_process, verbose)
+    },
     error = function(e) e
   )
 
@@ -272,6 +278,75 @@ get_marginalmeans <- function(
   out
 }
 
+# process subsequential comparisons ---------------------------
+
+.post_process_comparisons <- function(out, post_process = NULL, verbose = TRUE) {
+  if (!is.null(post_process)) {
+    # coerce to list
+    if (!is.list(post_process)) {
+      post_process <- list(post_process)
+    }
+    # save temporary parameter names
+    temp_params <- .safe(out$hypothesis)
+    # parameter names from groups
+    group_params <- .safe(out[seq_len(which(colnames(out) == "estimate") - 1)[-1]])
+    # check
+    for (i in seq_along(post_process)) {
+      if (verbose) {
+        insight::format_alert(paste0(
+          "Post-processing ",
+          deparse(post_process[[i]]),
+          "..."
+        ))
+      }
+      out <- suppressWarnings(marginaleffects::hypotheses(
+        out,
+        hypothesis = post_process[[i]]
+      ))
+      temp_params <- .update_post_process_params(out, temp_params, group_params)
+    }
+    # check if extractinb parameters worked
+    if (length(temp_params) == nrow(out)) {
+      out$hypothesis <- temp_params
+    }
+  }
+  out
+}
+
+.update_post_process_params <- function(out, temp_params, group_params) {
+  if (is.null(temp_params)) {
+    return(NULL)
+  }
+  .safe(
+    {
+      # clean current parameter names
+      temp_params <- gsub("(", "", gsub(")", "", temp_params, fixed = TRUE), fixed = TRUE)
+      # extract rownumbers (from b-cofficients) from current output
+      lhs <- as.numeric(gsub("\\(b(\\d+)\\) - \\(b(\\d+)\\)", "\\1", out$hypothesis))
+      rhs <- as.numeric(gsub("\\(b(\\d+)\\) - \\(b(\\d+)\\)", "\\2", out$hypothesis))
+      # update current parameter names
+      if (all(temp_params[lhs] == temp_params[rhs])) {
+        temp_params <- temp_params[lhs]
+      } else {
+        temp_params <- paste(temp_params[lhs], "-", temp_params[rhs])
+      }
+      if (!is.null(group_params)) {
+        for (i in colnames(group_params)) {
+          temp_params <- paste0(
+            temp_params,
+            ", ",
+            group_params[lhs, i],
+            " - ",
+            group_params[rhs, i]
+          )
+        }
+      }
+      temp_params
+    },
+    quietly = TRUE
+  )
+}
+
 
 # create data grid for marginal means ----------------------------------------
 #

R/get_marginaltrends.R

@@ -84,7 +84,7 @@ get_marginaltrends <- function(
   }
 
   # remove user-arguments from "..." that will be used when calling marginaleffects
-  dots[c("by", "conf_level", "digits")] <- NULL
+  dots[c("by", "conf_level", "digits", "post_process")] <- NULL
 
   # handle weights - argument is named "wts" in marginal effects
   if (!is.null(dots$weights)) {
@@ -136,6 +136,10 @@ get_marginaltrends <- function(
   estimated <- suppressWarnings(do.call(marginaleffects::avg_slopes, fun_args))
   vcov_slopes <- .safe(stats::vcov(estimated))
 
+  # any subsequent comparisons? -----------------------------------------------
+  # ---------------------------------------------------------------------------
+  estimated <- .post_process_comparisons(estimated, myargs$post_process)
+
   # Fourth step: back-transform response --------------------------------------
   # ---------------------------------------------------------------------------
 
@@ -199,12 +203,15 @@ get_marginaltrends <- function(
 # =========================================================================
 
 #' @keywords internal
-.guess_marginaltrends_arguments <- function(model,
-                                            trend = NULL,
-                                            by = NULL,
-                                            verbose = TRUE,
-                                            ...) {
+.guess_marginaltrends_arguments <- function(
+  model,
+  trend = NULL,
+  by = NULL,
+  verbose = TRUE,
+  ...
+) {
   # Gather info
+  dots <- list(...)
   model_data <- insight::get_data(model, verbose = FALSE)
   predictors <- intersect(
     colnames(model_data),
@@ -215,19 +222,33 @@ get_marginaltrends <- function(
   if (is.null(trend)) {
     trend <- predictors[sapply(model_data[predictors], is.numeric)][1]
     if (!length(trend) || is.na(trend)) {
-      insight::format_error("Model contains no numeric predictor. Please specify `trend`.")
+      insight::format_error(
+        "Model contains no numeric predictor. Please specify `trend`."
+      )
     }
     if (verbose) {
-      insight::format_alert(paste0("No numeric variable was specified for slope estimation. Selecting `trend = \"", trend, "\"`.")) # nolint
+      insight::format_alert(paste0(
+        "No numeric variable was specified for slope estimation. Selecting `trend = \"",
+        trend,
+        "\"`."
+      )) # nolint
     }
   }
 
   # check that we have only one predictor
   if (length(trend) > 1) {
     if (verbose) {
       insight::format_alert(paste0(
-        "More than one numeric variable was selected for slope estimation. Keeping only `", trend[1], "`. ", # nolint
-        "If you want to estimate the slope of `", trend[1], "` at different values of `", trend[2], "`, use `by=\"", trend[2], "\"` instead." # nolint
+        "More than one numeric variable was selected for slope estimation. Keeping only `",
+        trend[1],
+        "`. ", # nolint
+        "If you want to estimate the slope of `",
+        trend[1],
+        "` at different values of `",
+        trend[2],
+        "`, use `by=\"",
+        trend[2],
+        "\"` instead." # nolint
       ))
     }
     trend <- trend[1]
@@ -246,18 +267,28 @@ get_marginaltrends <- function(
   if (!is.null(by) && !is.null(range) && startsWith(by, trend)) {
     insight::format_error(
       paste0(
-        "To calculate average marginal effects over a range of `", trend, "` ",
-        "values, use `trend=\"", trend, "=seq(1, 3, 0.1)\"` (or similar) and omit `",
-        trend, "` from the `by` argument."
+        "To calculate average marginal effects over a range of `",
+        trend,
+        "` ",
+        "values, use `trend=\"",
+        trend,
+        "=seq(1, 3, 0.1)\"` (or similar) and omit `",
+        trend,
+        "` from the `by` argument."
       ),
       paste0(
-        "To get marginal effects at specific `", trend, "` values, use `trend=\"",
-        trend, "\"` along with `by=\"", trend, "=c(1, 3, 5)\"`."
+        "To get marginal effects at specific `",
+        trend,
+        "` values, use `trend=\"",
+        trend,
+        "\"` along with `by=\"",
+        trend,
+        "=c(1, 3, 5)\"`."
       )
     )
   }
 
-  list(trend = trend, range = range, by = NULL)
+  list(trend = trend, range = range, by = NULL, post_process = dots$post_process)
 }
 
 

R/utils.R

@@ -138,9 +138,11 @@
 
 #' @keywords internal
 #' @noRd
-.safe <- function(code, on_error = NULL) {
+.safe <- function(code, on_error = NULL, quietly = FALSE) {
   if (isTRUE(getOption("easystats_errors", FALSE)) && is.null(on_error)) {
     code
+  } else if (quietly) {
+    suppressWarnings(tryCatch(code, error = function(e) on_error))
   } else {
     tryCatch(code, error = function(e) on_error)
   }