From af6dfbabb3b003b3c80ebfda172f120f64cf949b Mon Sep 17 00:00:00 2001 From: "Mattan S. Ben-Shachar" <35330040+mattansb@users.noreply.github.com> Date: Mon, 29 Jan 2024 15:38:46 +0200 Subject: [PATCH] Better ncp ci optimization (#629) * use q* instead of p* * v bump * update docs --- DESCRIPTION | 2 +- NEWS.md | 1 + R/docs_extra.R | 9 ++------- R/utils_ci.R | 41 ++++++++++++++++----------------------- man/cohens_d.Rd | 2 +- man/effectsize-package.Rd | 1 + man/effectsize_CIs.Rd | 12 ++---------- 7 files changed, 25 insertions(+), 43 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 45c7a283d..6446df71f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: effectsize Title: Indices of Effect Size -Version: 0.8.6.5 +Version: 0.8.6.6 Authors@R: c(person(given = "Mattan S.", family = "Ben-Shachar", diff --git a/NEWS.md b/NEWS.md index 4fa674822..8fe69bac8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -11,6 +11,7 @@ ## Bug fixes +- Minor stability fix to ncp-based CI methods ( #628 ) - `nnt()` now properly accepts the `y` argument. # effectsize 0.8.7 diff --git a/R/docs_extra.R b/R/docs_extra.R index 39249d9ff..cb133e6bc 100644 --- a/R/docs_extra.R +++ b/R/docs_extra.R @@ -124,13 +124,8 @@ #' smaller than the tolerance of the optimizer, resulting in CIs of width 0. #' This can also result in the estimated CIs excluding the point estimate. #' -#' For example: -#' ```{r} -#' t_to_d(80, df_error = 4555555) -#' ``` -#' -#' In these cases, consider an alternative optimizer, or an alternative method -#' for computing CIs, such as the bootstrap. +#' In these cases, consider an alternative method for computing CIs, such as the +#' bootstrap. #' #' #' @references diff --git a/R/utils_ci.R b/R/utils_ci.R index 20950705b..9f36db721 100644 --- a/R/utils_ci.R +++ b/R/utils_ci.R @@ -15,10 +15,8 @@ ncp <- suppressWarnings(stats::optim( par = 1.1 * rep(lambda, 2), fn = function(x) { - p <- stats::pf(q = f, df, df_error, ncp = x) - - abs(max(p) - probs[2]) + - abs(min(p) - probs[1]) + q <- stats::qf(p = probs, df, df_error, ncp = x) + sum(abs(q - f)) }, control = list(abstol = 1e-09) )) @@ -37,10 +35,9 @@ #' @keywords internal .get_ncp_t <- function(t, df_error, conf.level = 0.95) { - # # Note: these aren't actually needed - all t related functions would fail earlier - # if (!is.finite(t) || !is.finite(df_error)) { - # return(c(NA, NA)) - # } + if (!is.finite(t) || !is.finite(df_error)) { + return(c(NA, NA)) + } alpha <- 1 - conf.level probs <- c(alpha / 2, 1 - alpha / 2) @@ -48,45 +45,41 @@ ncp <- suppressWarnings(stats::optim( par = 1.1 * rep(t, 2), fn = function(x) { - p <- stats::pt(q = t, df = df_error, ncp = x) - - abs(max(p) - probs[2]) + - abs(min(p) - probs[1]) + q <- stats::qt(p = probs, df = df_error, ncp = x) + sum(abs(q - t)) }, control = list(abstol = 1e-09) )) + t_ncp <- unname(sort(ncp$par)) return(t_ncp) } #' @keywords internals -.get_ncp_chi <- function(chi, df, conf.level = 0.95) { - # # Note: these aren't actually needed - all chisq related functions would fail earlier - # if (!is.finite(chi) || !is.finite(df)) { - # return(c(NA, NA)) - # } +.get_ncp_chi <- function(chisq, df, conf.level = 0.95) { + if (!is.finite(chisq) || !is.finite(df)) { + return(c(NA, NA)) + } alpha <- 1 - conf.level probs <- c(alpha / 2, 1 - alpha / 2) ncp <- suppressWarnings(stats::optim( - par = 1.1 * rep(chi, 2), + par = 1.1 * rep(chisq, 2), fn = function(x) { - p <- stats::pchisq(q = chi, df, ncp = x) - - abs(max(p) - probs[2]) + - abs(min(p) - probs[1]) + q <- stats::qchisq(p = probs, df, ncp = x) + sum(abs(q - chisq)) }, control = list(abstol = 1e-09) )) chi_ncp <- sort(ncp$par) - if (chi <= stats::qchisq(probs[1], df)) { + if (chisq <= stats::qchisq(probs[1], df)) { chi_ncp[2] <- 0 } - if (chi <= stats::qchisq(probs[2], df)) { + if (chisq <= stats::qchisq(probs[2], df)) { chi_ncp[1] <- 0 } diff --git a/man/cohens_d.Rd b/man/cohens_d.Rd index 60236654e..6a67d0b3b 100644 --- a/man/cohens_d.Rd +++ b/man/cohens_d.Rd @@ -38,7 +38,7 @@ glass_delta( y = NULL, data = NULL, mu = 0, - adjust = FALSE, + adjust = TRUE, ci = 0.95, alternative = "two.sided", verbose = TRUE, diff --git a/man/effectsize-package.Rd b/man/effectsize-package.Rd index 30c9a4414..6cd2eb230 100644 --- a/man/effectsize-package.Rd +++ b/man/effectsize-package.Rd @@ -47,6 +47,7 @@ Authors: \item Indrajeet Patil \email{patilindrajeet.science@gmail.com} (\href{https://orcid.org/0000-0003-1995-6531}{ORCID}) (@patilindrajeets) \item Brenton M. Wiernik \email{brenton@wiernik.org} (\href{https://orcid.org/0000-0001-9560-6336}{ORCID}) (@bmwiernik) \item Rémi Thériault \email{remi.theriault@mail.mcgill.ca} (\href{https://orcid.org/0000-0003-4315-6788}{ORCID}) (@rempsyc) + \item Philip Waggoner \email{philip.waggoner@gmail.com} (\href{https://orcid.org/0000-0002-7825-7573}{ORCID}) [contributor] } Other contributors: diff --git a/man/effectsize_CIs.Rd b/man/effectsize_CIs.Rd index 8f848c601..db30b81ca 100644 --- a/man/effectsize_CIs.Rd +++ b/man/effectsize_CIs.Rd @@ -155,16 +155,8 @@ For very large sample sizes or effect sizes, the width of the CI can be smaller than the tolerance of the optimizer, resulting in CIs of width 0. This can also result in the estimated CIs excluding the point estimate. -For example: - -\if{html}{\out{
}}\preformatted{t_to_d(80, df_error = 4555555) -#> d | 95\% CI -#> ------------------- -#> 0.07 | [0.08, 0.08] -}\if{html}{\out{
}} - -In these cases, consider an alternative optimizer, or an alternative method -for computing CIs, such as the bootstrap. +In these cases, consider an alternative method for computing CIs, such as the +bootstrap. } \references{