fix

R/equivalence_test.R

@@ -340,6 +340,7 @@ equivalence_test.ggeffects <- function(x,
   l <- Map(
     function(ci_wide, ci_narrow) {
       .equivalence_test_numeric(
+        ci = ci,
         ci_wide,
         ci_narrow,
         range_rope = range,
@@ -433,11 +434,11 @@ equivalence_test.ggeffects <- function(x,
   l <- Map(
     function(ci_wide, ci_narrow) {
       .equivalence_test_numeric(
+        ci = ci,
         ci_wide,
         ci_narrow,
         range_rope = range,
         rule = rule,
-        ci = ci,
         verbose = verbose
       )
     }, conf_int, conf_int2
@@ -521,11 +522,11 @@ equivalence_test.ggeffects <- function(x,
     l <- Map(
       function(ci_wide, ci_narrow) {
         .equivalence_test_numeric(
+          ci = ci,
           ci_wide,
           ci_narrow,
           range_rope = range,
           rule = rule,
-          ci = ci,
           verbose = verbose
         )
       }, conf_int, conf_int2
@@ -543,7 +544,12 @@ equivalence_test.ggeffects <- function(x,
 
 
 #' @keywords internal
-.equivalence_test_numeric <- function(ci_wide, ci_narrow, range_rope, rule, ci = 0.95, verbose) {
+.equivalence_test_numeric <- function(ci = 0.95,
+                                      ci_wide,
+                                      ci_narrow,
+                                      range_rope,
+                                      rule,
+                                      verbose) {
   final_ci <- NULL
 
   # ==== HDI+ROPE decision rule, by Kruschke ====
@@ -594,7 +600,7 @@ equivalence_test.ggeffects <- function(x,
   data.frame(
     CI_low = final_ci[1],
     CI_high = final_ci[2],
-    SGPV = .rope_coverage(range_rope, final_ci),
+    SGPV = .rope_coverage(ci = ci, range_rope, ci_range = final_ci),
     ROPE_low = range_rope[1],
     ROPE_high = range_rope[2],
     ROPE_Equivalence = decision,
@@ -646,8 +652,8 @@ equivalence_test.ggeffects <- function(x,
 # same range / limits as the confidence interval, thus indeed representing a
 # normally distributed confidence interval. We then calculate the probability
 # mass of this interval that is inside the ROPE.
-.rope_coverage <- function(range_rope, ci_range) {
-  out <- .generate_posterior_from_ci(ci_range)
+.rope_coverage <- function(ci = 0.95, range_rope, ci_range) {
+  out <- .generate_posterior_from_ci(ci, ci_range)
   # compare: ci_range and range(out)
   # The SGPV refers to the proportion of the confidence interval inside the
   # full ROPE - thus, we set ci = 1 here
@@ -656,21 +662,43 @@ equivalence_test.ggeffects <- function(x,
 }
 
 
-.generate_posterior_from_ci <- function(ci_range, precision = 10000) {
+.generate_posterior_from_ci <- function(ci = 0.95, ci_range, precision = 10000) {
   # this function creates an approximate normal distribution that covers
   # the CI-range, i.e. we "simulate" a posterior distribution of a
   # frequentist CI
   diff_ci <- abs(diff(ci_range))
   bayestestR::distribution_normal(
     n = precision,
     mean = ci_range[2] - (diff_ci / 2),
-    # we divide the complete range by 2, the one-directional range for the SD
+    # we divide the complete range by 2, the one-directional range for the SD.
     # then, the range from mean value to lower/upper limit, for a normal
     # distribution is approximately 3.3 SD (3 SD cover 99.7% of the probability
-    # mass of the normal distribution). Thus, assuming that half of the ci_range
-    # refers to ~ 3.3 SD, we "normalize" the value (i.e. divide by 3.290525) to
-    # get the value for one SD, which we need to build the normal distribution.
-    sd = (diff_ci / 2) / 3.290525
+    # mass of the normal distribution, `1 - ((1 - pnorm(3)) * 2)`). Thus,
+    # assuming that half of the ci_range refers to ~ 3.3 SD, we "normalize" the
+    # value (i.e. divide by 3.3) to get the value for one SD, which we need
+    # to build the normal distribution. The SD itself varies by confidence level,
+    # therefore we have a multiplier based on the confidence level. I agree this
+    # looks *very* hacky, but it is tested against following code, which used
+    # this code to create a normal distribution with "full" coverage, based on
+    # the approximation of the SD related to the CI-level. From this normal-
+    # distribution, the CI-level % interval is drawn and the range of the
+    # simulated normal distribution equals the desired range.
+    # -------------------------------------------------------------------------
+    # m <- lm(mpg ~ gear + hp + wt + cyl + am, data = mtcars)
+    # ci <- 0.75
+    # mp <- model_parameters(m, ci = ci)
+    # ci_range <- c(mp$CI_low[2], mp$CI_high[2])
+    # diff_ci <- abs(diff(ci_range))
+    # out <- bayestestR::distribution_normal(
+    #   n = 10000,
+    #   mean = ci_range[2] - (diff_ci / 2),
+    #   sd = diff_ci / ((stats::qnorm((1+ci)/2) * (stats::qnorm(0.999975) / 2)))
+    # )
+    # # these to ranges are roughly the same
+    # ci(out, ci = ci)
+    # ci_range
+    # -------------------------------------------------------------------------
+    sd = diff_ci / ((stats::qnorm((1 + ci) / 2) * (stats::qnorm(0.999975) / 2)))
   )
 }
 

R/p_significance.R

@@ -63,7 +63,7 @@ p_significance.lm <- function(x, threshold = "default", ci = 0.95, verbose = TRU
   # distribution that covers the CI-range.
   posterior <- as.data.frame(lapply(seq_len(nrow(out)), function(i) {
     ci_range <- as.numeric(out[i, c("CI_low", "CI_high")])
-    .generate_posterior_from_ci(ci_range)
+    .generate_posterior_from_ci(ci, ci_range)
   }))
   colnames(posterior) <- out$Parameter
 

tests/testthat/_snaps/equivalence_test.md

@@ -10,8 +10,8 @@
       Parameter   |         90% CI |   SGPV | Equivalence |      p
       ------------------------------------------------------------
       (Intercept) | [26.52, 46.86] | < .001 |    Rejected | > .999
-      gear        | [-1.34,  2.07] | 0.648  |   Undecided | 0.578 
+      gear        | [-1.34,  2.07] | 0.480  |   Undecided | 0.578 
       wt          | [-4.47, -1.57] | < .001 |    Rejected | 0.996 
-      cyl         | [-1.94,  0.32] | 0.270  |   Undecided | 0.644 
+      cyl         | [-1.94,  0.32] | 0.350  |   Undecided | 0.644 
       hp          | [-0.05,  0.01] | > .999 |    Accepted | < .001
 

tests/testthat/_snaps/p_significance.md

@@ -8,8 +8,8 @@
       Parameter   |         95% CI |         ps
       -----------------------------------------
       (Intercept) | [24.44, 48.94] |  ps > .999
-      gear        | [-1.69,  2.41] | ps = 0.663
+      gear        | [-1.69,  2.41] | ps = 0.600
       wt          | [-4.77, -1.28] |  ps > .999
-      cyl         | [-2.17,  0.55] | ps = 0.958
+      cyl         | [-2.17,  0.55] | ps = 0.851
       hp          | [-0.05,  0.01] |  ps < .001
 

tests/testthat/test-equivalence_test.R

@@ -31,7 +31,7 @@ test_that("equivalence_test, unequal rope-range", {
     c("Rejected", "Accepted", "Undecided", "Rejected", "Accepted", "Undecided")
   )
 
-  # validate thet range of CI equals approximated  normal distribution
+  # validate that range of CI equals approximated normal distribution
   diff_ci <- abs(diff(c(rez$CI_low[3], rez$CI_high[3])))
   set.seed(123)
   out <- bayestestR::distribution_normal(
@@ -41,6 +41,13 @@ test_that("equivalence_test, unequal rope-range", {
   )
   expect_equal(range(out)[1], rez$CI_low[3], tolerance = 1e-4)
   expect_equal(range(out)[2], rez$CI_high[3], tolerance = 1e-4)
+  # need procedure for SGP here...
+  set.seed(123)
+  out <- bayestestR::distribution_normal(
+    n = 10000,
+    mean = rez$CI_high[3] - (diff_ci / 2),
+    sd = diff_ci / ((stats::qnorm((1 + 0.95) / 2) * (stats::qnorm(0.999975) / 2)))
+  )
   expect_equal(
     rez$SGPV[3],
     bayestestR::rope(out, range = c(-Inf, 0.5), ci = 1)$ROPE_Percentage,

tests/testthat/test-p_significance.R

@@ -12,9 +12,9 @@ test_that("p_significance", {
 
   set.seed(123)
   x <- p_significance(m, ci = 0.8)
-  expect_equal(x$ps, c(1, 0.7446, 1, 0.9964, 0), tolerance = 1e-4)
+  expect_equal(x$ps, c(1, 0.6025, 0.9997, 0.8561, 0), tolerance = 1e-4)
 
   set.seed(123)
   x <- p_significance(m, threshold = 0.5)
-  expect_equal(x$ps, c(1, 0.4128, 1, 0.7751, 0), tolerance = 1e-4)
+  expect_equal(x$ps, c(1, 0.4471, 0.998, 0.676, 0), tolerance = 1e-4)
 })