# How to obtain p-values of coefficients from bootstrap regression?

From Robert Kabacoff's Quick-R I have

# Bootstrap 95% CI for regression coefficients
library(boot)
# function to obtain regression weights
bs <- function(formula, data, indices) {
d <- data[indices,] # allows boot to select sample
fit <- lm(formula, data=d)
return(coef(fit))
}
# bootstrapping with 1000 replications
results <- boot(data=mtcars, statistic=bs,
R=1000, formula=mpg~wt+disp)

# view results
results
plot(results, index=1) # intercept
plot(results, index=2) # wt
plot(results, index=3) # disp

# get 95% confidence intervals
boot.ci(results, type="bca", index=1) # intercept
boot.ci(results, type="bca", index=2) # wt
boot.ci(results, type="bca", index=3) # disp


How can I obtain the p-values $H_0:\, b_j=0$ of the bootstrap regression coefficients?

• "the p values" means what? What specific test with what null hypothesis? Commented Jan 21, 2014 at 18:38
• Correction H0:bj=0
– ECII
Commented Jan 21, 2014 at 18:42
• You already get $p<0.05$/$p>0.05$ based on whether the confidence interval does not/does include 0. Any more detail is not possible since the distribution of the parameter from the bootstrap is not parametric (and thus you can not get a probability that the value is 0). Commented Jan 21, 2014 at 19:00
• If you cannot assume a distribution how do you know that p<0.05 if the CI do not include 0? This holds true for the z or t distrubutions.
– ECII
Commented Jan 21, 2014 at 19:52
• I get that but you can only say that p<0.05, you cannot attach a specific value right?
– ECII
Commented Jan 21, 2014 at 20:14

Just another variant that is somewhat simplistic but I think deliver the message without explicitly using the library boot that may confuse some people with the syntax it uses.

We have a linear model: $y = X \beta + \epsilon$, $\quad \epsilon \sim N(0,\sigma^2)$

The following is a parametric bootstrap for that linear model, that means that we do not resample our original data but actually we generate new data from our fitted model. Additionally we assume that the bootstrapped distribution of the regression coefficient $\beta$ is symmetric and that is translation invariant. (Very roughly speaking that we can move the axis of it with affecting its properties) The idea behind is that the fluctuations in the $\beta$ 's are due to $\epsilon$ and therefore with enough samples they should provide a good approximation of the true distribution of $\beta$ 's. As before we test again $H_0 : 0 = \beta_j$ and we defined our p-values as "the probability, given a null hypothesis for the probability distribution of the data, that the outcome would be as extreme as, or more extreme than, the observed outcome" (where the observed outcomes in this case are the $\beta$ 's we got for our original model). So here goes:

# Sample Size
N           <- 2^12;
# Linear Model to Boostrap
Model2Boot  <- lm( mpg ~ wt + disp, mtcars)
# Values of the model coefficients
Betas       <- coefficients(Model2Boot)
# Number of coefficents to test against
M           <- length(Betas)
# Matrix of M columns to hold Bootstraping results
BtStrpRes   <- matrix( rep(0,M*N), ncol=M)

for (i in 1:N) {
# Simulate data N times from the model we assume be true
# and save the resulting coefficient in the i-th row of BtStrpRes
BtStrpRes[i,] <-coefficients(lm(unlist(simulate(Model2Boot)) ~wt + disp, mtcars))
}

#Get the p-values for coefficient
P_val1 <-mean( abs(BtStrpRes[,1] - mean(BtStrpRes[,1]) )> abs( Betas[1]))
P_val2 <-mean( abs(BtStrpRes[,2] - mean(BtStrpRes[,2]) )> abs( Betas[2]))
P_val3 <-mean( abs(BtStrpRes[,3] - mean(BtStrpRes[,3]) )> abs( Betas[3]))

#and some parametric bootstrap confidence intervals (2.5%, 97.5%)
ConfInt1 <- quantile(BtStrpRes[,1], c(.025, 0.975))
ConfInt2 <- quantile(BtStrpRes[,2], c(.025, 0.975))
ConfInt3 <- quantile(BtStrpRes[,3], c(.025, 0.975))


As mentioned the whole idea is that you have the bootstrapped distribution of $\beta$ 's approximates their true one. (Clearly this code is optimized for speed but for readability. :) )

• I am not familiar with R but why do rep(0,M*N) when creating the BtStrpRes matrix only to fill the N first rows in the loop?
– LucG
Commented Nov 7, 2022 at 9:51
• rep(0,M*N) creates a vector of $M*N$ of all zeros, matrix(rep(0,M*N), ncol=M) makes that vector into a matrix of $M$ columns (and therefore $N$ rows). The for loop goes over $N$ rows so it actually fills all the rows of the matrix. :) Commented Nov 7, 2022 at 10:34

The community and @BrianDiggs may correct me if I am wrong, but I believe you can get a p-value for your problem as follows. A p-value for a two sided test is defined as

$$2*\text{min}[P(X \le x|H_0),P(X \ge x|H_0)]$$

So if you order the bootstrapped coefficients by size and then determine the proportions larger and smaller zero, the minimum proportion times two should give you a p-value.

I normally use the following function in such a situation:

twosidep<-function(data){
p1<-sum(data>0)/length(data)
p2<-sum(data<0)/length(data)
p<-min(p1,p2)*2
return(p)
}

• Why ordering the coefficients by size is useful? I don't see such ordering in the code you provided. I don't understand either why |H0 appears in your formula? We don't know if H0 is true. Thank you in advance :)
– LucG
Commented Nov 7, 2022 at 10:08
• @LucG see e.g. en.wikipedia.org/wiki/P-value . The ordering is not required, but only used to describe in words what's happening. Commented Dec 1, 2022 at 17:25
• Here is a reference for why this works: stat.umn.edu/geyer/5601/examp/tests.html. In my limited empirical testing, this answer was more accurate than the accepted one, but your mileage my vary. Commented May 29 at 14:27

The bootstrap can be used to compute $p$-values, but it would need a substantial change to your code. As I am not familiar with R I can only give you a reference in which you can look up what you would need to do: chapter 4 of (Davison and Hinkley 1997).

Davison, A.C. and Hinkley, D.V. 1997. Bootstrap methods and their application. Cambridge: Cambridge University Press.