R/fwb.array.R

In fwb: Fractional Weighted Bootstrap

#' Recover Bootstrap Weights
#'
#' @description `fwb.array()` returns the bootstrap weights generated by [fwb()].
#'
#' @inheritParams fwb.ci
#'
#' @returns
#' A matrix with `R` rows and `n` columns, where `R` is the number of bootstrap replications and `n` is the number of observations in `boot.out$data`.
#'
#' @details
#' The original seed is used to recover the bootstrap weights before being reset.
#'
#' Bootstrap weights are used in computing BCa confidence intervals by approximating the empirical influence function for each unit with respect to each parameter (see Examples).
#'
#' @seealso
#' * [fwb()] for performing the fractional weighted bootstrap
#' * \pkgfun{boot}{boot.array} for the equivalent function in \pkg{boot}
#'
#' See `vignette("fwb-rep")` for information on replicability.
#'
#'
#' @examples
#' set.seed(123, "L'Ecuyer-CMRG")
#' data("infert")
#'
#' fit_fun <- function(data, w) {
#'   fit <- glm(case ~ spontaneous + induced, data = data,
#'              family = "quasibinomial", weights = w)
#'   coef(fit)
#' }
#'
#' fwb_out <- fwb(infert, fit_fun, R = 300,
#'                verbose = FALSE)
#'
#' fwb_weights <- fwb.array(fwb_out)
#'
#' dim(fwb_weights)
#'
#' # Recover computed estimates:
#' est1 <- fit_fun(infert, fwb_weights[1, ])
#'
#' stopifnot(all.equal(est1, fwb_out$t[1, ]))
#'
#' # Compute empirical influence function:
#' empinf <- lm.fit(x = fwb_weights / ncol(fwb_weights),
#'                  y = fwb_out$t)$coefficients
#'
#' empinf <- sweep(empinf, 2L, colMeans(empinf))

#' @export
fwb.array <- function(fwb.out) {
  arg::arg_supplied(fwb.out)
  arg::arg_is(fwb.out, "boot")

  if (identical(.attr(fwb.out, "boot_type"), "boot")) {
    rlang::check_installed("boot")
    return(boot::boot.array(fwb.out))
  }

  arg::arg_is(fwb.out, "fwb")

  gen_weights <- make_gen_weights(fwb.out[["wtype"]])

  n <- nrow(fwb.out[["data"]])
  R <- fwb.out[["R"]]

  if (!isTRUE(.attr(fwb.out, "simple")) || is_null(.attr(fwb.out, "cl"))) {
    with_seed_preserved({
      return(gen_weights(n, R, fwb.out[["strata"]]))
    }, new_seed = fwb.out[["seed"]])
  }

  if (isTRUE(.attr(fwb.out, "simple")) &&
      isTRUE(.attr(fwb.out, "random_statistic"))) {
    arg::wrn('bootstrap weights cannot be reliably re-generated when there is randomness in {.arg statistic} and {.code simple = TRUE} in the call to {.fun fbw}. See {.vignette [vignette("fwb-rep")]{fwb::fwb-rep} for details')
  }

  FUN <- function(i) {
    drop(gen_weights(n, 1L, fwb.out[["strata"]]))
  }

  opb <- pbapply::pboptions(type = "none")
  on.exit(pbapply::pboptions(opb))

  #Run bootstrap
  with_seed_preserved({
    if (identical(.attr(fwb.out, "cl"), "future"))
      do.call("rbind", pbapply::pblapply(seq_len(R), FUN, cl = "future", future.seed = TRUE))
    else
      do.call("rbind", pbapply::pblapply(seq_len(R), FUN, cl = .attr(fwb.out, "cl")))
  }, new_seed = fwb.out[["seed"]])
}