There are several "flavours" or forms of the bootstrap (e.g. non-parametric, parametric, residual resampling and many more). The bootstrap in the example is called a non-parametric bootstrap, or case resampling (see here, here, here and here for applications in regression). The basic idea is that you treat your sample as population and repeatedly draw new samples from it with replacement. All original observations have equal probability of being drawn into the new sample. After drawing the new samples from your sample, you calculate the statistic of interest, this may be a mean, a median or regression coefficients. This is repeated $n$ times. In each iteration, some observations from your original sample are drawn multiple times while some observations may not be sampled at all. After the iterations, the mean, median and standard deviations of the $n$ bootstrap-estimates are calculated.
Bootstrapping is often used for:
- Calculation of confidence intervals (and estimation of the standard errors)
- Estimation of the bias of the point estimates
There are several methods for calculating confidence intervals based on the bootstrap samples (this paper provides explanation and guidance). One very simple method is just calculating the 2.5 and 97.5% percentiles of the bootstrap samples (see code below).
The bias is simply estimated as the difference between the bootstrap samples and the original estimates.
Let's replicate the example from the website but using our own loop:
#-----------------------------------------------------------------------------
# Load packages
#-----------------------------------------------------------------------------
require(ggplot2)
require(pscl)
require(MASS)
require(boot)
#-----------------------------------------------------------------------------
# Load data
#-----------------------------------------------------------------------------
zinb <- read.csv("http://www.ats.ucla.edu/stat/data/fish.csv")
zinb <- within(zinb, {
nofish <- factor(nofish)
livebait <- factor(livebait)
camper <- factor(camper)
})
#-----------------------------------------------------------------------------
# Calculate zero-inflated regression
#-----------------------------------------------------------------------------
m1 <- zeroinfl(count ~ child + camper | persons, data = zinb,
dist = "negbin", EM = TRUE)
#-----------------------------------------------------------------------------
# Store the original regression coefficients
#-----------------------------------------------------------------------------
original.estimates <- as.vector(t(do.call(rbind, coef(summary(m1)))[, 1:2]))
#-----------------------------------------------------------------------------
# Set the number of replications
#-----------------------------------------------------------------------------
n.sim <- 2000
#-----------------------------------------------------------------------------
# Set up a matrix to store the results
#-----------------------------------------------------------------------------
store.matrix <- matrix(NA, nrow=n.sim, ncol=12)
#-----------------------------------------------------------------------------
# The loop
#-----------------------------------------------------------------------------
set.seed(123)
for(i in 1:n.sim) {
#-----------------------------------------------------------------------------
# Draw the observations WITH replacement
#-----------------------------------------------------------------------------
data.new <- zinb[sample(1:dim(zinb)[1], dim(zinb)[1], replace=TRUE),]
#-----------------------------------------------------------------------------
# Calculate the model with this "new" data
#-----------------------------------------------------------------------------
m <- zeroinfl(count ~ child + camper | persons,
data = data.new, dist = "negbin",
start = list(count = c(1.3711, -1.5152, 0.879),
zero = c(1.6028, -1.6663)))
#-----------------------------------------------------------------------------
# Store the results
#-----------------------------------------------------------------------------
store.matrix[i, ] <- as.vector(t(do.call(rbind, coef(summary(m)))[, 1:2]))
}
#-----------------------------------------------------------------------------
# Save the means, medians and SDs of the bootstrapped statistics
#-----------------------------------------------------------------------------
boot.means <- colMeans(store.matrix, na.rm=T)
boot.medians <- apply(store.matrix,2,median, na.rm=T)
boot.sds <- apply(store.matrix,2,sd, na.rm=T)
#-----------------------------------------------------------------------------
# The bootstrap bias is the difference between the mean bootstrap estimates
# and the original estimates
#-----------------------------------------------------------------------------
boot.bias <- colMeans(store.matrix, na.rm=T) - original.estimates
#-----------------------------------------------------------------------------
# Quantile confidence intervals
#-----------------------------------------------------------------------------
conf.mat <- matrix(apply(store.matrix, 2 ,quantile, c(0.025, 0.975), na.rm=T),ncol=2)
colnames(conf.mat) <- c("95%-CI Lower", "95%-CI Upper")
#-----------------------------------------------------------------------------
# Set up summary data frame
#-----------------------------------------------------------------------------
summary.frame <- data.frame(mean=boot.means, median=boot.medians, sd=boot.sds, bias=boot.bias, "CI_lower"=conf.mat[,1], "CI_upper"=conf.mat[,2])
summary.frame
mean median sd bias CI_lower CI_upper
1 1.2997903 1.3013079 0.39673804 -0.0712912250 0.51960031 -1.3066351
2 0.2526688 0.2485721 0.03207790 -0.0034460509 2.06049015 -0.4379761
3 -1.5661686 -1.5571989 0.26220013 -0.0509238534 0.19897759 0.1449354
4 0.2004957 0.1986454 0.01949437 0.0049018955 0.32294130 0.2140118
5 0.9543834 0.9252035 0.48914523 0.0753405152 -2.12899783 0.4414957
6 0.2702361 0.2688095 0.02042504 0.0009582592 -1.09195998 8.0471143
7 -0.8996836 -0.9082379 0.22173921 0.0856792712 0.16744198 0.5811394
8 0.1788934 0.1781161 0.01667250 0.0029512565 0.24179778 57.6416680
9 2.0682618 1.7718903 1.59102322 0.4654897865 0.03492612 -8.4344038
10 4.0208757 0.8269775 13.23433539 3.1845709680 1.90245108 -1.1156127
11 -2.0969405 -1.6717102 1.56310762 -0.4306843771 0.23271566 0.3363101
12 3.8660345 0.6434859 13.27525033 3.1870641941 0.31370593 57.6061993
Some explanations
- The difference between the mean of the bootstrap estimates and the original estimates is what is called "bias" in the output of
boot - What the output of
bootcalls "std. error" is the standard deviation of the bootstrapped estimates
Compare it with the output from boot:
#-----------------------------------------------------------------------------
# Compare with boot output
#-----------------------------------------------------------------------------
set.seed(10)
res <- boot(zinb, f, R = 2000, parallel = "snow", ncpus = 4)
res
Bootstrap Statistics :
original bias std. error
t1* 1.3710504 -0.076735010 0.39842905
t2* 0.2561136 -0.003127401 0.03172301
t3* -1.5152609 -0.064110745 0.26554358
t4* 0.1955916 0.005819378 0.01933571
t5* 0.8790522 0.083866901 0.49476780
t6* 0.2692734 0.001475496 0.01957823
t7* -0.9853566 0.083186595 0.22384444
t8* 0.1759504 0.002507872 0.01648298
t9* 1.6031354 0.482973831 1.58603356
t10* 0.8365225 3.240981223 13.86307093
t11* -1.6665917 -0.453059768 1.55143344
t12* 0.6793077 3.247826469 13.90167954
#-----------------------------------------------------------------------------
# Our summary table
#-----------------------------------------------------------------------------
summary.frame
mean median sd bias CI_lower CI_upper
1 1.2997903 1.3013079 0.39673804 -0.0712912250 0.51960031 -1.3066351
2 0.2526688 0.2485721 0.03207790 -0.0034460509 2.06049015 -0.4379761
3 -1.5661686 -1.5571989 0.26220013 -0.0509238534 0.19897759 0.1449354
4 0.2004957 0.1986454 0.01949437 0.0049018955 0.32294130 0.2140118
5 0.9543834 0.9252035 0.48914523 0.0753405152 -2.12899783 0.4414957
6 0.2702361 0.2688095 0.02042504 0.0009582592 -1.09195998 8.0471143
7 -0.8996836 -0.9082379 0.22173921 0.0856792712 0.16744198 0.5811394
8 0.1788934 0.1781161 0.01667250 0.0029512565 0.24179778 57.6416680
9 2.0682618 1.7718903 1.59102322 0.4654897865 0.03492612 -8.4344038
10 4.0208757 0.8269775 13.23433539 3.1845709680 1.90245108 -1.1156127
11 -2.0969405 -1.6717102 1.56310762 -0.4306843771 0.23271566 0.3363101
12 3.8660345 0.6434859 13.27525033 3.1870641941 0.31370593 57.6061993
Compare the "bias" columns and the "std. error" with the "sd" column of our own summary table.