plot() for p_significance works with asymmetrical threshold (#361)

* plot() for p_significance works with asymmetrical threshold * minor fix * add more tests * update snapshots * tests only work with future updates * lintr * add test * fix * test * update snapshot * news * styler, lintr
easystats · Sep 5, 2024 · bdc2c62 · bdc2c62
1 parent 57ac150
commit bdc2c62
Show file tree

Hide file tree

Showing 17 changed files with 747 additions and 53 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -4,6 +4,9 @@
 
 - New `plot()` method for `performance::check_dag()`.
 
+- Minor improvements to `plot()` for methods `p_direction()` and `p_significance()`,
+  which also support forthcoming changes in the *parameters* package.
+
 ## Bug fixes
 
 - Fixed issue in `plot()` for `performance::check_model()` when package *qqplotr*

diff --git a/R/plot.check_predictions.R b/R/plot.check_predictions.R
@@ -536,7 +536,7 @@ plot.see_performance_pp_check <- function(x,
       size = ggplot2::guide_legend(reverse = TRUE)
     )
 
-  return(p)
+  p
 }
 
 

diff --git a/R/plot.equivalence_test.R b/R/plot.equivalence_test.R
@@ -110,7 +110,7 @@ plot.see_equivalence_test <- function(x,
   }
 
   # get labels
-  labels <- .clean_parameter_names(tmp$predictor, grid = !is.null(n_columns))
+  axis_labels <- .clean_parameter_names(tmp$predictor, grid = !is.null(n_columns))
 
   tmp <- .fix_facet_names(tmp)
 
@@ -126,11 +126,11 @@ plot.see_equivalence_test <- function(x,
 
   fill.color <- fill.color[sort(unique(match(x$ROPE_Equivalence, c("Accepted", "Rejected", "Undecided"))))]
 
-  add.args <- lapply(match.call(expand.dots = FALSE)$`...`, function(x) x)
+  add.args <- lapply(match.call(expand.dots = FALSE)[["..."]], function(x) x)
   if ("colors" %in% names(add.args)) fill.color <- eval(add.args[["colors"]])
   if ("x.title" %in% names(add.args)) x.title <- eval(add.args[["x.title"]])
   if ("legend.title" %in% names(add.args)) legend.title <- eval(add.args[["legend.title"]])
-  if ("labels" %in% names(add.args)) labels <- eval(add.args[["labels"]])
+  if ("labels" %in% names(add.args)) axis_labels <- eval(add.args[["labels"]])
 
   rope.line.alpha <- 1.25 * rope_alpha
   if (rope.line.alpha > 1) rope.line.alpha <- 1
@@ -170,7 +170,7 @@ plot.see_equivalence_test <- function(x,
     ) +
     scale_fill_manual(values = fill.color) +
     labs(x = x.title, y = NULL, fill = legend.title) +
-    scale_y_discrete(labels = labels) +
+    scale_y_discrete(labels = axis_labels) +
     theme(legend.position = "bottom")
 
   if (!is.null(n_columns)) {
@@ -193,10 +193,8 @@ plot.see_equivalence_test <- function(x,
         p <- p + facet_wrap(~Component, scales = "free", ncol = n_columns)
       }
     }
-  } else {
-    if (length(unique(tmp$HDI)) > 1L) {
-      p <- p + facet_wrap(~HDI, scales = "free", ncol = n_columns)
-    }
+  } else if (length(unique(tmp$HDI)) > 1L) {
+    p <- p + facet_wrap(~HDI, scales = "free", ncol = n_columns)
   }
 
   p
