Add vcov option to docs and examples

easystats · Sep 16, 2024 · 35915a7 · 35915a7
1 parent c4f5abe
commit 35915a7
Show file tree

Hide file tree

Showing 17 changed files with 144 additions and 28 deletions.
diff --git a/R/2_ci.R b/R/2_ci.R
@@ -24,13 +24,25 @@
 #' @param iterations The number of bootstrap replicates. Only applies to models
 #'   of class `merMod` when `method=boot`.
 #' @param verbose Toggle warnings and messages.
-#' @param ... Additional arguments
+#' @param ... Additional arguments passed down to the underlying functions.
+#' E.g., arguments like `vcov` or `vcov_args` can be used to compute confidence
+#' intervals using a specific variance-covariance matrix for the standard
+#' errors.
 #'
 #' @return A data frame containing the CI bounds.
 #'
 #' @inheritSection model_parameters Confidence intervals and approximation of degrees of freedom
 #'
-#' @examplesIf require("glmmTMB")
+#' @examplesIf require("glmmTMB") && requireNamespace("sandwich")
+#' data(qol_cancer)
+#' model <- lm(QoL ~ time + age + education, data = qol_cancer)
+#'
+#' # regular confidence intervals
+#' ci(model)
+#'
+#' # using heteroscedasticity-robust standard errors
+#' ci(model, vcov = "HC3")
+#'
 #' \donttest{
 #' library(parameters)
 #' data(Salamanders, package = "glmmTMB")

diff --git a/R/equivalence_test.R b/R/equivalence_test.R
@@ -19,7 +19,10 @@ bayestestR::equivalence_test
 #' See [`?ggeffects::test_predictions`](https://strengejacke.github.io/ggeffects/reference/test_predictions.html)
 #' for details.
 #' @param verbose Toggle warnings and messages.
-#' @param ... Arguments passed to or from other methods.
+#' @param ... Arguments passed to or from other methods, e.g. `ci()`. Arguments
+#' like `vcov` or `vcov_args` can be used to compute confidence intervals or
+#' p-values using a specific variance-covariance matrix for the standard
+#' errors..
 #' @inheritParams model_parameters.merMod
 #' @inheritParams p_value
 #'
@@ -218,13 +221,16 @@ bayestestR::equivalence_test
 #'     Synthese 200, 89 (2022). \doi{10.1007/s11229-022-03560-x}
 #'
 #' @return A data frame.
-#' @examples
+#' @examplesIf requireNamespace("sandwich")
 #' data(qol_cancer)
 #' model <- lm(QoL ~ time + age + education, data = qol_cancer)
 #'
 #' # default rule
 #' equivalence_test(model)
 #'
+#' # using heteroscedasticity-robust standard errors
+#' equivalence_test(model, vcov = "HC3")
+#'
 #' # conditional equivalence test
 #' equivalence_test(model, rule = "cet")
 #'
