Correct creation of the null distribution for bootstrapped $p$-values

Question

Let's say we want to calculate the two-sided $p$-value for a linear regression coefficient using bootstrapping. The null hypothesis is $\beta_1 = 0$. For the non-parametric (case-resampling) bootstrap, we would do the following steps:

1. Calculate the original statistic $T_0$, e.g. the $t$-statistic of the coefficient of interest in the original sample.
2. Sample the observations with replacement to create a bootstrap sample of the same size as the original sample.
3. Calculate the regression on this bootstrap sample and store the bootstrap statistic $T^*$, i.e. the coefficient of interest.
4. Repeat steps 2 and 3 a large number of times, say $R$ (e.g. $R = 1000$). This gives us $R$ bootstrapped statistics $T_r^*$.
5. Shift the distribution of $T_r^*$ to generate the null distribution.
6. Calculate the two-sided $p$-value by calculating: $p_{\text{boot}} = \dfrac{\text{#}\{|T_r^*| \geq |T_0|\} + 1}{R + 1}$.

My question concerns step 5: Should we shift the bootstrapped statistics by $T_r^*-\bar{T}_r^*$ so that its distribution is centered around $0$ or by $T_r^*-T_0$ so that it's centered around at the bias $\hat{B}^*=\bar{T}_r^* - T_0$?

The first version respects the bootstrap analogy that the sample is our population regarding the bootstrap samples. The second version makes sense because it corresponds exactly to the null value specified in the null hypothesis.

Here is a small R script that illustrates the procedures. The $p$-values calculated by the two versions are virtually identical here (the bias is only $0.12$) but I suspect that this is not always the case, especially when the bias is large.

# Load data
data(swiss)
Calculate original model
mod <- lm(Infant.Mortality~., data = swiss)
Store sample size
N <- nobs(mod)
Store original t-statistic for "Fertility"
stat_orig <- summary(mod)$coefficients["Fertility", "t value"]
Start the bootstrap
R <- 10000 # Number of bootstrap samples
Start resampling
set.seed(142857) # Reproducibility
stat_boot <- replicate(R, {
  boot_ind <- sample.int(N, replace = TRUE) # Sample indices with replacement
  tmp_mod <- lm(Infant.Mortality~., data = swiss[boot_ind, ])
  summary(tmp_mod)$coefficients["Fertility", "t value"]
})
Create null distribution by shifting bootstrap distribution
stat_boot_1 <- stat_boot - mean(stat_boot) # Mean centered
stat_boot_2 <- stat_boot - stat_orig # Centered on the bias
Calculate p-values
(pval_boot_1 <- (sum(abs(stat_boot_1) >= abs(stat_orig)) + 1)/(R + 1))
[1] 0.03019698
(pval_boot_1 <- (sum(abs(stat_boot_2) >= abs(stat_orig)) + 1)/(R + 1))
[1] 0.03149685

Answer 1

If I understand Hall and Wilson correctly,* then the problem is in your Step 3. You should not be calculating t-statistics with respect to a null hypothesis of $\beta_1=0$ in your bootstrap samples. Instead, for the bootstrap samples you should be calculating t-statistics with respect to a null hypothesis $\beta_1= \hat \beta_1$, where $\hat \beta_1$ is your coefficient estimate from the original model.

They summarize their general recommendations for bootstrap-based tests, with $\hat \theta$ and $\hat\sigma$ being the original sample estimates of location and scale and asterisks representing corresponding estimates from bootstrap samples, in their Second Guideline:

Base the test on the bootstrap distribution of $(\hat \theta^* -\hat \theta)/\hat \sigma^*$ , not on the bootstrap distribution of $(\hat \theta^* -\hat \theta)/\hat \sigma$ or of $(\hat \theta^* -\hat \theta)$.

That's called "bootstrap pivoting" and is what you quote from Hinkley and Davison in a comment. Hall and Wilson explain in their First Guideline that evaluating bootstrap differences from the original sample estimate $(\hat\theta^*-\hat \theta)$ instead of from the original null hypothesis, $(\hat\theta^*-\theta_0)$, increases the power of the test, potentially by a lot if the original null hypothesis is far from correct.

In your example, you should do the shift $(\hat\beta_1^*-\hat\beta_1)$ before you divide by $\hat\sigma^*$ to get each of your bootstrapped t-statistics. That nicely brings the center of the distribution to 0 in the best way while respecting the bootstrap principle.

---

*reference suggested by the late Michael Chernick in this answer