Bootstrap

From 太極
Jump to navigation Jump to search

General

  • Bootstrap from Wikipedia.
    • This contains an overview of different methods for computing bootstrap confidence intervals.
    • boot.ci() from the 'boot' package provides a short explanation for different methods for computing bootstrap confidence intervals.
    • Typically not useful for correlated (dependent) data.
  • bootstrap package. "An Introduction to the Bootstrap" by B. Efron and R. Tibshirani, 1993

Nonparametric bootstrap

This is the most common bootstrap method

The upstrap Crainiceanu & Crainiceanu, Biostatistics 2018

Parametric bootstrap

  • Parametric bootstraps resample a known distribution function, whose parameters are estimated from your sample
  • In small samples, a parametric bootstrap approach might be preferred according to the wikipedia
  • http://www.math.ntu.edu.tw/~hchen/teaching/LargeSample/notes/notebootstrap.pdf#page=3 No package is used
  • A parametric or non-parametric bootstrap?
  • https://www.stat.cmu.edu/~cshalizi/402/lectures/08-bootstrap/lecture-08.pdf#page=11
  • simulatorZ Bioc package
  • Bootstrap R tutorial : Learn about parametric and non-parametric bootstrap through simulation 
    df <- scan(textConnection("3331.608 3549.913 2809.252 2671.465 3383.177 3945.721 3672.601 3922.416 4647.278 4088.246 3718.874 3724.443 3925.457 3112.920 3495.621 3651.779 3240.194 3867.347 3431.015 4163.725"), sep = " ")
sd(df)
# [1] 464.5509
mean(df)
# [1] 3617.653
set.seed(1)
theta_star <- vector() 
for (i in 1:1000){ theta_star[i] <- mean(rnorm(length(df),mean = 3617.7 ,sd = 464.6)) }
theta_boot <- mean(theta_star) # bootstrap estimate of theta_hat theta_boot
theta_boot
# [1] 3615.208
boot_se <- sd(theta_star) # standard eorrs of the estimate boot_se
boot_se
# [1] 100.106

bias <- theta_boot - mean(df) 
bias
# -2.444507
CI <-c(theta_boot-1.96*boot_se,theta_boot +1.96*boot_se) 
CI
# 3419.000 3811.416

    Using the boot package 

    # tell the boot() which distribution random sample will be generated
gen_function <- function(x,mle) { 
  # Input:
  #   x: observed data
  #   mle: parameter estimates
  data <- rnorm(length(x),mle[1],mle[2]) 
  return(data)
}

# function to calculate sample statistic 
theta_star_function <- function(x,i) mean(x[i])

set.seed(1)
B <- boot(data = df, 
          sim = "parametric", 
          ran.gen = gen_function, 
          mle = c(3617.7,464.6), 
          statistic = theta_star_function, 
          R=1000) 
B
# PARAMETRIC BOOTSTRAP
# Call:
# boot(data = df, statistic = theta_star_function, R = 1000, sim = "parametric", 
#     ran.gen = gen_function, mle = c(3617.7, 464.6))
# 
# Bootstrap Statistics :
#     original    bias    std. error
# t1* 3617.653 -2.444507     100.106
plot(B)

Examples

Standard error

  • Standard error from a mean 
    foo <- function() mean(sample(x, replace = TRUE))
set.seed(1234)
x <- rnorm(300)
set.seed(1)
sd(replicate(10000, foo()))
# [1] 0.05717679
sd(x)/sqrt(length(x)) # The se of mean is s/sqrt(n)
# [1] 0.05798401

set.seed(1234)
x <- rpois(300, 2)
set.seed(1)
sd(replicate(10000, foo()))
# [1] 0.08038607
sd(x)/sqrt(length(x)) # The se of mean is s/sqrt(n)
# [1] 0.08183151
  • Difference of means from two samples (cf 8.3 The two-sample problem from the book "An introduction to Bootstrap" by Efron & Tibshirani) 
    # Define the two samples
sample1 <- 1:10
sample2 <- 11:20

# Define the number of bootstrap replicates
nboot <- 100000

# Initialize a vector to store the bootstrap estimates
boot_estimates <- numeric(nboot)

# Run the bootstrap
set.seed(123)
for (i in seq_len(nboot)) {
  # Resample the data with replacement
  resample1 <- sample(sample1, replace = TRUE)
  resample2 <- sample(sample2, replace = TRUE)
  
  # Compute the difference of means
  boot_estimates[i] <- mean(resample1) - mean(resample2)
}

