From 8866833f2c957ef6185d14d5ff543dda932412ec Mon Sep 17 00:00:00 2001 From: Daniel Date: Wed, 4 Sep 2024 15:09:24 +0200 Subject: [PATCH] draft p_direction (#1010) * draft p_direction * typos --- DESCRIPTION | 2 +- NAMESPACE | 19 +++ NEWS.md | 4 + R/equivalence_test.R | 4 +- R/p_direction.R | 198 ++++++++++++++++++++++++ R/p_significance.R | 51 ++++-- inst/WORDLIST | 2 + man/equivalence_test.lm.Rd | 4 +- man/p_direction.lm.Rd | 196 +++++++++++++++++++++++ man/p_significance.lm.Rd | 25 +-- man/reexports.Rd | 5 +- pkgdown/_pkgdown.yml | 1 + tests/testthat/_snaps/p_direction.md | 15 ++ tests/testthat/_snaps/p_significance.md | 14 +- tests/testthat/test-p_direction.R | 41 +++++ tests/testthat/test-p_significance.R | 21 +++ 16 files changed, 567 insertions(+), 35 deletions(-) create mode 100644 R/p_direction.R create mode 100644 man/p_direction.lm.Rd create mode 100644 tests/testthat/_snaps/p_direction.md create mode 100644 tests/testthat/test-p_direction.R diff --git a/DESCRIPTION b/DESCRIPTION index f42ee69cb..6cd3d46ec 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: parameters Title: Processing of Model Parameters -Version: 0.22.2.2 +Version: 0.22.2.3 Authors@R: c(person(given = "Daniel", family = "Lüdecke", diff --git a/NAMESPACE b/NAMESPACE index 411498706..f7bbef008 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -350,6 +350,22 @@ S3method(model_parameters,zeroinfl) S3method(model_parameters,zoo) S3method(p_calibrate,default) S3method(p_calibrate,numeric) +S3method(p_direction,coxph) +S3method(p_direction,feis) +S3method(p_direction,felm) +S3method(p_direction,gee) +S3method(p_direction,glm) +S3method(p_direction,glmmTMB) +S3method(p_direction,gls) +S3method(p_direction,hurdle) +S3method(p_direction,lm) +S3method(p_direction,lme) +S3method(p_direction,merMod) +S3method(p_direction,mixed) +S3method(p_direction,rma) +S3method(p_direction,svyglm) +S3method(p_direction,wbm) +S3method(p_direction,zeroinfl) S3method(p_significance,coxph) S3method(p_significance,feis) S3method(p_significance,felm) @@ -563,6 +579,7 @@ S3method(print,n_clusters_hclust) S3method(print,n_clusters_silhouette) S3method(print,n_factors) S3method(print,p_calibrate) +S3method(print,p_direction_lm) S3method(print,p_significance_lm) S3method(print,parameters_brms_meta) S3method(print,parameters_clusters) @@ -930,6 +947,7 @@ export(n_components) export(n_factors) export(n_parameters) export(p_calibrate) +export(p_direction) export(p_function) export(p_significance) export(p_value) @@ -969,6 +987,7 @@ export(supported_models) export(visualisation_recipe) importFrom(bayestestR,ci) importFrom(bayestestR,equivalence_test) +importFrom(bayestestR,p_direction) importFrom(bayestestR,p_significance) importFrom(datawizard,demean) importFrom(datawizard,describe_distribution) diff --git a/NEWS.md b/NEWS.md index 62095656d..e6a90255e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,10 @@ * Used more accurate analytic approach to calculate normal distributions for the SGPV in `equivalence_test()` and used in `p_significance()`. +* Added `p_direction()` methods for frequentist models. This is a convenient + way to test the direction of the effect, which formerly was already (and still + is) possible with `pd = TRUE` in `model_parameters()`. + # parameters 0.22.2 ## New supported models diff --git a/R/equivalence_test.R b/R/equivalence_test.R index 83020a587..86371a0dc 100644 --- a/R/equivalence_test.R +++ b/R/equivalence_test.R @@ -105,7 +105,9 @@ bayestestR::equivalence_test #' represents _the proportion of data-supported hypotheses that are also null #' hypotheses_ _(Blume et al. 2018, Lakens and Delacre 2020)_. It represents the #' proportion of the _full_ confidence interval range (assuming a normally -#' distributed, equal-tailed interval) that is inside the ROPE. +#' distributed, equal-tailed interval) that is inside the ROPE. The SGPV ranges +#' from zero to one. Higher values indicate that the effect is more likely to be +#' practically equivalent ("not of interest"). #' #' Note that the assumed interval, which is used to calculate the SGPV, is an #' estimation of the _full interval_ based on the chosen confidence level. For diff --git a/R/p_direction.R b/R/p_direction.R new file mode 100644 index 000000000..6c509cb3e --- /dev/null +++ b/R/p_direction.R @@ -0,0 +1,198 @@ +#' @importFrom bayestestR p_direction +#' @export +bayestestR::p_direction + + +#' @title Probability of Direction (pd) +#' @name p_direction.lm +#' +#' @description Compute the **Probability of Direction** (*pd*, also known as +#' the Maximum Probability of Effect - *MPE*). This can be interpreted as the +#' probability that a parameter (described by its full confidence, or +#' "compatibility" interval) is strictly positive or negative (whichever is the +#' most probable). Although differently expressed, this index is fairly similar +#' (i.e., is strongly correlated) to the frequentist *p-value* (see 'Details'). +#' +#' @param x A statistical model. +#' @inheritParams bayestestR::p_direction +#' @inheritParams model_parameters.default +#' @param ... Arguments passed to other methods, e.g. `ci()`. +#' +#' @seealso See also [`equivalence_test()`], [`p_function()`] and +#' [`p_significance()`] for functions related to checking effect existence and +#' significance. +#' +#' @inheritSection bayestestR::p_direction What is the *pd*? +#' +#' @inheritSection bayestestR::p_direction Relationship with the p-value +#' +#' @inheritSection bayestestR::p_direction Possible Range of Values +#' +#' @inheritSection model_parameters Statistical inference - how to quantify evidence +#' +#' @references +#' +#' - Amrhein, V., Korner-Nievergelt, F., and Roth, T. (2017). The earth is +#' flat (p > 0.05): Significance thresholds and the crisis of unreplicable +#' research. PeerJ, 5, e3544. \doi{10.7717/peerj.3544} +#' +#' - Greenland S, Rafi Z, Matthews R, Higgs M. To Aid Scientific Inference, +#' Emphasize Unconditional Compatibility Descriptions of Statistics. (2022) +#' https://arxiv.org/abs/1909.08583v7 (Accessed November 10, 2022) +#' +#' - Lakens, D. (2024). Improving Your Statistical Inferences (Version v1.5.1). +#' Retrieved from https://lakens.github.io/statistical_inferences/. +#' \doi{10.5281/ZENODO.6409077} +#' +#' - Lakens, D., Scheel, A. M., and Isager, P. M. (2018). Equivalence Testing +#' for Psychological Research: A Tutorial. Advances in Methods and Practices +#' in Psychological Science, 1(2), 259–269. \doi{10.1177/2515245918770963} +#' +#' - Makowski, D., Ben-Shachar, M. S., Chen, S. H. A., and Lüdecke, D. (2019). +#' Indices of Effect Existence and Significance in the Bayesian Framework. +#' Frontiers in Psychology, 10, 2767. \doi{10.3389/fpsyg.2019.02767} +#' +#' - Rafi Z, Greenland S. Semantic and cognitive tools to aid statistical +#' science: replace confidence and significance by compatibility and surprise. +#' BMC Medical Research Methodology (2020) 20:244. +#' +#' - Schweder T. Confidence is epistemic probability for empirical science. +#' Journal of Statistical Planning and Inference (2018) 195:116–125. +#' \doi{10.1016/j.jspi.2017.09.016} +#' +#' - Schweder T, Hjort NL. Frequentist analogues of priors and posteriors. +#' In Stigum, B. (ed.), Econometrics and the Philosophy of Economics: Theory +#' Data Confrontation in Economics, pp. 285-217. Princeton University Press, +#' Princeton, NJ, 2003 +#' +#' - Vos P, Holbert D. Frequentist statistical inference without repeated sampling. +#' Synthese 200, 89 (2022). \doi{10.1007/s11229-022-03560-x} +#' +#' @return A data frame. +#' +#' @examplesIf requireNamespace("bayestestR") && require("see", quietly = TRUE) +#' data(qol_cancer) +#' model <- lm(QoL ~ time + age + education, data = qol_cancer) +#' p_direction(model) +#' +#' result <- p_direction(model) +#' plot(result) +#' @export +p_direction.lm <- function(x, + ci = 0.95, + method = "direct", + null = 0, + ...) { + # first, we need CIs + out <- ci(x, ci = ci, ...) + # we now iterate all confidence intervals and create an approximate normal + # 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, ci_range) + })) + colnames(posterior) <- out$Parameter + + out$pd <- as.numeric(bayestestR::p_direction( + posterior, + method = method, + null = null, + ... + )) + + # check we don't have duplicated columns in "posterior" we need this for + # plotting + if (anyDuplicated(colnames(posterior)) > 0 && !is.null(out$Component)) { + comps <- .rename_values(out$Component, "zero_inflated", "zi") + comps <- .rename_values(comps, "conditional", "cond") + colnames(posterior) <- paste0(out$Parameter, "_", comps) + out$Parameter <- paste0(out$Parameter, "_", comps) + } + + # we need these for plotting + if (!"Effects" %in% colnames(out)) { + out$Effects <- "fixed" + } + if (!"Component" %in% colnames(out)) { + out$Component <- "conditional" + } + + # reorder + out <- out[intersect(c("Parameter", "CI", "CI_low", "CI_high", "pd", "Effects", "Component"), colnames(out))] + + attr(out, "data") <- posterior + attr(out, "null") <- null + class(out) <- c("p_direction_lm", "p_direction", "see_p_direction", "data.frame") + out +} + + +# methods --------------------------------------------------------------------- + +#' @export +print.p_direction_lm <- function(x, digits = 2, p_digits = 3, ...) { + null <- attributes(x)$null + caption <- sprintf( + "Probability of Direction (null: %s)", + insight::format_value(null, digits = digits, protect_integer = TRUE) + ) + + # deal with Effects and Component columns + if ("Effects" %in% colnames(x) && insight::n_unique(x$Effects) == 1) { + x$Effects <- NULL + } + if ("Component" %in% colnames(x) && insight::n_unique(x$Component) == 1) { + x$Component <- NULL + } + + x <- insight::format_table(x, digits = digits, p_digits = p_digits) + cat(insight::export_table(x, title = caption, ...)) +} + + +# other classes -------------------------------------------------------------- + +#' @export +p_direction.glm <- p_direction.lm + +#' @export +p_direction.coxph <- p_direction.lm + +#' @export +p_direction.svyglm <- p_direction.lm + +#' @export +p_direction.glmmTMB <- p_direction.lm + +#' @export +p_direction.merMod <- p_direction.lm + +#' @export +p_direction.wbm <- p_direction.lm + +#' @export +p_direction.lme <- p_direction.lm + +#' @export +p_direction.gee <- p_direction.lm + +#' @export +p_direction.gls <- p_direction.lm + +#' @export +p_direction.feis <- p_direction.lm + +#' @export +p_direction.felm <- p_direction.lm + +#' @export +p_direction.mixed <- p_direction.lm + +#' @export +p_direction.hurdle <- p_direction.lm + +#' @export +p_direction.zeroinfl <- p_direction.lm + +#' @export +p_direction.rma <- p_direction.lm diff --git a/R/p_significance.R b/R/p_significance.R index f332a8f7b..3946abc1b 100644 --- a/R/p_significance.R +++ b/R/p_significance.R @@ -11,15 +11,18 @@ bayestestR::p_significance #' the probability that an effect is above a given threshold corresponding to a #' negligible effect in the median's direction, considering a parameter's _full_ #' confidence interval. In other words, it returns the probability of a clear -#' direction of an effect, which is larger than the smalles effect size of -#' interest (e.g., a minimal important difference). In comparison the the -#' [`equivalence_test()`] function, where the *SGPV* (second generation p-value) -#' describes the proportion of the _full_ confidence interval that is _inside_ -#' the ROPE, the value returned by `p_significance()` describes the _larger_ -#' proportion of the _full_ confidence interval that is _outside_ the ROPE. This -#' makes `p_significance()` comparable to [`bayestestR::p_direction()`], -#' however, while `p_direction()` compares to a point-null by default, -#' `p_significance()` compares to a range-null. +#' direction of an effect, which is larger than the smallest effect size of +#' interest (e.g., a minimal important difference). Its theoretical range is +#' from zero to one, but the *ps* is typically larger than 0.5 (to indicate +#' practical significance). +#' +#' In comparison the the [`equivalence_test()`] function, where the *SGPV* +#' (second generation p-value) describes the proportion of the _full_ confidence +#' interval that is _inside_ the ROPE, the value returned by `p_significance()` +#' describes the _larger_ proportion of the _full_ confidence interval that is +#' _outside_ the ROPE. This makes `p_significance()` comparable to +#' [`bayestestR::p_direction()`], however, while `p_direction()` compares to a +#' point-null by default, `p_significance()` compares to a range-null. #' #' @param x A statistical model. #' @inheritParams bayestestR::p_significance @@ -118,7 +121,9 @@ bayestestR::p_significance #' - Vos P, Holbert D. Frequentist statistical inference without repeated sampling. #' Synthese 200, 89 (2022). \doi{10.1007/s11229-022-03560-x} #' -#' @return A data frame. +#' @return A data frame with columns for the parameter names, the confidence +#' 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" #' data(qol_cancer) @@ -144,6 +149,23 @@ p_significance.lm <- function(x, threshold = "default", ci = 0.95, verbose = TRU })) colnames(posterior) <- out$Parameter + # deal with Effects and Component columns + if ("Effects" %in% colnames(out) && insight::n_unique(out$Effects) == 1) { + out$Effects <- NULL + } + if ("Component" %in% colnames(out) && insight::n_unique(out$Component) == 1) { + out$Component <- NULL + } + + # check we don't have duplicated columns in "posterior" we need this for + # plotting + if (anyDuplicated(colnames(posterior)) > 0 && !is.null(out$Component)) { + comps <- .rename_values(out$Component, "zero_inflated", "zi") + comps <- .rename_values(comps, "conditional", "cond") + colnames(posterior) <- paste0(out$Parameter, "_", comps) + out$Parameter <- paste0(out$Parameter, "_", comps) + } + # calculate the ROPE range if (all(threshold == "default")) { threshold <- bayestestR::rope_range(x, verbose = verbose) @@ -160,6 +182,9 @@ p_significance.lm <- function(x, threshold = "default", ci = 0.95, verbose = TRU threshold <- 0.1 } + # reorder + out <- out[intersect(c("Parameter", "CI", "CI_low", "CI_high", "ps", "Effects", "Component"), colnames(out))] + attr(out, "data") <- posterior attr(out, "threshold") <- threshold class(out) <- c("p_significance_lm", "p_significance", "see_p_significance", "data.frame") @@ -170,7 +195,7 @@ p_significance.lm <- function(x, threshold = "default", ci = 0.95, verbose = TRU # methods --------------------------------------------------------------------- #' @export -print.p_significance_lm <- function(x, digits = 2, p_digits = 3, ...) { +print.p_significance_lm <- function(x, digits = 2, ...) { threshold <- attributes(x)$threshold # make sure it's numeric if (!is.numeric(threshold)) { @@ -184,8 +209,8 @@ print.p_significance_lm <- function(x, digits = 2, p_digits = 3, ...) { "Practical Significance (threshold: %s)", toString(insight::format_value(threshold, digits = 2)) ) - x$ps <- insight::format_p(x$ps, name = "ps", digits = p_digits) - x <- insight::format_table(x, digits = digits, p_digits = p_digits) + x$ps <- insight::format_pd(x$ps, name = NULL) + x <- insight::format_table(x, digits = digits) cat(insight::export_table(x, title = caption, ...)) } diff --git a/inst/WORDLIST b/inst/WORDLIST index 6d8a576d3..6ab528def 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -14,6 +14,7 @@ Bayarri BayesFM BayesFactor Bentler +Bergh Biometrics Biometrika Blume @@ -95,6 +96,7 @@ Liu MADs MCMCglmm MLM +MPE MSA Maechler Malo diff --git a/man/equivalence_test.lm.Rd b/man/equivalence_test.lm.Rd index 4989d236d..70d440c8f 100644 --- a/man/equivalence_test.lm.Rd +++ b/man/equivalence_test.lm.Rd @@ -150,7 +150,9 @@ Second generation p-values (SGPV) were proposed as a statistic that represents \emph{the proportion of data-supported hypotheses that are also null hypotheses} \emph{(Blume et al. 2018, Lakens and Delacre 2020)}. It represents the proportion of the \emph{full} confidence interval range (assuming a normally -distributed, equal-tailed interval) that is inside the ROPE. +distributed, equal-tailed interval) that is inside the ROPE. The SGPV ranges +from zero to one. Higher values indicate that the effect is more likely to be +practically equivalent ("not of interest"). Note that the assumed interval, which is used to calculate the SGPV, is an estimation of the \emph{full interval} based on the chosen confidence level. For diff --git a/man/p_direction.lm.Rd b/man/p_direction.lm.Rd new file mode 100644 index 000000000..e6f0d36ea --- /dev/null +++ b/man/p_direction.lm.Rd @@ -0,0 +1,196 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/p_direction.R +\name{p_direction.lm} +\alias{p_direction.lm} +\title{Probability of Direction (pd)} +\usage{ +\method{p_direction}{lm}(x, ci = 0.95, method = "direct", null = 0, ...) +} +\arguments{ +\item{x}{A statistical model.} + +\item{ci}{Confidence Interval (CI) level. Default to \code{0.95} (\verb{95\%}).} + +\item{method}{Can be \code{"direct"} or one of methods of \code{\link[bayestestR:estimate_density]{estimate_density()}}, +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()}.} +} +\value{ +A data frame. +} +\description{ +Compute the \strong{Probability of Direction} (\emph{pd}, also known as +the Maximum Probability of Effect - \emph{MPE}). This can be interpreted as the +probability that a parameter (described by its full confidence, or +"compatibility" interval) is strictly positive or negative (whichever is the +most probable). Although differently expressed, this index is fairly similar +(i.e., is strongly correlated) to the frequentist \emph{p-value} (see 'Details'). +} +\section{What is the \emph{pd}?}{ + + + +The Probability of Direction (pd) is an index of effect existence, representing +the certainty with which an effect goes in a particular direction (i.e., is +positive or negative / has a sign), typically ranging from 0.5 to 1 (but see +next section for cases where it can range between 0 and 1). Beyond +its simplicity of interpretation, understanding and computation, this index +also presents other interesting properties: +\itemize{ +\item Like other posterior-based indices, \emph{pd} is solely based on the posterior +distributions and does not require any additional information from the data +or the model (e.g., such as priors, as in the case of Bayes factors). +\item It is robust to the scale of both the response variable and the predictors. +\item It is strongly correlated with the frequentist p-value, and can thus +be used to draw parallels and give some reference to readers non-familiar +with Bayesian statistics (Makowski et al., 2019). +} + +} + +\section{Relationship with the p-value}{ + + + +In most cases, it seems that the \emph{pd} has a direct correspondence with the +frequentist one-sided \emph{p}-value through the formula (for two-sided \emph{p}): +\deqn{p = 2 \times (1 - p_d)}{p = 2 * (1 - pd)} +Thus, a two-sided p-value of respectively \code{.1}, \code{.05}, \code{.01} and \code{.001} would +correspond approximately to a \emph{pd} of \verb{95\%}, \verb{97.5\%}, \verb{99.5\%} and \verb{99.95\%}. +See \code{\link[bayestestR:pd_to_p]{pd_to_p()}} for details. + +} + +\section{Possible Range of Values}{ + + + +The largest value \emph{pd} can take is 1 - the posterior is strictly directional. +However, the smallest value \emph{pd} can take depends on the parameter space +represented by the posterior. + +\strong{For a continuous parameter space}, exact values of 0 (or any point null +value) are not possible, and so 100\% of the posterior has \emph{some} sign, some +positive, some negative. Therefore, the smallest the \emph{pd} can be is 0.5 - +with an equal posterior mass of positive and negative values. Values close to +0.5 \emph{cannot} be used to support the null hypothesis (that the parameter does +\emph{not} have a direction) is a similar why to how large p-values cannot be used +to support the null hypothesis (see \code{\link[bayestestR:pd_to_p]{pd_to_p()}}; Makowski et al., 2019). + +\strong{For a discrete parameter space or a parameter space that is a mixture +between discrete and continuous spaces}, exact values of 0 (or any point +null value) \emph{are} possible! Therefore, the smallest the \emph{pd} can be is 0 - +with 100\% of the posterior mass on 0. Thus values close to 0 can be used to +support the null hypothesis (see van den Bergh et al., 2021). + +Examples of posteriors representing discrete parameter space: +\itemize{ +\item When a parameter can only take discrete values. +\item When a mixture prior/posterior is used (such as the spike-and-slab prior; +see van den Bergh et al., 2021). +\item When conducting Bayesian model averaging (e.g., \code{\link[bayestestR:weighted_posteriors]{weighted_posteriors()}} or +\code{brms::posterior_average}). +} + +} + +\section{Statistical inference - how to quantify evidence}{ + +There is no standardized approach to drawing conclusions based on the +available data and statistical models. A frequently chosen but also much +criticized approach is to evaluate results based on their statistical +significance (\emph{Amrhein et al. 2017}). + +A more sophisticated way would be to test whether estimated effects exceed +the "smallest effect size of interest", to avoid even the smallest effects +being considered relevant simply because they are statistically significant, +but clinically or practically irrelevant (\emph{Lakens et al. 2018, Lakens 2024}). + +A rather unconventional approach, which is nevertheless advocated by various +authors, is to interpret results from classical regression models in terms of +probabilities, similar to the usual approach in Bayesian statistics +(\emph{Greenland et al. 2022; Rafi and Greenland 2020; Schweder 2018; Schweder and +Hjort 2003; Vos 2022}). + +The \strong{parameters} package provides several options or functions to aid +statistical inference. These are, for example: +\itemize{ +\item \code{\link[=equivalence_test.lm]{equivalence_test()}}, to compute the (conditional) +equivalence test for frequentist models +\item \code{\link[=p_significance.lm]{p_significance()}}, to compute the probability of +\emph{practical significance}, which can be conceptualized as a unidirectional +equivalence test +\item \code{\link[=p_function]{p_function()}}, or \emph{consonance function}, to compute p-values and +compatibility (confidence) intervals for statistical models +\item the \code{pd} argument (setting \code{pd = TRUE}) in \code{model_parameters()} includes +a column with the \emph{probability of direction}, i.e. the probability that a +parameter is strictly positive or negative. See \code{\link[bayestestR:p_direction]{bayestestR::p_direction()}} +for details. +\item the \code{s_value} argument (setting \code{s_value = TRUE}) in \code{model_parameters()} +replaces the p-values with their related \emph{S}-values (\emph{Rafi and Greenland 2020}) +\item finally, it is possible to generate distributions of model coefficients by +generating bootstrap-samples (setting \code{bootstrap = TRUE}) or simulating +draws from model coefficients using \code{\link[=simulate_model]{simulate_model()}}. These samples +can then be treated as "posterior samples" and used in many functions from +the \strong{bayestestR} package. +} + +Most of the above shown options or functions derive from methods originally +implemented for Bayesian models (\emph{Makowski et al. 2019}). However, assuming +that model assumptions are met (which means, the model fits well to the data, +the correct model is chosen that reflects the data generating process +(distributional model family) etc.), it seems appropriate to interpret +results from classical frequentist models in a "Bayesian way" (more details: +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} +data(qol_cancer) +model <- lm(QoL ~ time + age + education, data = qol_cancer) +p_direction(model) + +result <- p_direction(model) +plot(result) +\dontshow{\}) # examplesIf} +} +\references{ +\itemize{ +\item Amrhein, V., Korner-Nievergelt, F., and Roth, T. (2017). The earth is +flat (p > 0.05): Significance thresholds and the crisis of unreplicable +research. PeerJ, 5, e3544. \doi{10.7717/peerj.3544} +\item Greenland S, Rafi Z, Matthews R, Higgs M. To Aid Scientific Inference, +Emphasize Unconditional Compatibility Descriptions of Statistics. (2022) +https://arxiv.org/abs/1909.08583v7 (Accessed November 10, 2022) +\item Lakens, D. (2024). Improving Your Statistical Inferences (Version v1.5.1). +Retrieved from https://lakens.github.io/statistical_inferences/. +\doi{10.5281/ZENODO.6409077} +\item Lakens, D., Scheel, A. M., and Isager, P. M. (2018). Equivalence Testing +for Psychological Research: A Tutorial. Advances in Methods and Practices +in Psychological Science, 1(2), 259–269. \doi{10.1177/2515245918770963} +\item Makowski, D., Ben-Shachar, M. S., Chen, S. H. A., and Lüdecke, D. (2019). +Indices of Effect Existence and Significance in the Bayesian Framework. +Frontiers in Psychology, 10, 2767. \doi{10.3389/fpsyg.2019.02767} +\item Rafi Z, Greenland S. Semantic and cognitive tools to aid statistical +science: replace confidence and significance by compatibility and surprise. +BMC Medical Research Methodology (2020) 20:244. +\item Schweder T. Confidence is epistemic probability for empirical science. +Journal of Statistical Planning and Inference (2018) 195:116–125. +\doi{10.1016/j.jspi.2017.09.016} +\item Schweder T, Hjort NL. Frequentist analogues of priors and posteriors. +In Stigum, B. (ed.), Econometrics and the Philosophy of Economics: Theory +Data Confrontation in Economics, pp. 285-217. Princeton University Press, +Princeton, NJ, 2003 +\item Vos P, Holbert D. Frequentist statistical inference without repeated sampling. +Synthese 200, 89 (2022). \doi{10.1007/s11229-022-03560-x} +} +} +\seealso{ +See also \code{\link[=equivalence_test]{equivalence_test()}}, \code{\link[=p_function]{p_function()}} and +\code{\link[=p_significance]{p_significance()}} for functions related to checking effect existence and +significance. +} diff --git a/man/p_significance.lm.Rd b/man/p_significance.lm.Rd index 835dc79ea..fb0858f0b 100644 --- a/man/p_significance.lm.Rd +++ b/man/p_significance.lm.Rd @@ -28,7 +28,9 @@ asymmetric intervals. \item{...}{Arguments passed to other methods, e.g. \code{ci()}.} } \value{ -A data frame. +A data frame with columns for the parameter names, the confidence +intervals and the values for practical significance. Higher values indicate +more practical significance (upper bound is one). } \description{ Compute the probability of \strong{Practical Significance} (\emph{ps}), @@ -36,15 +38,18 @@ which can be conceptualized as a unidirectional equivalence test. It returns the probability that an effect is above a given threshold corresponding to a negligible effect in the median's direction, considering a parameter's \emph{full} confidence interval. In other words, it returns the probability of a clear -direction of an effect, which is larger than the smalles effect size of -interest (e.g., a minimal important difference). In comparison the the -\code{\link[=equivalence_test]{equivalence_test()}} function, where the \emph{SGPV} (second generation p-value) -describes the proportion of the \emph{full} confidence interval that is \emph{inside} -the ROPE, the value returned by \code{p_significance()} describes the \emph{larger} -proportion of the \emph{full} confidence interval that is \emph{outside} the ROPE. This -makes \code{p_significance()} comparable to \code{\link[bayestestR:p_direction]{bayestestR::p_direction()}}, -however, while \code{p_direction()} compares to a point-null by default, -\code{p_significance()} compares to a range-null. +direction of an effect, which is larger than the smallest effect size of +interest (e.g., a minimal important difference). Its theoretical range is +from zero to one, but the \emph{ps} is typically larger than 0.5 (to indicate +practical significance). + +In comparison the the \code{\link[=equivalence_test]{equivalence_test()}} function, where the \emph{SGPV} +(second generation p-value) describes the proportion of the \emph{full} confidence +interval that is \emph{inside} the ROPE, the value returned by \code{p_significance()} +describes the \emph{larger} proportion of the \emph{full} confidence interval that is +\emph{outside} the ROPE. This makes \code{p_significance()} comparable to +\code{\link[bayestestR:p_direction]{bayestestR::p_direction()}}, however, while \code{p_direction()} compares to a +point-null by default, \code{p_significance()} compares to a range-null. } \details{ \code{p_significance()} returns the proportion of the \emph{full} confidence diff --git a/man/reexports.Rd b/man/reexports.Rd index 8825104df..04c2b5645 100644 --- a/man/reexports.Rd +++ b/man/reexports.Rd @@ -1,12 +1,13 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/equivalence_test.R, R/methods_bayestestR.R, -% R/n_parameters.R, R/p_significance.R, R/reexports.R +% R/n_parameters.R, R/p_direction.R, R/p_significance.R, R/reexports.R \docType{import} \name{reexports} \alias{reexports} \alias{equivalence_test} \alias{ci} \alias{n_parameters} +\alias{p_direction} \alias{p_significance} \alias{standardize_names} \alias{supported_models} @@ -26,7 +27,7 @@ These objects are imported from other packages. Follow the links below to see their documentation. \describe{ - \item{bayestestR}{\code{\link[bayestestR]{ci}}, \code{\link[bayestestR]{equivalence_test}}, \code{\link[bayestestR]{p_significance}}} + \item{bayestestR}{\code{\link[bayestestR]{ci}}, \code{\link[bayestestR]{equivalence_test}}, \code{\link[bayestestR]{p_direction}}, \code{\link[bayestestR]{p_significance}}} \item{datawizard}{\code{\link[datawizard]{demean}}, \code{\link[datawizard]{describe_distribution}}, \code{\link[datawizard:skewness]{kurtosis}}, \code{\link[datawizard]{rescale_weights}}, \code{\link[datawizard]{skewness}}, \code{\link[datawizard]{visualisation_recipe}}} diff --git a/pkgdown/_pkgdown.yml b/pkgdown/_pkgdown.yml index 1dd7d3dce..01bc03c9d 100644 --- a/pkgdown/_pkgdown.yml +++ b/pkgdown/_pkgdown.yml @@ -60,6 +60,7 @@ reference: contents: - equivalence_test.lm - p_calibrate + - p_direction.lm - p_function - p_significance.lm diff --git a/tests/testthat/_snaps/p_direction.md b/tests/testthat/_snaps/p_direction.md new file mode 100644 index 000000000..14e9c1fa9 --- /dev/null +++ b/tests/testthat/_snaps/p_direction.md @@ -0,0 +1,15 @@ +# p_direction + + Code + print(x) + Output + Probability of Direction (null: 0) + + Parameter | 95% CI | pd + ------------------------------------- + (Intercept) | [24.44, 48.94] | 100% + gear | [-1.69, 2.41] | 63.55% + wt | [-4.77, -1.28] | 99.97% + cyl | [-2.17, 0.55] | 87.91% + hp | [-0.05, 0.01] | 90.61% + diff --git a/tests/testthat/_snaps/p_significance.md b/tests/testthat/_snaps/p_significance.md index e3c090876..0678b1b85 100644 --- a/tests/testthat/_snaps/p_significance.md +++ b/tests/testthat/_snaps/p_significance.md @@ -5,11 +5,11 @@ Output Practical Significance (threshold: -0.60, 0.60) - Parameter | 95% CI | ps - ----------------------------------------- - (Intercept) | [24.44, 48.94] | ps > .999 - gear | [-1.69, 2.41] | ps = 0.409 - wt | [-4.77, -1.28] | ps = 0.997 - cyl | [-2.17, 0.55] | ps = 0.619 - hp | [-0.05, 0.01] | ps < .001 + Parameter | 95% CI | ps + ------------------------------------- + (Intercept) | [24.44, 48.94] | 100% + gear | [-1.69, 2.41] | 40.93% + wt | [-4.77, -1.28] | 99.67% + cyl | [-2.17, 0.55] | 61.88% + hp | [-0.05, 0.01] | 0.00% diff --git a/tests/testthat/test-p_direction.R b/tests/testthat/test-p_direction.R new file mode 100644 index 000000000..09e7a4540 --- /dev/null +++ b/tests/testthat/test-p_direction.R @@ -0,0 +1,41 @@ +skip_on_cran() +skip_if_not_installed("bayestestR") + +test_that("p_direction", { + data(mtcars) + m <- lm(mpg ~ gear + wt + cyl + hp, data = mtcars) + set.seed(123) + x <- p_direction(m) + expect_identical(c(nrow(x), ncol(x)), c(5L, 7L)) + expect_named(x, c("Parameter", "CI", "CI_low", "CI_high", "pd", "Effects", "Component")) + expect_snapshot(print(x)) + + set.seed(123) + x <- p_direction(m, ci = 0.8) + expect_equal(x$pd, c(1, 0.6382, 0.9997, 0.884, 0.9107), tolerance = 1e-4) + + set.seed(123) + x <- p_direction(m, null = 0.2) + expect_equal(x$pd, c(1, 0.5617, 0.9999, 0.9276, 1), tolerance = 1e-4) +}) + +test_that("p_direction, glmmTMB", { + skip_if_not_installed("glmmTMB") + data(Salamanders, package = "glmmTMB") + m1 <- glmmTMB::glmmTMB(count ~ mined + cover + (1 | site), + zi = ~mined, + family = poisson, + data = Salamanders + ) + out <- p_direction(m1) + expect_identical(c(nrow(out), ncol(out)), c(5L, 7L)) + expect_named(out, c("Parameter", "CI", "CI_low", "CI_high", "pd", "Effects", "Component")) + expect_equal(out$pd, c(0.8245, 1, 0.9974, 1, 1), tolerance = 1e-4) + expect_identical( + out$Parameter, + c( + "(Intercept)_cond", "minedno_cond", "cover_cond", "(Intercept)_zi", + "minedno_zi" + ) + ) +}) diff --git a/tests/testthat/test-p_significance.R b/tests/testthat/test-p_significance.R index 1f9a02f07..c95756e01 100644 --- a/tests/testthat/test-p_significance.R +++ b/tests/testthat/test-p_significance.R @@ -18,3 +18,24 @@ test_that("p_significance", { x <- p_significance(m, threshold = 0.5) expect_equal(x$ps, c(1, 0.4478, 0.9977, 0.6737, 0), tolerance = 1e-4) }) + +test_that("p_significance, glmmTMB", { + skip_if_not_installed("glmmTMB") + data(Salamanders, package = "glmmTMB") + m1 <- glmmTMB::glmmTMB(count ~ mined + cover + (1 | site), + zi = ~mined, + family = poisson, + data = Salamanders + ) + out <- p_significance(m1) + expect_identical(c(nrow(out), ncol(out)), c(5L, 6L)) + expect_named(out, c("Parameter", "CI", "CI_low", "CI_high", "ps", "Component")) + expect_equal(out$ps, c(0.6451, 1, 0.9015, 1, 1), tolerance = 1e-4) + expect_identical( + out$Parameter, + c( + "(Intercept)_cond", "minedno_cond", "cover_cond", "(Intercept)_zi", + "minedno_zi" + ) + ) +})