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. :) )