# Bootstrap estimate for standard error in linear regression

I came across the bootstrap concept in the book Introduction to statistical learning, wherein the standard error for the linear regression coefficients is estimated both by the formula and the bootstrap process, and then following this paragraph (page 196):

Now although the formula for the standard errors do not rely on the linear model being correct, the estimate for $\sigma^2$ does. We see [...] that there is a non-linear relationship in the data, and so the residuals from a linear fit will be inflated, and so will $\hat\sigma^2$. [...] The bootstrap approach does not rely on any of these assumptions, and so it is likely giving a more accurate estimate of the standard errors of $\hat\beta_0$ and $\hat\beta_1$ than is the summary() function. [Emphasis added.]

So my question is why are the assumptions reducing the reliability or accuracy of the estimates? Because more is the residual error more is your standard error of estimates (indirectly due to inflated $\hat\sigma^2$) this assumption doesn't seem wrong to me.

• The bootstrap goes back to Efron's articles and his monograph in 1979. In regression there are mainly two ways to bootstrap. One involves resampling the model residuals and the other is to resample the n+1 dimensional vector of the response variable (usually denoted as y) and the n predictor variables. In the case when the correct model is nonlinear bootstrapping the vector is less dependent on linearity and hence generally leads to a better fit and more accurate standard errors. – Michael Chernick Apr 5 '17 at 17:22

On the surface your assessment isn't wrong: a worse fit gives larger residuals, which gives a larger standard error estimate. But this is not the whole story!

The standard error of the estimate $\hat \beta$ isn't to assess model fit. If all you care about is how well the line fits your data, you don't need to bother with the standard error of $\hat \beta$ at all. Just use $\hat\sigma^2$ directly, which is the MSE of your model. We want to use $\text{s.e}(\hat \beta)$ for inference. We use it to say something about how well we have estimated $\beta$. We want to for eg construct confidence intervals around our point estimate.

If we are to have any sort of confidence that these intervals relate correctly to the population $\beta$, $\text{s.e.}(\hat\beta)$ should be correct. It should be the standard deviation of whatever distribution $\hat \beta$ has. If the model is specified correctly, we can use the standard formulas to estimate $\text{s.e.}(\hat\beta)$ directly, and we know exactly which properties $\hat\beta$ has. These properties and formulas — and hence our inferences about $\hat\beta$ — are directly derived from the model assumptions.

If the model isn't correctly specified, all our beautiful and clean theory goes out the window. Elsewhere in the ISLR book you can read that an approximate 95% confidence interval for $\hat \beta$ is $$[\hat\beta - 2\cdot\text{s.e.}(\hat\beta), \hat\beta + 2\cdot\text{s.e.}(\hat\beta)].$$ This is true if you have a good standard error estimate, but if the estimate is too small/large, the confidence interval is too tight/wide.

Below is an illustration in code and figures. The examples in ISLR use some data that look quadratic. I will simulate some data so that we know what the truth is. My true model is $y = 4 + 5x -3x^2 + \epsilon$, where $\epsilon \sim N(0,1)$. This is classic linear regression, the big assumption is that errors should be iid from a normal distribution with mean zero. A quadratic model is the perfect fit. A linear model is a poor fit.

library(boot)
library(plyr)

set.seed(22042017) # for reproducibility
generate_data <- function(nsamples=100) {
x <- rnorm(nsamples,mean=.75, sd=.5)
y <-  4 + 5*x -3*x^2 + rnorm(length(x), mean=0, sd = 1)
data.frame(x,y)
}

dat <- generate_data()
plot(dat, pch=20, col="grey")
curve(4 + 5*x -3*x^2, add=T, lwd=1.5)
abline(lm(y~x, data=dat), col="red", lwd=1.5)


The figure shows some simulated data in grey, the true model in black, and a fitted linear regression in red. Now clearly the residuals aren't iid. They are below zero to the right and to the left, and above zero in the middle. There's no real reason to think $\hat\sigma^2 = \frac{1}{n}\sum \hat\epsilon_1$ makes sense here, why take the variance of some stuff that's not from the same distribution?

Let's estimate $\text{s.e}(\hat \beta)$ first from the linear model and then by bootstrap:

# estimates from assumptions
summary(lm(y~x, data=dat))
#>
#> Call:
#> lm(formula = y ~ x, data = dat)
#>
#> Residuals:
#>     Min      1Q  Median      3Q     Max
#> -4.8490 -0.7340  0.1801  0.9807  3.5730
#>
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)
#> (Intercept)   4.8882     0.2466  19.826   <2e-16 ***
#> x             0.5879     0.2758   2.132   0.0355 *
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 1.455 on 98 degrees of freedom
#> Multiple R-squared:  0.04432,    Adjusted R-squared:  0.03456
#> F-statistic: 4.544 on 1 and 98 DF,  p-value: 0.03553

# bootstrap estimates, no assumptions
bootfun <- function(data, index) {
coef(lm(y~x, data=dat, subset = index))
}
boot(dat, bootfun, 10000)
#>
#> ORDINARY NONPARAMETRIC BOOTSTRAP
#>
#>
#> Call:
#> boot(data = dat, statistic = bootfun, R = 10000)
#>
#>
#> Bootstrap Statistics :
#>      original      bias    std. error
#> t1* 4.8881549 0.006852531   0.3099025
#> t2* 0.5878975 0.012233797   0.4304412


We see that the bootstrap estimates, which aren't derived from the linear regression assumptions are larger. It's hard to say that this is necessarily better without knowing the truth. But we know the true model. Since I'm just making up data, I can draw samples directly from the distribution of $\hat \beta$. Below I do this and estimate the true $\text{s.e}(\hat \beta)$ directly:

# The "real" s.e. of beta estimates
sampling <- raply(10000, function() {
dat <- generate_data()
coef(lm(y~x, data=dat))
})
aaply(sampling, 2, sd)
#> (Intercept)           x
#>   0.4079006   0.5069441


We see that even the bootstrap underestimated $\text{s.e}(\hat \beta)$, but that it does a better job of it than the good old formulas.

If we look at the residuals real quick we see that they are skewed, which is another reason not to put too much trust in $\hat \sigma^2$. It's a nice summary statistic, but it's better for symmetrical distributions.

# The residuals are clearly not normal
hist(resid(lm(y~x, data=dat)), prob=T, nclass=20)
sdev <- sd(resid(lm(y~x, data=dat)))

• I am not questioning the correctness of your answer, but I would just like to point out that when talking of i.i.d.ness, there might a difference between conditional vs. unconditional distributions. For example, your examples illustrates a dependence between residuals and the regressor but perhaps not the i.i.d.ness (or lack thereof) of residuals alone. Your statement Now clearly the residuals aren't iid is valid for the conditional distribution of residuals given the regressor $x$; but is it valid for the unconditional distribution? Or is raising such a question problematic in itself? – Richard Hardy Oct 22 '18 at 8:45