10
$\begingroup$

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.

enter image description here

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.

$\endgroup$
8
  • $\begingroup$ Is this a homework problem? $\endgroup$
    – Jon
    Commented Dec 21, 2016 at 19:54
  • 3
    $\begingroup$ Here's a good paper to read, "Theoretical Comparison of Bootstrap Confidence Intervals": projecteuclid.org/download/pdf_1/euclid.aos/1176350933 $\endgroup$
    – Jon
    Commented Dec 21, 2016 at 19:55
  • 2
    $\begingroup$ 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. $\endgroup$ Commented Dec 21, 2016 at 20:25
  • 1
    $\begingroup$ Did you consider BC and BCa? $\endgroup$ Commented Dec 21, 2016 at 20:27
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Dec 21, 2016 at 20:38

2 Answers 2

5
$\begingroup$

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.
Nerve dataset downloaded from https://www.stat.cmu.edu/~larry/all-of-nonpar/data.html

$\endgroup$
4
$\begingroup$

As Michael Chernick notes, it would be useful to also look at the bias-corrected (BC) and bias-corrected and accelerated (BCa) bootstrap.

The BCa variant in particular attempts to deal with skewness in data, as you apparently have. DiCiccio & Efron (1996, Statistical Science) found that it performs well, as do Davison & Hinkley, Bootstrap Methods and their Applications (1997).

Why does my bootstrap interval have terrible coverage? is related, and I would especially recommend the article by Canto et al. (2006) that I cite there. And in the end, I agree that the answer is likely related to sample size, as well as your underlying distribution, and the pivotality or non of the statistic you are bootstrapping.

$\endgroup$
1
  • $\begingroup$ It is very important to simulate the non-coverage probabilities on both sides of the confidence interval. When I did that for the log-normal distribution, all bootstrap intervals have terrible coverage, except for bootstrap t which I didn't try. It would also be worth using a standard bootstrap package to check your bootstrap t result. $\endgroup$ Commented Nov 7, 2017 at 13:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.