Hypothesis testing with non-parametric bootstrap on beta parameter of linear model [closed]

Question

I'd like to test the value of a $\beta$ parameter estimated on a linear model. To be specific I want to test $H_0: \beta_1 = -3$ but I'm not sure if the procedure I used in R is correct.

This is the code I used:

attach(mtcars)
require(boot)
set.seed(5)

# Define the linear model
fit = lm(mpg ~ hp)
sum = summary(fit)
sum

# Define the null hypothesis for beta1
H0 = 1

Then, I calculate the test statistic with

$$T_{oss} = \frac{\hat{\beta_1} -\beta_1}{\sqrt{\hat{V}(\hat{\beta_1})}}$$

toss = (coefficients(fit)[2]-H0)/sum$coefficients[4]; toss
pval = 2*(1-pt(toss, length(hp)-2)); pval

# bootstrap procedure
e = rstandard(fit)
fitted = fit$fitted.values
toss.boot <- function(x, idx){
  y = fitted + e[idx]
  fit = lm(y~hp)
  sum = summary(fit)
  toss1 = (coefficients(fit)[2]-H0)/sum$coefficients[4]
  return (toss1)
}

out.boot = boot(mtcars, statistic = toss.boot, R=1000)

# pvalue as proportion of values more extreme than the one observed
pval = sum(out.boot$t>toss)/length(out.boot$t)
pval

I'm pretty sure there is something wrong but I would like to get your feedback and hopefully suggestions on other/better approaches.

Answer 1

Thanks for your edits. There are still some issues with your code, as the cyl variable continues to appear throughout the code, even though you replaced it with the hp variable.

I suspect the reason your code doesn't produce the intended results is because you are not simulating data according to the null hypothesis for the computation of your bootstrapped p-value. If you look at the post on http://dwoll.de/rexrepos/posts/resamplingBootALM.html, you will find an example of such simulation in the context of performing an ANOVA F-test rather than a t-test. However, the principle is similar and should look more or less like this:

# Test null hypothesis that coefficient of hp is equal to 1 
# using the original data
modelBase <- lm(I(mpg-1*hp) ~ hp, data=mtcars)
tBase  <- summary(modelBase)$coefficients["hp","t value"]
tBase  ## test statistic
pBase <- summary(modelBase)$coefficients["hp","Pr(>|t|)"]
pBase  ## p-value

fit0 <- lm(I(mpg-1*hp) ~ 1, data=mtcars)  ## fit null model
E    <- residuals(fit0)               ## residuals from null model 
Er   <- E / sqrt(1-hatvalues(fit0))   ## rescaled residuals from null model
Yhat <- fitted(fit0)                  ## fitted values from null model 

## function for performing model-based resampling 
gett <- function(dat, idx) {
        Ystar <- Yhat + Er[idx] 
        model  <- lm(Ystar ~ hp, data=dat)  
        summary(model)$coefficients["hp","t value"]
}

library(boot)
nR <- 999
bootres <- boot(mtcars, statistic=gett, R=nR)
bootres

bootres$t

hist(bootres$t)
abline(v=tBase, col="red", lty=2)

## Be careful with the computation of the bootstrap p-value, as it 
## needs to reflect what alternative hypothesis you are testing: 
## Ha: beta_hp > 1   OR
## Ha: beta_hp < 1   OR 
## Ha: beta_hp != 1
## 
## pvalue as proportion of values more extreme than the one observed
## pval <- sum(bootres$t > tBase)/length(bootres$t)
## pval

In the above, I used the fact that it is possible to test a null hypothesis like "coefficient of hp is equal to 1" by fitting the model:

lm(I(mpg-1*hp) ~ hp, data=mtcars)

and then looking at the t-test corresponding to the hp predictor variable. (See here for details: Test model coefficient (regression slope) against some value.)