Better ncp ci optimization

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",

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

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
@@ -140,7 +135,7 @@
 #' \doi{10.1093/biomet/83.4.934}
 #'
 #' Rafi, Z., & Greenland, S. (2020).
 #' Semantic and cognitive tools to aid statistical science: Replace confidence and significance by compatibility and surprise.
 #' _BMC Medical Research Methodology, 20_(1), Article 244.
 #' \doi{10.1186/s12874-020-01105-9}
 #'
@@ -150,7 +145,7 @@
 #' \doi{10.1017/CBO9781139046671}
 #'
 #' Steiger, J. H. (2004).
 #' Beyond the _F_ test: Effect size confidence intervals and tests of close fit in the analysis of variance and contrast analysis.
 #' _Psychological Methods, 9_(2), 164--182.
 #' \doi{10.1037/1082-989x.9.2.164}
 #'
@@ -166,7 +161,7 @@
 
 #' `effectsize` API
 #'
 #' Read the [*Support functions for model extensions*](https://easystats.github.io/effectsize/articles/effectsize_API.html) vignette.
 #'
 #' @rdname effectsize_API
 #' @name effectsize_API
@@ -177,7 +172,7 @@
 #' `effectsize` options
 #'
 #' Currently, the following global options are supported:
 #' - `es.use_symbols` [logical]: Should proper symbols be printed (`TRUE`) instead of transliterated effect size names (`FALSE`; default).
 #'
 #' @rdname effectsize_options
 #' @name effectsize_options

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,56 +35,51 @@
 
 #' @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)
 
   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
   }
 
@@ -107,7 +100,7 @@
     ci > 1) {
     insight::format_error("ci must be a single numeric value between (0, 1)")
   }
   return(TRUE)
 }
 
 #' @keywords internal