# Compute the standard error of the bootstrap estimates
se_boot <- sd(boot_estimates)

# Print the result
cat("Bootstrap SE estimate of difference of means:", se_boot, "\n")
# 1.283541

sd1 <- sd(sample1)
sd2 <- sd(sample2)

# Calculate the sample sizes
n1 <- length(sample1)
n2 <- length(sample2)

# Calculate the true standard error of the difference of means
se_true <- sqrt((sd1^2/n1) + (sd2^2/n2))

# Print the result
cat("True SE of difference of means:", se_true, "\n") \
# 1.354006
  • Difference of means from two samples, another example. It is generally observed that the bootstrap estimate of the standard error (SE) can be more conservative compared to the true SE. See Bootstrap Standard Error Estimates and Inference. 
    library(boot)

# Set a seed for reproducibility
set.seed(123)

# Generate two samples from normal distribution
sample1 <- rnorm(100, mean = 0, sd = 1)
sample2 <- rnorm(100, mean = 1, sd = 1)

# Define a function to compute the difference in means
diff_in_means <- function(data, indices) {
  sample1 <- data$sample1[indices]
  sample2 <- data$sample2[indices]
  # Return the difference in means
  return(mean(sample1) - mean(sample2))
}

# Combine the samples into a list
samples <- list(sample1 = sample1, sample2 = sample2)

# Use the boot() function to perform bootstrap
results <- boot(data = samples, statistic = diff_in_means, R = 1000)

# Print the bootstrap results
print(results)
# Bootstrap Statistics :
#      original       bias    std. error
# t1* -1.168565 -0.001910976   0.2193785

# Compute the true standard error
true_se <- sqrt(var(sample1)/length(sample1) + var(sample2)/length(sample2))

# Print the true standard error
print(paste("True SE: ", true_se))
# [1] "True SE:  0.132977288582124"
  • SE of the sample median 
    x <- c(1, 5, 8, 3, 7)
median(x) # 5
set.seed(1)
foo <- function(x) median(sample(x, replace = T))
y <- replicate(100, foo(x=x))
sd(y)
# 1.917595
c(median(x)-1.96*sd(y), median(x)+1.96*sd(y))
# [1] 1.241513 8.758487
  • How to Perform Bootstrapping in R (With Examples). boot package was used to estimate the SE of R-squared.
  • SE for sample relative risk

Bootstrapping Extreme Value Estimators

Bootstrapping Extreme Value Estimators de Haan, 2022

R-squared from a regression

R in Action Nonparametric bootstrapping 

# Compute the bootstrapped 95% confidence interval for R-squared in the linear regression
rsq <- function(data, indices, formula) {
  d <- data[indices,] # allows boot to select sample
  fit <- lm(formula, data=d)
  return(summary(fit)$r.square)
} # 'formula' is optional depends on the problem

# bootstrapping with 1000 replications
set.seed(1234)
bootobject <- boot(data=mtcars, statistic=rsq, R=1000, 
                   formula=mpg~wt+disp)
plot(bootobject) # or plot(bootobject, index = 1) if we have multiple statistics
ci <- boot.ci(bootobject, conf = .95, type=c("perc", "bca") ) 
    # default type is "all" which contains c("norm","basic", "stud", "perc", "bca"). 
    # 'bca' (Bias Corrected and Accelerated) by Efron 1987 uses 
    # percentiles but adjusted to account for bias and skewness.
# Level     Percentile            BCa          
# 95%   ( 0.6838,  0.8833 )   ( 0.6344,  0.8549 )
# Calculations and Intervals on Original Scale
# Some BCa intervals may be unstable
ci$bca[4:5]  
# [1] 0.6343589 0.8549305
# the mean is not the same
mean(c(0.6838,  0.8833 ))
# [1] 0.78355
mean(c(0.6344,  0.8549 ))
# [1] 0.74465
summary(lm(mpg~wt+disp, data = mtcars))$r.square
# [1] 0.7809306

Coefficients in linear regression

Coefficients in a multiple linear regression

Correlation

A Practical Guide to Bootstrap in R

Stability of clinical prediction models

Stability of clinical prediction models developed using statistical or machine learning methods 2023

Retrieved from "https://wiki.taichimd.us/index.php?title=Bootstrap&oldid=38098"