Draft marginalcontrasts

R/estimate_contrasts.R

@@ -7,9 +7,9 @@
 #' @inheritParams estimate_means
 #' @inheritParams get_emcontrasts
 #' @param p_adjust The p-values adjustment method for frequentist multiple
-#'   comparisons. Can be one of `"holm"` (default), `"tukey"`, `"hochberg"`,
-#'   `"hommel"`, `"bonferroni"`, `"BH"`, `"BY"`, `"fdr"` or `"none"`. See the
-#'   p-value adjustment section in the `emmeans::test` documentation.
+#' comparisons. Can be one of `"holm"` (default), `"tukey"`, `"hochberg"`,
+#' `"hommel"`, `"bonferroni"`, `"BH"`, `"BY"`, `"fdr"` or `"none"`. See the
+#' p-value adjustment section in the `emmeans::test` documentation.
 #'
 #' @inherit estimate_slopes details
 #'
@@ -78,19 +78,67 @@
                                ci = 0.95,
                                p_adjust = "holm",
                                method = "pairwise",
+                               backend = "emmeans",
                                ...) {
-  # Run emmeans
-  estimated <- get_emcontrasts(model,
-    contrast = contrast,
-    by = by,
-    transform = transform,
-    method = method,
-    adjust = p_adjust,
-    ...
+  if (backend == "emmeans") {
+    # Emmeans ------------------------------------------------------------------
+    estimated <- get_emcontrasts(model,
+      contrast = contrast,
+      by = by,
+      transform = transform,
+      method = method,
+      adjust = p_adjust,
+      ...
+    )
+    out <- .format_emmeans_contrasts(model, estimated, ci, transform, p_adjust, ...)
+    info <- attributes(estimated)
+  } else {
+    # Marginalmeans ------------------------------------------------------------
+    estimated <- get_marginalcontrasts(model,
+      contrast = contrast,
+      by = by,
+      transform = transform,
+      method = method,
+      p_adjust = p_adjust,
+      ci = ci,
+      ...
+    )
+    out <- .format_marginaleffects_contrasts(model, estimated, p_adjust, method, ...)
+    ## TODO: needs to be fixed
+    info <- list(contrast = contrast, by = by)
+  }
+
+
+  # Table formatting
+  attr(out, "table_title") <- c("Marginal Contrasts Analysis", "blue")
+  attr(out, "table_footer") <- .estimate_means_footer(
+    out,
+    info$contrast,
+    type = "contrasts",
+    p_adjust = p_adjust
   )
 
-  info <- attributes(estimated)
+  # Add attributes
+  attr(out, "model") <- model
+  attr(out, "response") <- insight::find_response(model)
+  attr(out, "ci") <- ci
+  attr(out, "transform") <- transform
+  attr(out, "at") <- info$by
+  attr(out, "by") <- info$by
+  attr(out, "contrast") <- info$contrast
+  attr(out, "p_adjust") <- p_adjust
+
+  # Output
+  class(out) <- c("estimate_contrasts", "see_estimate_contrasts", class(out))
+  out
+}
+
+
+
+# Table formatting emmeans ----------------------------------------------------
 
+
+.format_emmeans_contrasts <- function(model, estimated, ci, transform, p_adjust, ...) {
   # Summarize and clean
   if (insight::model_info(model)$is_bayesian) {
     out <- cbind(estimated@grid, bayestestR::describe_posterior(estimated, ci = ci, verbose = FALSE, ...))
@@ -120,30 +168,29 @@
 
   # Merge levels and rest
   out$contrast <- NULL
-  out <- cbind(level_cols, out)
+  cbind(level_cols, out)
+}
 
 
-  # Table formatting
-  attr(out, "table_title") <- c("Marginal Contrasts Analysis", "blue")
-  attr(out, "table_footer") <- .estimate_means_footer(
-    out,
-    info$contrast,
-    type = "contrasts",
-    p_adjust = p_adjust
-  )
 
-  # Add attributes
-  attr(out, "model") <- model
-  attr(out, "response") <- insight::find_response(model)
-  attr(out, "ci") <- ci
-  attr(out, "transform") <- transform
-  attr(out, "at") <- info$by
-  attr(out, "by") <- info$by
-  attr(out, "contrast") <- info$contrast
-  attr(out, "p_adjust") <- p_adjust
+# Table formatting marginal effects -------------------------------------------
 
 
-  # Output
-  class(out) <- c("estimate_contrasts", "see_estimate_contrasts", class(out))
-  out
+.format_marginaleffects_contrasts <- function(model, estimated, p_adjust, method, ...) {
+  groups <- attributes(estimated)$by
+  contrast <- attributes(estimated)$contrast
+  focal_terms <- attributes(estimated)$focal_terms
+
+  estimated <- .p_adjust(model, estimated, p_adjust, ...)
+
+  valid_methods <- c(
+    "pairwise", "reference", "sequential", "meandev","meanotherdev",
+    "revpairwise", "revreference", "revsequential"
+  )
+  ## TODO: split Parameter column into levels indicated in "contrast", and filter by "by"
+  if (!is.null(method) && is.character(method) && method %in% valid_methods) {
+
+  }
+
+  estimated
 }

R/get_marginalcontrasts.R

@@ -1,27 +1,86 @@
 #' @rdname get_marginalmeans
 #'
 #' @param method Contrast method.
+#' @inheritParams get_emcontrasts
 #' @export
 get_marginalcontrasts <- function(model,
+                                  contrast = NULL,
                                   by = NULL,
-                                  ci = 0.95,
+                                  transform = "none",
                                   method = "pairwise",
+                                  ci = 0.95,
                                   ...) {
+  # check if available
+  insight::check_if_installed("marginaleffects")
+
   out <- estimate_means(
     model = model,
-    by = by,
+    by = c(contrast, by),
     ci = ci,
     hypothesis = method,
+    transform = "response",
     backend = "marginaleffects",
     ...
   )
 
-  attr(out, "table_title") <- c("Marginal Contrasts Analysis", "blue")
-  attr(out, "table_footer") <- .estimate_means_footer(out, by = by, type = "contrasts", p_adjust = NULL)
-
+  attr(out, "contrast") <- contrast
   out
 }
 
 #' @rdname get_marginalmeans
 #' @export
 model_marginalcontrasts <- get_marginalcontrasts
+
+
+# p-value adjustment -------------------
+
+.p_adjust <- function(model, params, p_adjust, ...) {
+  # extract information
+  datagrid <- attributes(params)$datagrid
+  focal <- attributes(params)$focal_terms
+  statistic <- insight::get_statistic(model)$Statistic
+  df <- insight::get_df(model)
+  verbose <- isTRUE(list(...)$verbose)
+
+  # exit on NULL, or if no p-adjustment requested
+  if (is.null(p_adjust) || identical(p_adjust, "none")) {
+    return(params)
+  }
+
+  all_methods <- c(tolower(stats::p.adjust.methods), "tukey", "sidak")
+
+  # needed for rank adjustment
+  focal_terms <- datagrid[focal]
+  rank_adjust <- prod(vapply(focal_terms, insight::n_unique, numeric(1)))
+
+  # only proceed if valid argument-value
+  if (tolower(p_adjust) %in% all_methods) {
+    if (tolower(p_adjust) %in% tolower(stats::p.adjust.methods)) {
+      # base R adjustments
+      params[["p"]] <- stats::p.adjust(params[["p"]], method = p_adjust)
+    } else if (tolower(p_adjust) == "tukey") {
+      if (!is.null(statistic)) {
+        # tukey adjustment
+        params[["p"]] <- suppressWarnings(stats::ptukey(
+          sqrt(2) * abs(statistic),
+          rank_adjust,
+          df,
+          lower.tail = FALSE
+        ))
+        # for specific contrasts, ptukey might fail, and the tukey-adjustement
+        # could just be simple p-value calculation
+        if (all(is.na(params[["p"]]))) {
+          params[["p"]] <- 2 * stats::pt(abs(statistic), df = df, lower.tail = FALSE)
+        }
+      } else if (verbose) {
+        insight::format_alert("No test-statistic found. P-values were not adjusted.")
+      }
+    } else if (tolower(p_adjust) == "sidak") {
+      # sidak adjustment
+      params[["p"]] <- 1 - (1 - params[["p"]])^rank_adjust
+    }
+  } else if (verbose) {
+    insight::format_alert(paste0("`p_adjust` must be one of ", toString(all_methods)))
+  }
+  params
+}

R/get_marginalmeans.R

@@ -102,6 +102,7 @@ get_marginalmeans <- function(model,
   attr(means, "at") <- my_args$by
   attr(means, "by") <- my_args$by
   attr(means, "focal_terms") <- at_specs$varname
+  attr(means, "datagrid") <- datagrid
   means
 }
 
@@ -119,19 +120,30 @@ model_marginalmeans <- get_marginalmeans
   model_data <- insight::get_data(model)
   info <- insight::model_info(model, verbose = FALSE)
   non_focal <- setdiff(colnames(model_data), attr(means, "focal_terms"))
+  is_contrast_analysis <- !is.null(list(...)$hypothesis)
 
-  # estimate name
-  if (!identical(transform, "none") && (info$is_binomial || info$is_bernoulli)) {
-    estimate_name <- "Probability"
+  # do we have contrasts? For contrasts, we want to keep p-values
+  if (is_contrast_analysis) {
+    remove_column <- "SE"
+    estimate_name <- "Difference"
   } else {
-    estimate_name <- "Mean"
+    remove_column <- "p"
+    # estimate name
+    if (!identical(transform, "none") && (info$is_binomial || info$is_bernoulli)) {
+      estimate_name <- "Probability"
+    } else {
+      estimate_name <- "Mean"
+    }
   }
 
   # Format
   params <- suppressWarnings(parameters::model_parameters(means, verbose = FALSE))
   params <- datawizard::data_relocate(params, c("Predicted", "SE", "CI_low", "CI_high"), after = -1, verbose = FALSE) # nolint
+  # move p to the end
+  params <- datawizard::data_relocate(params, "p", after = -1, verbose = FALSE)
   params <- datawizard::data_rename(params, "Predicted", estimate_name)
-  params <- datawizard::data_remove(params, c("p", "Statistic", "s.value", "S", "CI", "df", "rowid_dedup", non_focal), verbose = FALSE) # nolint
+  # remove redundant columns
+  params <- datawizard::data_remove(params, c(remove_column, "Statistic", "s.value", "S", "CI", "df", "rowid_dedup", non_focal), verbose = FALSE) # nolint
   params <- datawizard::data_restoretype(params, model_data)
 
   # Store info