From a0688510334412687702c69a64a49879210caf6b Mon Sep 17 00:00:00 2001 From: Daniel Date: Mon, 16 Sep 2024 12:35:01 +0200 Subject: [PATCH] Add `vcov` option to docs and examples (#1014) * Add `vcov` option to docs and examples * news desc * docs --- DESCRIPTION | 2 +- NEWS.md | 4 ++++ R/2_ci.R | 16 ++++++++++++++-- R/equivalence_test.R | 20 +++++++++++++------- R/p_direction.R | 9 +++++++-- R/p_function.R | 8 +++++--- R/p_significance.R | 10 +++++++--- man/ci.default.Rd | 16 ++++++++++++++-- man/cluster_analysis.Rd | 5 ++++- man/equivalence_test.lm.Rd | 10 +++++++++- man/p_direction.lm.Rd | 9 +++++++-- man/p_function.Rd | 8 +++++--- man/p_significance.lm.Rd | 10 +++++++--- man/p_value_betwithin.Rd | 5 ++++- man/p_value_ml1.Rd | 5 ++++- man/p_value_satterthwaite.Rd | 5 ++++- tests/testthat/_snaps/equivalence_test.md | 17 +++++++++++++++++ tests/testthat/_snaps/p_significance.md | 15 +++++++++++++++ tests/testthat/test-equivalence_test.R | 8 ++++++++ tests/testthat/test-p_function.R | 3 +-- tests/testthat/test-p_significance.R | 9 +++++++++ 21 files changed, 159 insertions(+), 35 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 5aec96ffa..a2bf90f45 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: parameters Title: Processing of Model Parameters -Version: 0.22.2.8 +Version: 0.22.2.9 Authors@R: c(person(given = "Daniel", family = "Lüdecke", diff --git a/NEWS.md b/NEWS.md index 39b9a541d..d8cda798b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -12,6 +12,10 @@ * `p_function()` gets a `vcov` and `vcov_args` argument to compute robust standard errors for the confidence curves. +* Functions `p_significance()` and `equivalence_test()` now pass arguments + `vcov` and `vcov_args` to `p_value()` and `ci()`, hence, tests can be based + on robust standard errors. + * Revision / enhancement of some documentation. # parameters 0.22.2 diff --git a/R/2_ci.R b/R/2_ci.R index 79e3e59da..c010dcf37 100644 --- 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 index 3d2b44ea5..71f07d11f 100644 --- 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 index 4030cbc85..c0d975663 100644 --- 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_function.R b/R/p_function.R index 54596f701..cfb5aeda3 100644 --- a/R/p_function.R +++ b/R/p_function.R @@ -119,9 +119,11 @@ #' whereas values values far from the observed point estimate (where _p_ #' approaches 0) are least compatible with the data and model assumptions #' (_Schweder and Hjort 2016, pp. 60-61; Amrhein and Greenland 2022_). In this -#' regards, _p_-values are are statements about confidence/compatibility, not -#' probability per se, but still the interpretation of _p_-values might be -#' guided using [`bayestestR::p_to_pd()`] +#' regards, _p_-values are statements about _confidence_ or _compatibility_ and +#' can be considered as _epistemic probability_ - "not necessarily of the +#' hypothesis being true, but of it _possibly_ being true" (_Schweder 2018_). +#' Hence, the interpretation of _p_-values might be guided using +#' [`bayestestR::p_to_pd()`] #' #' ## Compatibility intervals - is their interpretation conditional or not? #' diff --git a/R/p_significance.R b/R/p_significance.R index 2629d793f..d22e06e7d 100644 --- 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 index 2e86093f1..e1ca2e1eb 100644 --- a/man/ci.default.Rd +++ b/man/ci.default.Rd @@ -40,7 +40,10 @@ options (which vary depending on the model class): \code{"residual"}, \emph{Confidence intervals and approximation of degrees of freedom} in \code{\link[=model_parameters]{model_parameters()}} for further details.} -\item{...}{Additional arguments} +\item{...}{Additional arguments passed down to the underlying functions. +E.g., arguments like \code{vcov} or \code{vcov_args} can be used to compute confidence +intervals using a specific variance-covariance matrix for the standard +errors.} \item{component}{Model component for which parameters should be shown. See the documentation for your object's class in \code{\link[=model_parameters]{model_parameters()}} or @@ -224,7 +227,16 @@ which is converted into a p-value using \code{\link[bayestestR:pd_to_p]{bayestes } \examples{ -\dontshow{if (require("glmmTMB")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (require("glmmTMB") && requireNamespace("sandwich")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +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/man/cluster_analysis.Rd b/man/cluster_analysis.Rd index 7d2e9cd81..c7659aeeb 100644 --- a/man/cluster_analysis.Rd +++ b/man/cluster_analysis.Rd @@ -69,7 +69,10 @@ this argument.} \item{iterations}{The number of replications.} -\item{...}{Arguments passed to or from other methods.} +\item{...}{Arguments passed to or from other methods, e.g. \code{ci()}. Arguments +like \code{vcov} or \code{vcov_args} can be used to compute confidence intervals or +p-values using a specific variance-covariance matrix for the standard +errors..} } \value{ The group classification for each observation as vector. The diff --git a/man/equivalence_test.lm.Rd b/man/equivalence_test.lm.Rd index 8796960d9..9e8ec0b0c 100644 --- a/man/equivalence_test.lm.Rd +++ b/man/equivalence_test.lm.Rd @@ -48,7 +48,10 @@ equivalence. Can be \code{"bayes"}, \code{"classic"} or \code{"cet"}. See 'Detai \item{verbose}{Toggle warnings and messages.} -\item{...}{Arguments passed to or from other methods.} +\item{...}{Arguments passed to or from other methods, e.g. \code{ci()}. Arguments +like \code{vcov} or \code{vcov_args} can be used to compute confidence intervals or +p-values using a specific variance-covariance matrix for the standard +errors..} \item{effects}{Should parameters for fixed effects (\code{"fixed"}), random effects (\code{"random"}), or both (\code{"all"}) be returned? Only applies @@ -249,12 +252,16 @@ documentation in \code{\link[=p_function]{p_function()}}). } \examples{ +\dontshow{if (requireNamespace("sandwich")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} 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") @@ -263,6 +270,7 @@ if (require("see", quietly = TRUE)) { result <- equivalence_test(model) plot(result) } +\dontshow{\}) # examplesIf} } \references{ \itemize{ diff --git a/man/p_direction.lm.Rd b/man/p_direction.lm.Rd index c057861be..2337eefe2 100644 --- a/man/p_direction.lm.Rd +++ b/man/p_direction.lm.Rd @@ -17,7 +17,9 @@ such as \code{"kernel"}, \code{"logspline"} or \code{"KernSmooth"}. See details. \item{null}{The value considered as a "null" effect. Traditionally 0, but could also be 1 in the case of ratios of change (OR, IRR, ...).} -\item{...}{Arguments passed to other methods, e.g. \code{ci()}.} +\item{...}{Arguments passed to other methods, e.g. \code{ci()}. Arguments like +\code{vcov} or \code{vcov_args} can be used to compute confidence intervals using a +specific variance-covariance matrix for the standard errors.} } \value{ A data frame. @@ -150,11 +152,14 @@ documentation in \code{\link[=p_function]{p_function()}}). } \examples{ -\dontshow{if (requireNamespace("bayestestR") && require("see", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (requireNamespace("bayestestR") && require("see", quietly = TRUE) && requireNamespace("sandwich")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} 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) \dontshow{\}) # examplesIf} diff --git a/man/p_function.Rd b/man/p_function.Rd index b7b2be9c0..77a587a17 100644 --- a/man/p_function.Rd +++ b/man/p_function.Rd @@ -248,9 +248,11 @@ estimated to be \emph{most compatible} with the data and model assumptions, whereas values values far from the observed point estimate (where \emph{p} approaches 0) are least compatible with the data and model assumptions (\emph{Schweder and Hjort 2016, pp. 60-61; Amrhein and Greenland 2022}). In this -regards, \emph{p}-values are are statements about confidence/compatibility, not -probability per se, but still the interpretation of \emph{p}-values might be -guided using \code{\link[bayestestR:pd_to_p]{bayestestR::p_to_pd()}} +regards, \emph{p}-values are statements about \emph{confidence} or \emph{compatibility} and +can be considered as \emph{epistemic probability} - "not necessarily of the +hypothesis being true, but of it \emph{possibly} being true" (\emph{Schweder 2018}). +Hence, the interpretation of \emph{p}-values might be guided using +\code{\link[bayestestR:pd_to_p]{bayestestR::p_to_pd()}} } \subsection{Compatibility intervals - is their interpretation conditional or not?}{ diff --git a/man/p_significance.lm.Rd b/man/p_significance.lm.Rd index 4c598e487..37a0adb36 100644 --- a/man/p_significance.lm.Rd +++ b/man/p_significance.lm.Rd @@ -25,7 +25,9 @@ asymmetric intervals. \item{verbose}{Toggle warnings and messages.} -\item{...}{Arguments passed to other methods, e.g. \code{ci()}.} +\item{...}{Arguments passed to other methods, e.g. \code{ci()}. Arguments like +\code{vcov} or \code{vcov_args} can be used to compute confidence intervals using a +specific variance-covariance matrix for the standard errors.} } \value{ A data frame with columns for the parameter names, the confidence @@ -151,14 +153,16 @@ documentation in \code{\link[=p_function]{p_function()}}). } \examples{ -\dontshow{if (requireNamespace("bayestestR") && packageVersion("bayestestR") > "0.14.0") (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (requireNamespace("bayestestR") && packageVersion("bayestestR") > "0.14.0" && requireNamespace("sandwich")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} 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/p_value_betwithin.Rd b/man/p_value_betwithin.Rd index 17d814646..9f2c7069d 100644 --- a/man/p_value_betwithin.Rd +++ b/man/p_value_betwithin.Rd @@ -18,7 +18,10 @@ p_value_betwithin(model, dof = NULL, ...) \item{ci}{Confidence Interval (CI) level. Default to \code{0.95} (\verb{95\%}).} -\item{...}{Additional arguments} +\item{...}{Additional arguments passed down to the underlying functions. +E.g., arguments like \code{vcov} or \code{vcov_args} can be used to compute confidence +intervals using a specific variance-covariance matrix for the standard +errors.} \item{dof}{Degrees of Freedom.} } diff --git a/man/p_value_ml1.Rd b/man/p_value_ml1.Rd index 0fb9229f6..c272fe8d9 100644 --- a/man/p_value_ml1.Rd +++ b/man/p_value_ml1.Rd @@ -17,7 +17,10 @@ p_value_ml1(model, dof = NULL, ...) \item{ci}{Confidence Interval (CI) level. Default to \code{0.95} (\verb{95\%}).} -\item{...}{Additional arguments} +\item{...}{Additional arguments passed down to the underlying functions. +E.g., arguments like \code{vcov} or \code{vcov_args} can be used to compute confidence +intervals using a specific variance-covariance matrix for the standard +errors.} \item{dof}{Degrees of Freedom.} } diff --git a/man/p_value_satterthwaite.Rd b/man/p_value_satterthwaite.Rd index 5ac41192c..ebe3586b3 100644 --- a/man/p_value_satterthwaite.Rd +++ b/man/p_value_satterthwaite.Rd @@ -21,7 +21,10 @@ se_satterthwaite(model) \item{ci}{Confidence Interval (CI) level. Default to \code{0.95} (\verb{95\%}).} -\item{...}{Additional arguments} +\item{...}{Additional arguments passed down to the underlying functions. +E.g., arguments like \code{vcov} or \code{vcov_args} can be used to compute confidence +intervals using a specific variance-covariance matrix for the standard +errors.} \item{dof}{Degrees of Freedom.} } diff --git a/tests/testthat/_snaps/equivalence_test.md b/tests/testthat/_snaps/equivalence_test.md index 8bec29040..325858efa 100644 --- 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 index b80688ac1..81b4cb7ce 100644 --- 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 index 3c95f7625..e358b690d 100644 --- 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 index 9d937aef4..c1090a2cf 100644 --- 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, diff --git a/tests/testthat/test-p_significance.R b/tests/testthat/test-p_significance.R index 7ec769851..092abf8bf 100644 --- a/tests/testthat/test-p_significance.R +++ b/tests/testthat/test-p_significance.R @@ -39,3 +39,12 @@ test_that("p_significance, glmmTMB", { ) ) }) + +test_that("p_significance, robust", { + skip_if_not_installed("sandwich") + data(mtcars) + m <- lm(mpg ~ gear + wt + cyl + hp, data = mtcars) + set.seed(123) + x <- p_significance(m, vcov = "HC3") + expect_snapshot(print(x)) +})