Skip to content

Commit

Permalink
plot() for p_significance works with asymmetrical threshold (#361)
Browse files Browse the repository at this point in the history
* 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
  • Loading branch information
strengejacke authored Sep 5, 2024
1 parent 57ac150 commit bdc2c62
Show file tree
Hide file tree
Showing 17 changed files with 747 additions and 53 deletions.
3 changes: 3 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -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*
Expand Down
2 changes: 1 addition & 1 deletion R/plot.check_predictions.R
Original file line number Diff line number Diff line change
Expand Up @@ -536,7 +536,7 @@ plot.see_performance_pp_check <- function(x,
size = ggplot2::guide_legend(reverse = TRUE)
)

return(p)
p
}


Expand Down
25 changes: 11 additions & 14 deletions R/plot.equivalence_test.R
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand All @@ -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
Expand Down Expand Up @@ -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)) {
Expand All @@ -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
Expand Down Expand Up @@ -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

Expand All @@ -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

Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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
Expand Down
35 changes: 19 additions & 16 deletions R/plot.p_direction.R
Original file line number Diff line number Diff line change
Expand Up @@ -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")))
}
}

Expand All @@ -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(
Expand All @@ -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
}
)
)
Expand All @@ -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
}
)
)
Expand All @@ -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))
Expand All @@ -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)
Expand Down Expand Up @@ -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")

Expand Down Expand Up @@ -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)
}


Expand Down
60 changes: 38 additions & 22 deletions R/plot.p_significance.R
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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]
)
)
}
}
Expand All @@ -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(
Expand All @@ -78,7 +77,7 @@ data_plot.p_significance <- function(x,
list(dataplot$y, dataplot$fill),
function(df) {
df$n <- nrow(df)
return(df)
df
}
)
)
Expand All @@ -89,7 +88,7 @@ data_plot.p_significance <- function(x,
dataplot$y,
function(df) {
df$prop <- df$n / nrow(df)
return(df)
df
}
)
)
Expand All @@ -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))
Expand All @@ -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
}

Expand Down Expand Up @@ -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")

Expand Down Expand Up @@ -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)
}


Expand Down
2 changes: 2 additions & 0 deletions R/utils.R
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Loading

0 comments on commit bdc2c62

Please sign in to comment.