@@ -504,14 +510,14 @@ equivalence_test.ggeffects <- function(x,
 
   # ==== requested confidence intervals ====
 
-  params <- conf_int <- .ci_generic(x, ci = ci)
+  params <- conf_int <- .ci_generic(x, ci = ci, ...)
   conf_int <- as.data.frame(t(conf_int[, c("CI_low", "CI_high")]))
 
 
   # ==== the "narrower" intervals (1-2*alpha) for CET-rules. ====
 
   alpha <- 1 - ci
-  conf_int2 <- .ci_generic(x, ci = (ci - alpha))
+  conf_int2 <- .ci_generic(x, ci = (ci - alpha), ...)
   conf_int2 <- as.data.frame(t(conf_int2[, c("CI_low", "CI_high")]))
 
 
@@ -544,7 +550,7 @@ equivalence_test.ggeffects <- function(x,
 
   # ==== (adjusted) p-values for tests ====
 
-  out$p <- .add_p_to_equitest(x, ci, range)
+  out$p <- .add_p_to_equitest(x, ci, range, ...)
 
   attr(out, "rope") <- range
   out
@@ -786,7 +792,7 @@ equivalence_test.ggeffects <- function(x,
 }
 
 
-.add_p_to_equitest <- function(model, ci, range) {
+.add_p_to_equitest <- function(model, ci, range, ...) {
   tryCatch(
     {
       params <- insight::get_parameters(model)
@@ -798,7 +804,7 @@ equivalence_test.ggeffects <- function(x,
       params$mu <- params$Estimate * -1
 
       # se
-      se <- standard_error(model)
+      se <- standard_error(model, ...)
 
       stats::pt((range[1] - params$mu) / se$SE, df = dof, lower.tail = TRUE) +
         stats::pt((range[2] - params$mu) / se$SE, df = dof, lower.tail = FALSE)

diff --git a/R/p_direction.R b/R/p_direction.R
@@ -16,7 +16,9 @@ bayestestR::p_direction
 #' @param x A statistical model.
 #' @inheritParams bayestestR::p_direction
 #' @inheritParams model_parameters.default
-#' @param ... Arguments passed to other methods, e.g. `ci()`.
+#' @param ... Arguments passed to other methods, e.g. `ci()`. Arguments like
+#' `vcov` or `vcov_args` can be used to compute confidence intervals using a
+#' specific variance-covariance matrix for the standard errors.
 #'
 #' @seealso See also [`equivalence_test()`], [`p_function()`] and
 #' [`p_significance()`] for functions related to checking effect existence and
@@ -70,11 +72,14 @@ bayestestR::p_direction
 #'
 #' @return A data frame.
 #'
-#' @examplesIf requireNamespace("bayestestR") && require("see", quietly = TRUE)
+#' @examplesIf requireNamespace("bayestestR") && require("see", quietly = TRUE) && requireNamespace("sandwich")
 #' data(qol_cancer)
 #' model <- lm(QoL ~ time + age + education, data = qol_cancer)
 #' p_direction(model)
 #'
+#' # based on heteroscedasticity-robust standard errors
+#' p_direction(model, vcov = "HC3")
+#'
 #' result <- p_direction(model)
 #' plot(result)
 #' @export

diff --git a/R/p_significance.R b/R/p_significance.R
@@ -28,7 +28,9 @@ bayestestR::p_significance
 #' @inheritParams bayestestR::p_significance
 #' @inheritParams model_parameters.default
 #' @param verbose Toggle warnings and messages.
-#' @param ... Arguments passed to other methods, e.g. `ci()`.
+#' @param ... Arguments passed to other methods, e.g. `ci()`. Arguments like
+#' `vcov` or `vcov_args` can be used to compute confidence intervals using a
+#' specific variance-covariance matrix for the standard errors.
 #'
 #' @seealso For more details, see [`bayestestR::p_significance()`]. See also
 #' [`equivalence_test()`], [`p_function()`] and [`bayestestR::p_direction()`]
@@ -126,14 +128,16 @@ bayestestR::p_significance
 #' intervals and the values for practical significance. Higher values indicate
 #' more practical significance (upper bound is one).
 #'
-#' @examplesIf requireNamespace("bayestestR") && packageVersion("bayestestR") > "0.14.0"
+#' @examplesIf requireNamespace("bayestestR") && packageVersion("bayestestR") > "0.14.0" && requireNamespace("sandwich")
 #' data(qol_cancer)
 #' model <- lm(QoL ~ time + age + education, data = qol_cancer)
 #'
 #' p_significance(model)
 #' p_significance(model, threshold = c(-0.5, 1.5))
 #'
-#' # plot method
+#' # based on heteroscedasticity-robust standard errors
+#' p_significance(model, vcov = "HC3")
+#'
 #' if (require("see", quietly = TRUE)) {
 #'   result <- p_significance(model)
 #'   plot(result)

diff --git a/man/ci.default.Rd b/man/ci.default.Rd
diff --git a/man/cluster_analysis.Rd b/man/cluster_analysis.Rd
diff --git a/man/equivalence_test.lm.Rd b/man/equivalence_test.lm.Rd
diff --git a/man/p_direction.lm.Rd b/man/p_direction.lm.Rd
diff --git a/man/p_significance.lm.Rd b/man/p_significance.lm.Rd
diff --git a/man/p_value_betwithin.Rd b/man/p_value_betwithin.Rd
diff --git a/man/p_value_ml1.Rd b/man/p_value_ml1.Rd
diff --git a/man/p_value_satterthwaite.Rd b/man/p_value_satterthwaite.Rd
diff --git a/tests/testthat/_snaps/equivalence_test.md b/tests/testthat/_snaps/equivalence_test.md
@@ -15,3 +15,20 @@
       cyl         | [-1.94,  0.32] | 0.351  |   Undecided | 0.644 
       hp          | [-0.05,  0.01] | > .999 |    Accepted | < .001
 
+# equivalence_test, robust
+
+    Code
+      print(x)
+    Output
+      # TOST-test for Practical Equivalence
+      
+        ROPE: [-0.60 0.60]
+      
+      Parameter   |         90% CI |   SGPV | Equivalence |      p
+      ------------------------------------------------------------
+      (Intercept) | [23.10, 50.28] | < .001 |    Rejected | > .999
+      gear        | [-1.63,  2.36] | 0.421  |   Undecided | 0.628 
+      wt          | [-4.59, -1.45] | 0.001  |    Rejected | 0.993 
+      cyl         | [-2.24,  0.62] | 0.361  |   Undecided | 0.649 
+      hp          | [-0.05,  0.01] | > .999 |    Accepted | < .001
+
diff --git a/tests/testthat/_snaps/p_significance.md b/tests/testthat/_snaps/p_significance.md
@@ -13,3 +13,18 @@
       cyl         | [-2.17,  0.55] | 61.88%
       hp          | [-0.05,  0.01] |  0.00%
 
+# p_significance, robust
+
+    Code
+      print(x)
+    Output
+      Practical Significance (threshold: -0.60, 0.60)
+      
+      Parameter   |         95% CI |     ps
+      -------------------------------------
+      (Intercept) | [20.32, 53.06] |   100%
+      gear        | [-2.04,  2.77] | 41.23%
+      wt          | [-4.91, -1.13] | 99.39%
+      cyl         | [-2.53,  0.91] | 59.51%
+      hp          | [-0.06,  0.01] |  0.00%
+
diff --git a/tests/testthat/test-equivalence_test.R b/tests/testthat/test-equivalence_test.R
@@ -10,6 +10,14 @@ test_that("equivalence_test", {
   expect_snapshot(print(x))
 })
 
+test_that("equivalence_test, robust", {
+  skip_if_not_installed("sandwich")
+  data(mtcars)
+  m <- lm(mpg ~ gear + wt + cyl + hp, data = mtcars)
+  x <- equivalence_test(m, vcov = "HC3")
+  expect_snapshot(print(x))
+})
+
 test_that("equivalence_test, unequal rope-range", {
   data(iris)
   m <- lm(Sepal.Length ~ Species, data = iris)

diff --git a/tests/testthat/test-p_function.R b/tests/testthat/test-p_function.R
@@ -1,5 +1,3 @@
-skip_if_not_installed("sandwich")
-
 data(iris)
 model <- lm(Sepal.Length ~ Species, data = iris)
 
@@ -43,6 +41,7 @@ test_that("p_function ci-levels", {
     tolerance = 1e-4
   )
 
+  skip_if_not_installed("sandwich")
   out <- p_function(model, vcov = "HC3")
   expect_equal(
     out$CI_low,