@@ -259,7 +257,7 @@ plot.see_equivalence_test_df <- function(x,
   tmp$predictor <- factor(tmp$predictor, levels = rev(unique(tmp$predictor)))
 
   # get labels
-  labels <- .clean_parameter_names(tmp$predictor, grid = !is.null(n_columns))
+  axis_labels <- .clean_parameter_names(tmp$predictor, grid = !is.null(n_columns))
 
   # check for user defined arguments
 
@@ -273,11 +271,11 @@ plot.see_equivalence_test_df <- function(x,
 
   fill.color <- fill.color[sort(unique(match(x$ROPE_Equivalence, c("Accepted", "Rejected", "Undecided"))))]
 
-  add.args <- lapply(match.call(expand.dots = FALSE)$`...`, function(x) x)
+  add.args <- lapply(match.call(expand.dots = FALSE)[["..."]], function(x) x)
   if ("colors" %in% names(add.args)) fill.color <- eval(add.args[["colors"]])
   if ("x.title" %in% names(add.args)) x.title <- eval(add.args[["x.title"]])
   if ("legend.title" %in% names(add.args)) legend.title <- eval(add.args[["legend.title"]])
-  if ("labels" %in% names(add.args)) labels <- eval(add.args[["labels"]])
+  if ("labels" %in% names(add.args)) axis_labels <- eval(add.args[["labels"]])
 
   rope.line.alpha <- 1.25 * rope_alpha
 
@@ -317,7 +315,7 @@ plot.see_equivalence_test_df <- function(x,
     ) +
     scale_fill_manual(values = fill.color) +
     labs(x = x.title, y = NULL, fill = legend.title) +
-    scale_y_discrete(labels = labels) +
+    scale_y_discrete(labels = axis_labels) +
     theme(legend.position = "bottom")
 
   if (length(unique(tmp$HDI)) > 1L) {
@@ -390,11 +388,10 @@ plot.see_equivalence_test_lm <- function(x,
 
   fill.color <- fill.color[sort(unique(match(x$ROPE_Equivalence, c("Accepted", "Rejected", "Undecided"))))]
 
-  add.args <- lapply(match.call(expand.dots = FALSE)$`...`, function(x) x)
+  add.args <- lapply(match.call(expand.dots = FALSE)[["..."]], function(x) x)
   if ("colors" %in% names(add.args)) fill.color <- eval(add.args[["colors"]])
   if ("x.title" %in% names(add.args)) x.title <- eval(add.args[["x.title"]])
   if ("legend.title" %in% names(add.args)) legend.title <- eval(add.args[["legend.title"]])
-  if ("labels" %in% names(add.args)) labels <- eval(add.args[["labels"]])
 
   rope.line.alpha <- 1.25 * rope_alpha
   if (rope.line.alpha > 1) rope.line.alpha <- 1

diff --git a/R/plot.p_direction.R b/R/plot.p_direction.R
@@ -27,17 +27,17 @@ data_plot.p_direction <- function(x, data = NULL, show_intercept = FALSE, ...) {
     data <- data[, x$Parameter, drop = FALSE]
     dataplot <- data.frame()
     for (i in names(data)) {
-      if (!is.null(params)) {
+      if (!is.null(params) && all(c("Effects", "Component") %in% colnames(params))) {
         dataplot <- rbind(
           dataplot,
           cbind(
-            .compute_densities_pd(data[[i]], name = i),
-            "Effects" = params$Effects[params$Parameter == i],
-            "Component" = params$Component[params$Parameter == i]
+            .compute_densities_pd(data[[i]], name = i, null = attr(x, "null")),
+            Effects = params$Effects[params$Parameter == i],
+            Component = params$Component[params$Parameter == i]
           )
         )
       } else {
-        dataplot <- rbind(dataplot, .compute_densities_pd(data[[i]], name = i))
+        dataplot <- rbind(dataplot, .compute_densities_pd(data[[i]], name = i, null = attr(x, "null")))
       }
     }
 
@@ -60,7 +60,7 @@ data_plot.p_direction <- function(x, data = NULL, show_intercept = FALSE, ...) {
     }
   } else {
     levels_order <- NULL
-    dataplot <- .compute_densities_pd(data[, 1], name = "Posterior")
+    dataplot <- .compute_densities_pd(data[, 1], name = "Posterior", null = attr(x, "null"))
   }
 
   dataplot <- do.call(
@@ -70,7 +70,7 @@ data_plot.p_direction <- function(x, data = NULL, show_intercept = FALSE, ...) {
       list(dataplot$y, dataplot$fill),
       function(df) {
         df$n <- nrow(df)
-        return(df)
+        df
       }
     )
   )
@@ -81,7 +81,7 @@ data_plot.p_direction <- function(x, data = NULL, show_intercept = FALSE, ...) {
       dataplot$y,
       function(df) {
         df$prop <- df$n / nrow(df)
-        return(df)
+        df
       }
     )
   )
@@ -108,10 +108,10 @@ data_plot.p_direction <- function(x, data = NULL, show_intercept = FALSE, ...) {
   dataplot <- .fix_facet_names(dataplot)
 
   attr(dataplot, "info") <- list(
-    "xlab" = "Possible parameter values",
-    "ylab" = ylab,
-    "legend_fill" = "Effect direction",
-    "title" = "Probability of Direction"
+    xlab = "Possible parameter values",
+    ylab = ylab,
+    legend_fill = "Effect direction",
+    title = "Probability of Direction"
   )
 
   class(dataplot) <- c("data_plot", "see_p_direction", class(dataplot))
@@ -121,11 +121,14 @@ data_plot.p_direction <- function(x, data = NULL, show_intercept = FALSE, ...) {
 
 
 #' @keywords internal
-.compute_densities_pd <- function(x, name = "Y") {
+.compute_densities_pd <- function(x, name = "Y", null = 0) {
   out <- .as.data.frame_density(
     stats::density(x)
   )
-  out$fill <- ifelse(out$x < 0, "Negative", "Positive")
+  if (is.null(null)) {
+    null <- 0
+  }
+  out$fill <- ifelse(out$x < null, "Negative", "Positive")
   out$height <- as.vector(
     (out$y - min(out$y, na.rm = TRUE)) /
       diff(range(out$y, na.rm = TRUE), na.rm = TRUE)
@@ -182,7 +185,7 @@ plot.see_p_direction <- function(x,
   params <- unique(x$y)
 
   # get labels
-  labels <- .clean_parameter_names(x$y, grid = !is.null(n_columns))
+  axis_labels <- .clean_parameter_names(x$y, grid = !is.null(n_columns))
 
   insight::check_if_installed("ggridges")
 
@@ -216,7 +219,7 @@ plot.see_p_direction <- function(x,
   if (length(unique(x$y)) == 1 && is.numeric(x$y)) {
     p <- p + scale_y_continuous(breaks = NULL, labels = NULL)
   } else {
-    p <- p + scale_y_discrete(labels = labels)
+    p <- p + scale_y_discrete(labels = axis_labels)
   }
 
 

diff --git a/R/plot.p_significance.R b/R/plot.p_significance.R
@@ -12,7 +12,6 @@ data_plot.p_significance <- function(x,
 
   if (inherits(data, "emmGrid")) {
     insight::check_if_installed("emmeans")
-
     data <- as.data.frame(as.matrix(emmeans::as.mcmc.emmGrid(data, names = FALSE)))
   } else if (inherits(data, c("stanreg", "brmsfit"))) {
     params <- insight::clean_parameters(data)
@@ -32,19 +31,19 @@ data_plot.p_significance <- function(x,
     data <- data[, x$Parameter, drop = FALSE]
     dataplot <- data.frame()
     for (i in names(data)) {
-      if (!is.null(params)) {
+      if (is.null(params) || !all(c("Effects", "Component") %in% colnames(params))) {
         dataplot <- rbind(
           dataplot,
-          cbind(
-            .compute_densities_ps(data[[i]], name = i, threshold = attr(x, "threshold")),
-            "Effects" = params$Effects[params$Parameter == i],
-            "Component" = params$Component[params$Parameter == i]
-          )
+          .compute_densities_ps(data[[i]], name = i, threshold = attr(x, "threshold"))
         )
       } else {
         dataplot <- rbind(
           dataplot,
-          .compute_densities_ps(data[[i]], name = i, threshold = attr(x, "threshold"))
+          cbind(
+            .compute_densities_ps(data[[i]], name = i, threshold = attr(x, "threshold")),
+            Effects = params$Effects[params$Parameter == i],
+            Component = params$Component[params$Parameter == i]
+          )
         )
       }
     }
@@ -68,7 +67,7 @@ data_plot.p_significance <- function(x,
     }
   } else {
     levels_order <- NULL
-    dataplot <- .compute_densities_pd(data[, 1], name = "Posterior")
+    dataplot <- .compute_densities_ps(data[, 1], name = "Posterior", threshold = attr(x, "threshold"))
   }
 
   dataplot <- do.call(
@@ -78,7 +77,7 @@ data_plot.p_significance <- function(x,
       list(dataplot$y, dataplot$fill),
       function(df) {
         df$n <- nrow(df)
-        return(df)
+        df
       }
     )
   )
@@ -89,7 +88,7 @@ data_plot.p_significance <- function(x,
       dataplot$y,
       function(df) {
         df$prop <- df$n / nrow(df)
-        return(df)
+        df
       }
     )
   )
@@ -116,10 +115,10 @@ data_plot.p_significance <- function(x,
   dataplot <- .fix_facet_names(dataplot)
 
   attr(dataplot, "info") <- list(
-    "xlab" = "Possible parameter values",
-    "ylab" = ylab,
-    "legend_fill" = "Probability",
-    "title" = "Practical Significance"
+    xlab = "Possible parameter values",
+    ylab = ylab,
+    legend_fill = "Probability",
+    title = "Practical Significance"
   )
 
   class(dataplot) <- c("data_plot", "see_p_significance", class(dataplot))
@@ -132,18 +131,35 @@ data_plot.p_significance <- function(x,
 .compute_densities_ps <- function(x, name = "Y", threshold = 0) {
   out <- .as.data.frame_density(stats::density(x))
 
-  fifty_cents <- sum(out$y[out$x > threshold]) > (sum(out$y) / 2)
+  # sanity check
+  if (is.null(threshold)) {
+    threshold <- 0
+  }
+
+  # make sure we have a vector of length 2
+  if (length(threshold) == 1) {
+    threshold <- c(-1 * threshold, threshold)
+  }
+
+  # find out the probability mass larger or lower than the ROPE (outside)
+  p_mass_ht_rope <- sum(out$y[out$x > threshold[2]])
+  p_mass_lt_rope <- sum(out$y[out$x < threshold[1]])
+
+  # find out whether probability mass "above" ROPE is larger than the probability
+  # mass that is on the left (negative) side of the ROPE
+  fifty_cents <- p_mass_ht_rope > p_mass_lt_rope
 
   out$fill <- "Less Probable"
-  out$fill[abs(out$x) < threshold] <- "ROPE"
-  out$fill[(out$x > threshold)] <- ifelse(fifty_cents, "Significant", "Less Probable")
-  out$fill[out$x < (-1 * threshold)] <- ifelse(fifty_cents, "Less Probable", "Significant")
+  out$fill[out$x > threshold[1] & out$x < threshold[2]] <- "ROPE"
+  out$fill[out$x > threshold[2]] <- ifelse(fifty_cents, "Significant", "Less Probable")
+  out$fill[out$x < threshold[1]] <- ifelse(fifty_cents, "Less Probable", "Significant")
 
   out$height <- out$y
   out$y <- name
 
   # normalize
-  out$height <- as.vector((out$height - min(out$height, na.rm = TRUE)) / diff(range(out$height, na.rm = TRUE), na.rm = TRUE))
+  range_diff <- diff(range(out$height, na.rm = TRUE), na.rm = TRUE)
+  out$height <- as.vector((out$height - min(out$height, na.rm = TRUE)) / range_diff)
   out
 }
 
@@ -194,7 +210,7 @@ plot.see_p_significance <- function(x,
   params <- unique(x$y)
 
   # get labels
-  labels <- .clean_parameter_names(x$y, grid = !is.null(n_columns))
+  axis_labels <- .clean_parameter_names(x$y, grid = !is.null(n_columns))
 
   insight::check_if_installed("ggridges")
 
@@ -235,7 +251,7 @@ plot.see_p_significance <- function(x,
   if (length(unique(x$y)) == 1L && is.numeric(x$y)) {
     p <- p + scale_y_continuous(breaks = NULL, labels = NULL)
   } else {
-    p <- p + scale_y_discrete(labels = labels)
+    p <- p + scale_y_discrete(labels = axis_labels)
   }
 
 

diff --git a/R/utils.R b/R/utils.R
@@ -42,6 +42,8 @@
 
   # clean parameters names
   params <- gsub("(b_|bs_|bsp_|bcs_)(.*)", "\\2", params, perl = TRUE)
+  params <- gsub("^cond_(.*)", "\\1 (Conditional)", params, perl = TRUE)
+  params <- gsub("(.*)_cond$", "\\1 (Conditional)", params, perl = TRUE)
   params <- gsub("^zi_(.*)", "\\1 (Zero-Inflated)", params, perl = TRUE)
   params <- gsub("(.*)_zi$", "\\1 (Zero-Inflated)", params, perl = TRUE)
   params <- gsub("(.*)_disp$", "\\1 (Dispersion)", params, perl = TRUE)