bcaboot: Perform BCa Bootstrap for a Given Statistic

In ddarmon/clp: R Package for Computing Confidence Functions to Accompany CLP

Description

Approximate the bootstrap distribution and estimate bias and acceleration parameters for the BCa bootstrap, given a statistic as would be supplied to boot from the boot package.

Usage

bcaboot(
  data,
  statistic,
  B = 2000,
  sim = "ordinary",
  stratified = FALSE,
  ran.gen = function(d, p) d,
  mle = NULL,
  formula = NULL
)

Arguments

data	the dataframe or data matrix
statistic	a function that computes the statistic from data
B	the number of bootstrap replicates
sim	either "ordinary" (for case resampling bootstrap) or "parametric" (for parametric bootstrap)
stratified	whether or not (default) to use the use stratified sampling for the bootstrapping.
ran.gen	a function returning a random sample, for the parametric bootstrap
mle	the maximum likelihood estimate from the original sample, for the parametric bootstrap
formula	for use with interfacing to lm, glm, etc.

Value

A list containing the sample estimate, bootstrap estimates, bootstrap CDF, bias, and acceleration parameters for the BCa bootstrap.

References

Bradley Efron. "Better bootstrap confidence intervals." Journal of the American Statistical Association 82.397 (1987): 171-185.

Thomas J. DiCiccio and Bradley Efron. "Bootstrap confidence intervals." Statistical Science (1996): 189-212.

Examples

# Bootstrap confidence functions for a single mean.

t.one.sample <- function(data, id = 1:length(data)){
  dat <- data[id]

  d <- mean(dat)

  return(d)
}

data(dietstudy)

bc <- bcaboot(data = dietstudy$weightchange[dietstudy$diet == 'Low Carb'],
        statistic = t.one.sample,
        B = 2000)

# Reproduce BCa confidence density from Figure 11.7
# of *Computer Age Statistical Inference*.

scor <-
read.table('https://web.stanford.edu/~hastie/CASI_files/DATA/student_score.txt',
header = TRUE)

statistic <- function(data, id = 1:nrow(data)){
  dat <- data[id, ]

  Sigma <- cov(dat)*((nrow(dat)-1)/nrow(dat))

  lams <- eigen(Sigma, symmetric = TRUE, only.values = TRUE)$values

  return(lams[1])
}

ran.gen <- function(data, mle){
  MASS::mvrnorm(n = nrow(data), mle$mu, mle$Sigma)
}

bc <- bcaboot(data = scor, statistic = statistic, B = 8000, sim = "parametric", ran.gen = ran.gen,
              mle = list(mu = colMeans(scor),
                         Sigma = cov(scor)*((nrow(scor)-1)/nrow(scor))))

ddarmon/clp documentation built on Jan. 25, 2021, 6:22 p.m.