Skip to content

Commit

Permalink
equi_test
Browse files Browse the repository at this point in the history
  • Loading branch information
strengejacke committed Sep 17, 2024
1 parent 9ed77e1 commit 983a719
Show file tree
Hide file tree
Showing 5 changed files with 83 additions and 20 deletions.
46 changes: 37 additions & 9 deletions R/equivalence_test.R
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,8 @@
#' \donttest{
#' model <- rstanarm::stan_glm(mpg ~ wt + cyl, data = mtcars)
#' equivalence_test(model)
#' # multiple ROPE ranges - asymmetric, symmetric, default
#' equivalence_test(model, range = list(c(10, 40), c(-5, -4), "default"))
#'
#' # plot result
#' test <- equivalence_test(model)
Expand Down Expand Up @@ -244,7 +246,7 @@ equivalence_test.BFBayesFactor <- function(x, range = "default", ci = 0.95, verb
verbose = TRUE) {
if (all(range == "default")) {
range <- rope_range(x, verbose = verbose)
} else if (!all(is.numeric(range)) || length(range) != 2L) {
} else if (!is.list(range) && (!all(is.numeric(range)) || length(range) != 2L)) {
insight::format_error("`range` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).")
}

Expand All @@ -257,14 +259,40 @@ equivalence_test.BFBayesFactor <- function(x, range = "default", ci = 0.95, verb
verbose = verbose
)

l <- sapply(
params,
equivalence_test,
range = range,
ci = ci,
verbose = verbose,
simplify = FALSE
)
if (is.list(range)) {
if (length(range) != ncol(params)) {
insight::format_error("Length of `range` (i.e. number of ROPE limits) should match the number of parameters.")
}
# check if list of values contains only valid values
checks <- vapply(range, function(r) {
!all(r == "default") || !all(is.numeric(r)) || length(r) != 2
}, logical(1))
if (!all(checks)) {
insight::format_error("`range` should be 'default' or a vector of 2 numeric values (e.g., c(-0.1, 0.1)).")
}
l <- mapply(
function(p, r) {
equivalence_test(
p,
range = r,
ci = ci,
verbose = verbose
)
},
params,
range,
SIMPLIFY = FALSE
)
} else {
l <- sapply(
params,
equivalence_test,
range = range,
ci = ci,
verbose = verbose,
simplify = FALSE
)
}

dat <- do.call(rbind, l)
out <- data.frame(
Expand Down
20 changes: 11 additions & 9 deletions R/print.equivalence_test.R
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,10 @@
print.equivalence_test <- function(x, digits = 2, ...) {
orig_x <- x
insight::print_color("# Test for Practical Equivalence\n\n", "blue")
cat(sprintf(" ROPE: [%.*f %.*f]\n\n", digits, x$ROPE_low[1], digits, x$ROPE_high[1]))
# print ROPE limits, if we just have one set of ROPE values
if (insight::n_unique(x$ROPE_low) == 1) {
cat(sprintf(" ROPE: [%.*f %.*f]\n\n", digits, x$ROPE_low[1], digits, x$ROPE_high[1]))
}

# fix "sd" pattern
model <- .retrieve_model(x)
Expand All @@ -23,22 +26,21 @@ print.equivalence_test <- function(x, digits = 2, ...) {
}
}

# find the longest HDI-value, so we can align the brackets in the ouput
x$HDI_low <- sprintf("%.*f", digits, x$HDI_low)
x$HDI_high <- sprintf("%.*f", digits, x$HDI_high)

maxlen_low <- max(nchar(x$HDI_low))
maxlen_high <- max(nchar(x$HDI_high))

x$ROPE_Percentage <- sprintf("%.*f %%", digits, x$ROPE_Percentage * 100)
x$HDI <- sprintf("[%*s %*s]", maxlen_low, x$HDI_low, maxlen_high, x$HDI_high)
x$HDI <- insight::format_ci(x$HDI_low, x$HDI_high, ci = NULL, digits = digits)

ci <- unique(x$CI)
keep.columns <- c(
attr(x, "idvars"), "Parameter", "Effects", "Component",
"ROPE_Equivalence", "ROPE_Percentage", "CI", "HDI"
)

# keep ROPE columns for multiple ROPE values
if (insight::n_unique(x$ROPE_low) > 1) {
keep.columns <- c(keep.columns, "ROPE")
x$ROPE <- insight::format_ci(x$ROPE_low, x$ROPE_high, ci = NULL, digits = digits)
}

x <- x[, intersect(keep.columns, colnames(x))]

colnames(x)[which(colnames(x) == "ROPE_Equivalence")] <- "H0"
Expand Down
2 changes: 2 additions & 0 deletions man/equivalence_test.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

32 changes: 32 additions & 0 deletions tests/testthat/test-equivalence_test.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,32 @@
test_that("equivalence test", {
skip_if_offline()
skip_if_not_or_load_if_installed("rstanarm")
m <- insight::download_model("stanreg_merMod_5")

out <- equivalence_test(m, verbose = FALSE)
expect_snapshot(print(out))

out <- equivalence_test(
m,
range = list(c(-1, 1), "default", c(0, 2), c(-2, 0), "default"),
verbose = FALSE
)
expect_snapshot(print(out))

expect_error(
equivalence_test(
m,
range = list(c(-1, 1), "default", c(0, 2), c(-2, 0)),
verbose = FALSE
),
regex = "Length of"
)
expect_error(
equivalence_test(
m,
range = list(c(-1, 1), "default", c(0, 2), c(-2, 0), "a"),
verbose = FALSE
),
regex = "should be 'default'"
)
})
3 changes: 1 addition & 2 deletions tests/testthat/test-p_rope.R
Original file line number Diff line number Diff line change
@@ -1,8 +1,7 @@
test_that("rope", {
test_that("p_rope", {
skip_if_offline()
skip_if_not_or_load_if_installed("rstanarm")
m <- insight::download_model("stanreg_merMod_5")
p <- insight::get_parameters(m, effects = "all")
expect_equal(
p_rope(as.data.frame(m)[2:4], range = list(c(0, 40), "default", c(-1, 0.8)))$p_ROPE,
c(0.598, 0.002, 0.396),
Expand Down

0 comments on commit 983a719

Please sign in to comment.