Calculating a statistic via linear regression for a bootstrap procedure in R

Question

I am interested in finding the interval estimate of the mean of my response variable when the predictor variable is equal to a certain value, that is $E[Y|X=x]$. The point estimate for this statistic I am bootstrapping is $\hat{\beta_0} + \hat{\beta_1}x$.

However, I am having trouble calculating this statistic in the bootstrap replication. Here is what I have so far:

# load the libraries
library(purrr)
library(caret)

# set the seed 
set.seed(872)

# load the data
boston = tibble::as_tibble(MASS::Boston)

# store the indices for use in 400 bootstrap resamples
index_resampling = createResample(boston$medv, times = 400)

# create a bootstrap resampling function
create_boot_resampling = function(data, idx) {
  data[idx]
}

# fit the model
mod = lm(medv ~ rm, data = boston)

# calculate the bootstrap resamples
bootstrap_resample = map(index_resampling, ~create_boot_resampling(data = boston$medv, idx = .x))

# calculate the bootstrap replicates
bootstrap_replicate = unname(map_dbl(bootstrap_resample, ~predict(lm(.x ~ rm, data = boston), data.frame(rm = 4.929))))

# find the 90% confidence interval
quantile(bootstrap_replicate, probs = c(0.05, 0.950))

The medv is the response variable, and the rm is the predictor variable from the Boston dataset. I would like to find the point estimates of medv when rm is equal to 4.929, that is $\hat{\beta_0} + \hat{\beta_1} *$ 4.929. The code runs fine, but the confidence interval it outputs is incorrect. I have a feeling my error lies in the predict function used in the bootstrap replicate code, but I am unsure how to get the correct answer. Any hints or advice would be helpful!

Answer 1

As already explained in my first comment I believe your problem is that you bootstrapped medv but not rm. I tried to change your code to work appropriately by conserving as much of your code as possible. Even though I myself would not have defined boston as a tibble copy of MASS::Boston when nothing is done to boston that is tibble specific and I myself would not have defined a create_boot_resampling function that does so little.

I hope you learn most from changes in your old code, otherwise Roland has posted a good solution as the first comment to the question.

set.seed(872)

# load the data
boston = tibble::as_tibble(MASS::Boston)

# store the indices for use in 400 bootstrap resamples
index_resampling = createResample(boston$medv, times = 400)

# create a bootstrap resampling function
create_boot_resampling = function(data, idx) {
  data[idx]
}

# fit the model
mod = lm(medv ~ rm, data = boston)

# calculate the bootstrap resamples
#bootstrap_resample = map(index_resampling, ~create_boot_resampling(data = boston$medv, idx = .x))
# instead: make a list of data.frames each containing bootstrap samples of medv and rm 
bootstrap_rs <- lapply(index_resampling, function(is) boston[is, c("medv", "rm")])

# calculate the bootstrap replicates
#bootstrap_replicate = unname(map_dbl(bootstrap_resample, ~predict(lm(.x ~ rm, data = boston), data.frame(rm = 4.929))))
bootstrap_replicate <- sapply(1:length(bootstrap_rs), 
                              function(i) predict(lm(medv ~ rm, bootstrap_rs[[i]]), data.frame(rm = 4.929)))

# find the 90% confidence interval
quantile(bootstrap_replicate, probs = c(0.05, 0.950))

# always look at your data to confirm everything worked fine
hist(bootstrap_replicate, freq = FALSE)
lines(density(bootstrap_replicate), lwd = 2)
rug(bootstrap_replicate)

This gives the same result as Rolands and my comment code so we can assume this is correct. I apologize for using lapply and sapply which are standard R function when obviously you prefer purrr but I assume there will be an easy way to translate that to purrr.