What is the value of bootstrapping residuals? I have come across a question in Gelman - Regression and Other Stories that asks;

Using data of interest to you, fit a model of interest.
(a) Simulate replicated datasets and visually compare to the actual data

Using information contained within the chapter I take this as being fit a model:
y = X β + error

compute residuals of this model
r = y − X β

and bootstrap these residuals to form new fitted values of form:
y_boot = X β + r_boot

I then compare y_boot to the y value of the original data.
I have done this with the following R code below using the BostonHousing data from package mlbench;
library(mlbench)
library(plyr)
library(dplyr)
library(stringr)
library(tidyr)
library(ggplot2)
library(forcats)

data("BostonHousing")
df <- BostonHousing %>% mutate(ID = 1:nrow(.))

gen_datasets <- function(n_sims){
  
  gen_boot.res <- function(){
    fit <- lm(medv ~ crim, df)
    boot_res <- sample(fit$residuals, replace = TRUE)
    fit_boot <- predict(fit) + boot_res
    return(fit_boot)
  }
  output <- tibble()
  for(i in 1:n_sims){
    output[1:506, paste('Sim_', i ,sep = '')] = gen_boot.res()
  }
  return(output)
  }
  
df.plot <- cbind(df %>% select(medv, ID, crim), gen_datasets(20)) %>% 
  pivot_longer(cols = 4:ncol(.),
               names_to = 'sims',
               values_to = 'vals') 

ggplot(df.plot, aes( crim, medv, group = 1)) + 
  geom_point() + 
  stat_smooth(aes(x = crim, y = vals), method = "lm",se = TRUE, formula = y ~ x) +
  facet_wrap(~fct_relevel(sims, str_c('Sim_' , c(1:20))), scales = 'free_y')

ggplot(df.plot, aes( crim, medv, group = sims)) + 
  geom_point() + 
  stat_smooth(aes(x = crim, y = vals),
              method = "lm",se = FALSE, 
              formula = y ~ x, geom = 'line',
              alpha = .5, colour = 'red'
              )





This reply is primarily comments aimed at creating a reliable workflow for such an analysis.
The suggested procedure is this:

*

*Prepare the data by creating a complete dataset with just the needed variables and no missing values.  This reduces the chances of errors when the software (perhaps silently) removes incomplete cases.


*Fit the model.  Encapsulate the fitting process in a procedure or function for reuse.


*Perform the bootstrapping in a loop by

*

*creating a bootstrap sample

*fitting the sample

*storing the sample and any essential results of the fit (coefficients, etc.)



*Examine the results both visually and with statistical summaries.

*

*Include the original data and fit for comparison.

*Consider a "data lineup" for comparing a small set of bootstrap samples to the data.



*Based on this examination return to step (2) or step (1) and repeat as needed.

The "data lineup" is a way to assess the "significance" of any result, visually and quickly.  If you present, say, 20 versions of the data and they all differ from the data in some clear way, you have visual evidence that the data differ from the procedure you used to generate those 20 versions.  This is the statistical analog of the familiar police lineup used to identify a crime suspect by viewing that suspect within a group of comparable people.
Some tips:

*

*Do not rely on the graphics software to perform your analysis.  Use it to display the results of your analysis.


*Encapsulate the basic operations of fitting a model, extracting its results, and plotting those results.  This helps ensure consistent, reliable processing of the data and the bootstrap samples.  See the analyze, predict., and residuals. functions in the R code below.
The R code below illustrates the foregoing points.  In particular, notice the ease of bootstrapping: it's just a loop over two lines of code,
    X$medv <- y.hat + sample(residuals.(fit), replace=TRUE)
    y <- predict.(analyze(X), newdata=W)

The first line performs the bootstrap sampling while the second redoes the fit and extracts the information of interest (in this case, y-values for the graph of the fitted values).
This code produced the following plots at step (4).

Notice that this model is not a simple linear relationship between medv and crim: I have fit a more suitable logarithmic one.  The opportunity to fit exactly the right model is the main reason not to rely on the graphics program to do the fitting for you: it often will not be flexible enough.

This is the same plot, displayed with one panel per iteration.

This is the same as the previous plot, displayed now on a logarithmic axis. Much more can be seen now.  Some outliers are evident in the upper right corner of the data (Panel 20) and there's some evidence of positive skewness in the responses -- perhaps even bimodality.  These both suggest avenues for further analysis.
Finally, R users might appreciate the effective use of the Grammar of Graphics approach to visualization (as implemented in ggplot2): the graphical display is described once (with almost no customization of scales) and then redrawn multiple times with added variations.  This economy of expression helps assure consistency in the output.
I hope this small example illustrates the merits of the approach suggested by the original exercise.
library(ggplot2)
#
# Access and prepare the data.
#
library(mlbench)
data(BostonHousing)
df <- subset(BostonHousing, select=c(medv, crim))
df <- subset(df, complete.cases(df)) # Essential for handling missing data!
#
# Perform the original analysis.
#
analyze <- function(df) lm(medv ~ log(crim), df)
predict. <- function(f, ...) predict(f,  ...)
residuals. <- function(f, ...) residuals(f, ...)
fit <- analyze(df)
#
# Prepare for bootstrapping.
#
n.boot <- 20
set.seed(17)
y.hat <- predict.(fit)
n <- length(residuals.(fit))
W <- data.frame(crim = exp(seq(log(min(df$crim)), log(max(df$crim)), length.out=101)))
#
# Obtain bootstrap samples.
#
B.boot <- lapply(1:n.boot, function(i) {
  X <- cbind(data.frame(Iteration = i), df)
  if(i < n.boot) {
    # This is one draw from the bootstrap distribution
    X$medv <- y.hat + sample(residuals.(fit), replace=TRUE)
    y <- predict.(analyze(X), newdata=W)
    status <- "Bootstrap"
  } else {
    # The last time through, just store the data and its analysis.
    y <- predict(fit, newdata=W)
    status <- "Data"
  }
  list(Data = cbind(X, data.frame(Status = status)), 
       Fit = cbind(data.frame(Iteration = i, Status = status, medv = y), W))
})
#
# Prepare for post-processing by assembling the results into data frames.
#
Points <- do.call(rbind, lapply(B.boot, function(l) l[["Data"]]))
Fits <- do.call(rbind, lapply(B.boot, function(l) l[["Fit"]]))
#
# Display versions of the plot.
#
G <- ggplot(Points, aes(crim, medv, group=Iteration)) + 
  geom_point(aes(fill=Status), shape=21, alpha=1/2) + 
  geom_line(aes(color=Status), data=Fits, size=0.95) + 
  ggtitle("All Data, Original and Bootstrapped")

print(G)
print(H <- G + facet_wrap(~ Iteration))
#
# More displays as suggested by review.
#
print (H + scale_x_log10() + 
         theme(axis.text.x = element_text(angle=90)) +
         ggtitle("Results on a Logarithmic X Scale"))

