I've been looking into the boot package in R and while I have found a number of good primers on how to use it, I have yet to find anything that describes exactly what is happening "behind the scenes". For instance, in this example, the guide shows how to use standard regression coefficients as a starting point for a bootstrap regression but doesn't explain what the bootstrap procedure is actually doing to derive the bootstrap regression coefficients. It appears there is some sort of iterative process that is happening but I can't seem to figure out exactly what is going on.
$\begingroup$
$\endgroup$
1
-
3$\begingroup$ It's a long time since I last opened it so I don't know if it will answer your question but the boot package is based in particular on the methods detailed in Davison, A. C., & Hinkley, D. V. (1997). Bootstrap methods and their application. Cambridge: Cambridge University Press. (The reference is also cited in the help file for the boot package.) $\endgroup$– GalaCommented Jul 8, 2013 at 14:32
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
4
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. Then you calculate and store the statistic(s) of interest, this may be the mean, the median or regression coefficients using the newly drawn sample. This is repeated $n$ times. In each iteration, some observations from your original sample are drawn multiple times while some observations may not be drawn at all. After $n$ iterations, you have $n$ stored bootstrap estimates of the statistic(s) of interest (e.g. if $n=1000$ and the statistic of interest is the mean, you have 1000 bootstrapped estimates of the mean). Lastly, summary statistics such as the mean, median and the standard deviation 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 for calculating a 95%-confidence interval is just calculating the empirical 2.5th and 97.5th percentiles of the bootstrap samples (this interval is called the bootstrap percentile interval; see code below). The simple percentile interval method is rarely used in practice as there are better methods, such as the bias-corrected and accelerated bootstrap (BCa). BCa intervals adjust for both bias and skewness in the bootstrap distribution.
The bias is simply estimated as the difference between the mean of the $n$ stored bootstrap samples and the original estimate(s).
Let's replicate the example from the website but using our own loop incorporating the ideas I've outlined above (drawing repeatedly with replacement):
#-----------------------------------------------------------------------------
# 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
#-----------------------------------------------------------------------------
# Basic bootstrap CIs based on the empirical quantiles
#-----------------------------------------------------------------------------
conf.mat <- matrix(apply(store.matrix, 2 ,quantile, c(0.025, 0.975), na.rm=T),
ncol=2, byrow=TRUE)
colnames(conf.mat) <- c("95%-CI Lower", "95%-CI Upper")
And here is our summary table:
#-----------------------------------------------------------------------------
# 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.2998 1.3013 0.39674 -0.0712912 0.51960 2.0605
2 0.2527 0.2486 0.03208 -0.0034461 0.19898 0.3229
3 -1.5662 -1.5572 0.26220 -0.0509239 -2.12900 -1.0920
4 0.2005 0.1986 0.01949 0.0049019 0.16744 0.2418
5 0.9544 0.9252 0.48915 0.0753405 0.03493 1.9025
6 0.2702 0.2688 0.02043 0.0009583 0.23272 0.3137
7 -0.8997 -0.9082 0.22174 0.0856793 -1.30664 -0.4380
8 0.1789 0.1781 0.01667 0.0029513 0.14494 0.2140
9 2.0683 1.7719 1.59102 0.4654898 0.44150 8.0471
10 4.0209 0.8270 13.23434 3.1845710 0.58114 57.6417
11 -2.0969 -1.6717 1.56311 -0.4306844 -8.43440 -1.1156
12 3.8660 0.6435 13.27525 3.1870642 0.33631 57.6062
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 and confidence intervals
#-----------------------------------------------------------------------------
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
perc.cis <- matrix(NA, nrow=dim(res$t)[2], ncol=2)
for( i in 1:dim(res$t)[2] ) {
perc.cis[i,] <- boot.ci(res, conf=0.95, type="perc", index=i)$percent[4:5]
}
colnames(perc.cis) <- c("95%-CI Lower", "95%-CI Upper")
perc.cis
95%-CI Lower 95%-CI Upper
[1,] 0.52240 2.1035
[2,] 0.19984 0.3220
[3,] -2.12820 -1.1012
[4,] 0.16754 0.2430
[5,] 0.04817 1.9084
[6,] 0.23401 0.3124
[7,] -1.29964 -0.4314
[8,] 0.14517 0.2149
[9,] 0.29993 8.0463
[10,] 0.57248 56.6710
[11,] -8.64798 -1.1088
[12,] 0.33048 56.6702
#-----------------------------------------------------------------------------
# Our summary table
#-----------------------------------------------------------------------------
summary.frame
mean median sd bias CI_lower CI_upper
1 1.2998 1.3013 0.39674 -0.0712912 0.51960 2.0605
2 0.2527 0.2486 0.03208 -0.0034461 0.19898 0.3229
3 -1.5662 -1.5572 0.26220 -0.0509239 -2.12900 -1.0920
4 0.2005 0.1986 0.01949 0.0049019 0.16744 0.2418
5 0.9544 0.9252 0.48915 0.0753405 0.03493 1.9025
6 0.2702 0.2688 0.02043 0.0009583 0.23272 0.3137
7 -0.8997 -0.9082 0.22174 0.0856793 -1.30664 -0.4380
8 0.1789 0.1781 0.01667 0.0029513 0.14494 0.2140
9 2.0683 1.7719 1.59102 0.4654898 0.44150 8.0471
10 4.0209 0.8270 13.23434 3.1845710 0.58114 57.6417
11 -2.0969 -1.6717 1.56311 -0.4306844 -8.43440 -1.1156
12 3.8660 0.6435 13.27525 3.1870642 0.33631 57.6062
Compare the "bias" columns and the "std. error" with the "sd" column of our own summary table. Our 95%-confidence intervals are very similar to the confidence intervals calculated by boot.ci using the percentile method (not all though: look at the lower limit of parameter with index 9).
answered Jul 8, 2013 at 14:40
-
$\begingroup$ Thanks for the detailed reply. Are you basically saying that the coefficients are the average of the 2000 sets of coefficients that were generated? $\endgroup$– zgall1Commented Jul 8, 2013 at 15:02
-
4$\begingroup$ @zgall1 Yes, basically (or the medians). But most often, the point estimates are of less interest in bootstrapping. Often, the reason to use bootstrap is to find confidence intervals that take the sample distribution into account and don't assume any distribution (e.g. normal). In that way, the bootstrapped confidence intervals are often asymmtetric, whereas confidence intervals based on the normal or $t$-distribution are symmetric. $\endgroup$ Commented Jul 8, 2013 at 15:04
-
$\begingroup$ 'The basic idea is that you treat your sample as population and repeatedly draw new samples from it with replacement' - how to determine what is the size of the new samples? $\endgroup$– SinusxCommented Jan 19, 2017 at 19:26
-
1$\begingroup$ @Sinusx Normally, you draw samples of the same size as the original sample. The crucial idea is to draw the sample with replacement. So some values from the original sample will get drawn multiple times and some values not at all. $\endgroup$ Commented Jan 19, 2017 at 22:16
$\begingroup$
$\endgroup$
You should focus on the function that is passed to boot as the "statistic" parameter and notice how it is constructed.
f <- function(data, i) {
require(pscl)
m <- zeroinfl(count ~ child + camper | persons,
data = data[i, ], dist = "negbin",
start = list(count = c(1.3711, -1.5152, 0.879), zero = c(1.6028, -1.6663)))
as.vector(t(do.call(rbind, coef(summary(m)))[, 1:2]))
}
The "data" argument is going to receive an entire data frame, but the "i" argument is going to receive a sample of row indices generated by the "boot" and taken from 1:NROW(data). As you can see from that code, "i" is then used to create a neo-sample which is passed to zeroinl and then only selected portions of it's results are returned.
Let's imagine that "i" is {1,2,3,3,3,6,7,7,10}. The "[" function will return just those rows with 3 copies of row 3 and 2 copies of row 7. That would be the basis for a single zeroinl() calculation and then the coefficients will be returned to boot as the result from that replicate of the process. The number of such replicates is controlled by the "R" parameter.
Since only the regression coefficients are returned from statistic in this case, the boot function will return these accumulated coefficients as the value of "t". Further comparisons can be performed by other boot-package functions.
