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
