Allow to return the leave-one-values of jackknife()

NEWS.md

@@ -3,6 +3,12 @@
 * Add `math_lcm()` and `math_gcd()` to compute the lowest common multiple and the greatest common divisor.
 * Add `interval_hdr()` and `interval_credible()` to compute the credible intervals.
 
+## Bugfixes & changes
+* `jackknife()` gained a new argument to apply a function on the leave-one-out values (`f`).
+
+## Internals
+* Add `with_seed()` to evaluate an expression with a temporarily seed.
+
 # arkhe 1.1.0
 ## New classes and methods
 * Add `needs()` to check for the availability of a package.

R/AllGenerics.R

@@ -494,17 +494,17 @@ setGeneric(
 #'  `do`) as argument.
 #' @param ... Extra arguments to be passed to `do`.
 #' @return
-#'  If `f` is not `NULL`, `bootstrap()` returns the result of `f` applied to
-#'  the `n` values of `do`.
-#'
-#'  If `f` is `NULL`, `bootstrap()` returns a named `numeric` `vector` with the
-#'  following elements:
+#'  If `f` is `NULL` (the default), `bootstrap()` returns a named `numeric`
+#'  vector with the following elements:
 #'  \describe{
 #'   \item{`original`}{The observed value of `do` applied to `object`.}
 #'   \item{`mean`}{The bootstrap estimate of mean of `do`.}
 #'   \item{`bias`}{The bootstrap estimate of bias of `do`.}
 #'   \item{`error`}{he bootstrap estimate of standard error of `do`.}
 #'  }
+#'
+#'  If `f` is a `function`, `bootstrap()` returns the result of `f` applied to
+#'  the `n` values of `do`.
 #' @example inst/examples/ex-resample.R
 #' @author N. Frerebeau
 #' @docType methods
@@ -523,14 +523,20 @@ setGeneric(
 #' @param do A [`function`] that takes `object` as an argument and returns a
 #'  single numeric value.
 #' @param ... Extra arguments to be passed to `do`.
+#' @param f A [`function`] that takes a single numeric vector (the leave-one-out
+#'  values of `do`) as argument.
 #' @return
-#'  Returns a named `numeric` `vector` with the following elements:
+#'  If `f` is `NULL` (the default), `jackknife()` returns a named `numeric`
+#'  vector with the following elements:
 #'  \describe{
 #'   \item{`original`}{The observed value of `do` applied to `object`.}
 #'   \item{`mean`}{The jackknife estimate of mean of `do`.}
 #'   \item{`bias`}{The jackknife estimate of bias of `do`.}
 #'   \item{`error`}{he jackknife estimate of standard error of `do`.}
 #'  }
+#'
+#'  If `f` is a `function`, `jackknife()` returns the result of `f` applied to
+#'  the leave-one-out values of `do`.
 #' @example inst/examples/ex-resample.R
 #' @author N. Frerebeau
 #' @docType methods

R/statistics.R

@@ -164,8 +164,9 @@ setMethod(
     spl <- sample(object, size = length(object) * n, replace = TRUE)
     replicates <- t(matrix(spl, nrow = n))
     values <- apply(X = replicates, MARGIN = 2, FUN = do, ...)
-    values <- if (is.function(f)) f(values) else summary_bootstrap(values, hat)
-    values
+
+    if (is.function(f)) return(f(values))
+    summary_bootstrap(values, hat)
   }
 )
 
@@ -187,7 +188,7 @@ summary_bootstrap <- function(x, hat) {
 setMethod(
   f = "jackknife",
   signature = c(object = "numeric"),
-  definition = function(object, do, ...) {
+  definition = function(object, do, ..., f = NULL) {
     n <- length(object)
     hat <- do(object, ...)
 
@@ -199,8 +200,9 @@ setMethod(
       FUN.VALUE = double(1),
       x = object, do = do, ...
     )
-    values <- summary_jackknife(values, hat)
-    values
+
+    if (is.function(f)) return(f(values))
+    summary_jackknife(values, hat)
   }
 )
 

R/utilities.R

@@ -1,5 +1,43 @@
 # HELPERS
 
+## https://stackoverflow.com/questions/56191862/where-do-i-specify-random-seed-for-tests-in-r-package
+#' Evaluate an Expression with a Temporarily Seed
+#'
+#' @param expr An [`expression`] to be evaluated.
+#' @param seed A single value to be passed to [set.seed()].
+#' @param envir	The [environment][environment()] in which `expr` should be
+#'  evaluated.
+#' @param rounding A [`logical`] scalar: should the default discrete uniform
+#'  generation method in \R versions prior to 3.6.0 be used? Usefull for unit
+#'  testing.
+#' @param ... Further arguments to be passed to [set.seed()].
+#' @return
+#'  The results of `expr` evaluated.
+#' @seealso [set.seed()]
+#' @keywords internal
+with_seed <- function(expr, seed, ..., envir = parent.frame(), rounding = TRUE) {
+  expr <- substitute(expr)
+  ## Save and restore the random number generator (RNG) state
+  env <- globalenv()
+  old_seed <- env$.Random.seed
+  on.exit({
+    if (is.null(old_seed)) {
+      rm(list = ".Random.seed", envir = env, inherits = FALSE)
+    } else {
+      assign(".Random.seed", value = old_seed, envir = env, inherits = FALSE)
+    }
+  })
+  ## Keep the results the same for R versions prior to 3.6
+  if (isTRUE(rounding) && getRversion() >= "3.6") {
+    ## Set sample.kind = "Rounding" to reproduce the old sampling
+    ## Suppress warning "non-uniform 'Rounding' sampler used"
+    suppressWarnings(set.seed(seed, sample.kind = "Rounding"))
+  } else {
+    set.seed(seed)
+  }
+  eval(expr, envir = envir)
+}
+
 # Helpers ======================================================================
 #' Helpers
 #'

inst/examples/ex-resample.R

@@ -7,5 +7,11 @@ bootstrap(x, do = mean, n = 100)
 quant <- function(x) { quantile(x, probs = c(0.25, 0.75)) }
 bootstrap(x, n = 100, do = mean, f = quant)
 
+## Get the n bootstrap values
+bootstrap(x, n = 100, do = mean, f = function(x) { x })
+
 ## Jackknife
 jackknife(x, do = mean) # Sample mean
+
+## Get the leave-one-out values instead of summary
+jackknife(x, do = mean, f = function(x) { x })

tests/testthat/_snaps/statistics.md

@@ -61,3 +61,57 @@
       [3,] 0.06413649 0.1558635
       [4,] 0.27661853 0.4133815
 
+# Bootstrap
+
+    Code
+      bootstrap_summary
+    Output
+         original        mean        bias       error 
+       0.07651681  0.05439040 -0.02212641  0.19517006 
+
+---
+
+    Code
+      bootstrap_values
+    Output
+        [1]  0.274803536  0.271761759  0.131918811 -0.108350186 -0.407722040
+        [6]  0.474370950  0.096520144  0.140542686  0.048469812  0.045673600
+       [11]  0.098911690  0.174660229  0.407774218  0.076149892  0.233883643
+       [16] -0.158214771  0.296334314 -0.236746678 -0.027278798  0.140680055
+       [21]  0.084317514  0.131248858 -0.126159227  0.351657449  0.130019851
+       [26] -0.064196657  0.049937514 -0.166241173 -0.085116371  0.003131881
+       [31]  0.159010387  0.126015372  0.175198395  0.127748767 -0.047995525
+       [36]  0.245223426 -0.029654335  0.153690734  0.372854685 -0.061229774
+       [41]  0.102818533  0.037641950  0.125330950 -0.082734254 -0.081537996
+       [46]  0.003646969  0.096332439  0.116993888 -0.176113259  0.345299218
+       [51]  0.595379927 -0.311485051  0.093283378  0.342125925  0.001140435
+       [56]  0.246825067 -0.173706701  0.113309462  0.338231445 -0.425915621
+       [61]  0.446095061  0.156535443 -0.266725129  0.175778042  0.201231074
+       [66] -0.065400705  0.015916353  0.257033280 -0.111331877  0.073891696
+       [71] -0.004623160 -0.326442795  0.048737537  0.005317883 -0.053255377
+       [76]  0.081287636 -0.013228710  0.344409439 -0.125705295 -0.181941858
+       [81]  0.118077490 -0.154583474  0.035114092  0.026580327 -0.119802041
+       [86] -0.197818974  0.107870115  0.253120777 -0.199949361  0.050689044
+       [91]  0.002231045  0.309152956 -0.390878169  0.027390498  0.026426616
+       [96]  0.223192535 -0.022840641 -0.146577129 -0.014807361  0.038401655
+
+# Jackknife
+
+    Code
+      jackknife(x, do = mean)
+    Output
+        original       mean       bias      error 
+      0.07651681 0.07651681 0.00000000 0.18647363 
+
+---
+
+    Code
+      jackknife(x, do = mean, f = function(x) {
+        x
+      })
+    Output
+       [1]  0.04972670  0.04320369  0.08629682  0.10441228  0.04865520  0.17622590
+       [7]  0.04738093  0.09508001  0.09549979  0.12892938  0.08666232 -0.01510399
+      [13]  0.06103728  0.05316420  0.12004569  0.03754928  0.12719441  0.09799546
+      [19]  0.02155913  0.06482171
+

tests/testthat/test-statistics.R

@@ -21,3 +21,19 @@ test_that("Confidence interval for binomial proportions", {
   expect_snapshot(confidence_multinomial(x))
   expect_snapshot(confidence_multinomial(x, corrected = TRUE))
 })
+test_that("Bootstrap", {
+  bootstrap_summary <- with_seed({
+    bootstrap(rnorm(20), n = 100, do = mean)
+  }, seed = 12345)
+  expect_snapshot(bootstrap_summary)
+
+  bootstrap_values <- with_seed({
+    bootstrap(rnorm(20), n = 100, do = mean, f = function(x) { x })
+  }, seed = 12345)
+  expect_snapshot(bootstrap_values)
+})
+test_that("Jackknife", {
+  x <- with_seed(rnorm(20), seed = 12345)
+  expect_snapshot(jackknife(x, do = mean))
+  expect_snapshot(jackknife(x, do = mean, f = function(x) { x }))
+})