# Which Bootstrap method is most preferred?

Maybe this question depends on the given data, but is there a "better" bootstrap method than the others? I'm simply using a one variable dataset (that consists of the differences between football scores (2 teams) over the past 15 weeks)..

First note the right skew of this data, I feel like this will take into consideration which bootstrap I would recommend as being "better" or most accurate for the representation of the data.

First here is the standard bootstrap interval

N <- 10^4
n <- length(Differences)
Differences.mean <- numeric(N)
for(i in 1:N)
{
x <- sample(Differences, n, replace = TRUE)
Differences.mean[i]<- mean(x)
}

lower = mean(Differences.mean)-1.96*sd(Differences.mean) #Lower CI
upper = mean(Differences.mean)+1.96*sd(Differences.mean) #Upper CI
= (8.875, 10.916)

mean(Differences.mean)-m  #The bias is fairly small also
= -.0019


Here is a bootstrap percentile interval

 quantile(Differences.mean,c(.025,.975)
= (8.893, 10.938)


Lastly here is the Bootstrap T interval

Tstar = numeric(N)
for(i in 1:N)
{
y =sample(Differences, size = n, replace = TRUE)
Tstar[i] = (mean(y)-m) / (sd(y)/sqrt(n))
}
q1 = quantile(Tstar,.025) #empirical quantiles for bootstrap t (lower)
q2 = quantile(Tstar,.975) #empirical quantiles for bootstrap t (upper)

mean(Differences)-(q2*sd(Differences/sqrt(n)))
mean(Differences)-(q1*sd(Differences/sqrt(n)))
= (8.925, 10.997)


In addition even the t confidence interval seems pretty accurate

 t.test(Differences, conf.level = .95, alternative = "two.sided")
= (8.867, 10.928)


My conclusion would be to choose the bootstrap t interval because it reflects the right skew of the data, it is stretched further to the right than any of the others. My sample size is 224. I think sample size plays a huge role in my conclusion, but my initial question was "is there a better bootstrap method than the others?".. Maybe it really does depend on the data and sample size. Hopefully this isn't too broad.

• Is this a homework problem?
– Jon
Dec 21, 2016 at 19:54
– Jon
Dec 21, 2016 at 19:55
• It looks like you have some good ideas. Look at Tim Hesterberg's work on bootstrap t. I don't quite understand your graph because it does not show how far on the negative side the curve goes. I am a little surprised at how close the intervals are for all methods. I would have expected the standard bootstrap to do much worse. Dec 21, 2016 at 20:25
• Did you consider BC and BCa? Dec 21, 2016 at 20:27
• You could look at these books by authors Efron and Tibshirani 1993. Davison and Hinkley 1997 and Chernick 2007. My book was published by Wiley. Efron Chapman and Hall and Davison and Hinkley Cambridge University Press I think. Dec 21, 2016 at 20:38

It's never too late to fix an error... At least I think there is an error in the implementation of the bootstrap $$t$$ interval, which requires a bootstrap within the bootstrap (double bootstrap).

Note: I refactor the R code for clarity. And since the OP didn't provide data or the code to reproduce the data, I use the nerve data from the book All of Nonparametric Statistics by L. Wasserman. The measurements are 799 waiting times between successive pulses along a nerve fiber and was originally reported in Cox and Lewis (1966).

The goal is to compute the studentized bootstrap interval defined as: \begin{aligned} \left( T - t^*_{1-\alpha/2}\widehat{\operatorname{se}}_{\text{boot}}, T - t^*_{\alpha/2}\widehat{\operatorname{se}}_{\text{boot}} \right) \end{aligned} where

• $$T$$ is the statistic of interest (here skewness);
• $$\widehat{\operatorname{se}}_{\text{boot}}$$ is the standard error of the bootstrapped statistics $$\tilde{T}_1, \ldots, \tilde{T}_B$$;
• $$t^*_q$$ is the $$q$$ sample quantile of $$T^*_1, \ldots, T^*_B$$; and
• $$T^*_b$$ is the bootstrapped pivotal quantity: \begin{aligned} T^*_b = \frac{\tilde{T}_b - T}{\widehat{\operatorname{se}}^*_b} \end{aligned}

The standard error $$\widehat{\operatorname{se}}^*_b$$ of $$T^*_b$$ is estimated with an inner bootstrap inside each bootstrap iteration $$b$$. That's computationally expensive.

set.seed(1234)

# Dataset used in the book All of Nonparametric Statistics by L. Wasserman.
# https://www.stat.cmu.edu/~larry/all-of-nonpar/data.html
x <- scan("nerve.dat")
n <- length(x)


The bootstrap is a general procedure, so let's make the implementation general as well by defining a function estimator to calculate a statistic $$T$$ (for the nerve data it's the skewness) from a sample $$x$$ as well as a function simulator to generate bootstrap resamples of $$x$$.

estimator <- function(x) {
# sample skewness
mean((x - mean(x))^3) / sd(x)^3
}
simulator <- function(x) {
sample(x, size = length(x), replace = TRUE)
}


So here is the bootstrap $$t$$ implementation outlined in the question. I've highlighted the errors in comments.

alpha <- 0.05
B <- 2000

# Warning: This code snippet doesn't implement the bootstrap t correctly.
Tstat <- estimator(x)
Tboot <- numeric(B)
Tstar <- numeric(B)

for (i in seq(B)) {
boot <- simulator(x)
Tboot[i] <- estimator(boot)
# > Error: Doesn't bootstrap the std. error of Tboot[i] correctly.
Tstar[i] <- (Tboot[i] - Tstat) / sd(boot)
}

q.lower <- quantile(Tstar, alpha / 2)
q.upper <- quantile(Tstar, 1 - alpha / 2)

# > Error: Uses the standard deviation of the sample x
#    instead of the standard error of the statistics T
c(
Tstat - q.upper * sd(x),
Tstat - q.lower * sd(x)
)
#>    97.5%     2.5%
#> 1.451426 2.102211


And here is the correct way to do the bootstrap $$t$$ method. Note that this procedure requires to do bootstrap inside the bootstrap. This is more computationally intensive but (we expect) more accurate.

Tstat <- estimator(x)
Tboot <- numeric(B)
Tstar <- numeric(B)

for (i in seq(B)) {
boot <- simulator(x)
Tboot[i] <- estimator(boot)
# > Bootstrap the bootstrap sample to estimate the std. error of Tboot.
boot_of_boot <- replicate(B, estimator(simulator(boot)))
Tstar[i] <- (Tboot[i] - Tstat) / sd(boot_of_boot)
}

q.lower <- quantile(Tstar, alpha / 2)
q.upper <- quantile(Tstar, 1 - alpha / 2)

# > Use the bootstrap estimate of the std. error of T.
c(
Tstat - q.upper * sd(Tboot),
Tstat - q.lower * sd(Tboot)
)
#>    97.5%     2.5%
#> 1.463658 2.292374


Example 3.17 of All of Nonparametric Statistics reports the studentized 95% interval for the skewness of the nerve data as (1.45, 2.28). This agrees well with the result obtain with the second/correct implementation of the bootstrap $$t$$ above.

References

L. Wasserman (2007). All of Nonparametric Statistics. Springer.
Cox, D. and Lewis, P. (1966). The Statistical Analysis of Series of Events. Chapman and Hall.