---
title: "Cars"
author: "Trevor Hastie and Robert Tibshirani"
output:
  pdf_document: default
---

Here, I am adapting the lab associated with Chapter 5 of the textbook. 

```{r}
library(ISLR2)
library(boot)
```

### Estimating the Accuracy of a Linear Regression Model

The bootstrap approach can be used  to assess the
variability of the coefficient estimates and predictions from a statistical learning method. Here we use the bootstrap approach in order to assess the variability of
the estimates for $\beta_0$ and $\beta_1$, the intercept and slope terms for the *simple* linear regression model
that uses  `horsepower` to predict `mpg` in the `Auto` data set. We will compare the estimates obtained using the bootstrap to those obtained using the formulas
for ${\rm SE}(\hat{\beta}_0)$ and ${\rm SE}(\hat{\beta}_1)$ described in Section 3.1.2 (*and the slides from class*).

*Let's make some plots of the data to begin with.*

```{r}
Auto
attach(Auto)
names(Auto)
#start with the scatterplot 
plot(horsepower, mpg, 
     main="Dependence of efficiency on engine power",
     pch=20, col="blue")
```

It looks suspiciously non-linear. So, let's add the least-squares line. 

```{r}
plot(horsepower, mpg, 
     main="Dependence of efficiency on engine power",
     pch=20, col="blue")
reg=lm(mpg ~ horsepower)
summary(reg)
abline(reg, col="red", lwd=2)
```

Now, what about the residuals? We want to see if the residuals have an association with the explanatory, i.e., engine power. 
```{r}
res=summary(reg)$residuals
plot(horsepower, res, 
     main="Residuals", 
     xlab="Horsepower", ylab="Residuals",
     pch=20, col="blue")
abline(0,0, col="red", lwd=2)
```

We first create a simple function, `boot.fn()`, which takes in the `Auto` data set as well as a set of indices for the observations, and
returns the intercept and slope estimates for the linear regression model. We then apply this function
to the full set of $n=392$ observations in order to compute the estimates of $\beta_0$ and $\beta_1$ on the entire data set using the usual linear regression coefficient estimate
formulas from Chapter 3. Note that we do not need the `{` and `}` at the beginning and end of the function because it is only one line long.

```{r chunk15}
boot.fn <- function(data, index)
  coef(lm(mpg ~ horsepower, data = data, subset = index))
boot.fn(Auto, 1:392)
```

The `boot.fn()` function can also be used in order to create
bootstrap estimates for the intercept and slope terms by randomly sampling from among the observations with replacement. Here we give two examples.

```{r chunk16}
set.seed(1)
boot.fn(Auto, sample(392, 392, replace = T))
boot.fn(Auto, sample(392, 392, replace = T))
```

Next, we use the `boot()` function to compute the standard errors of 1,000 bootstrap estimates for the intercept and slope terms.

```{r chunk17}
boot(Auto, boot.fn, R=1000)
```

This indicates that the bootstrap estimate for ${\rm SE}(\hat{\beta}_0)$ is $0.84$, and that the bootstrap estimate for ${\rm SE}(\hat{\beta}_1)$ is $0.0073$.
As discussed in Section 3.1.2, standard formulas can be used to compute the standard errors for the regression coefficients in a linear model. These can be obtained using the  `summary()` function.

```{r chunk18}
summary(lm(mpg ~ horsepower, data = Auto))$coef
```


The standard error estimates for $\hat{\beta}_0$ and
$\hat{\beta}_1$ obtained using the formulas from
Section 3.1.2 are $0.717$ for the intercept and $0.0064$
for the slope. Interestingly, these are somewhat different from the
estimates obtained using the bootstrap.  Does this indicate a problem
with the bootstrap? In fact, it suggests the opposite.  Recall that
the standard formulas for the standard errors rely on certain assumptions. For example, they depend
on the unknown parameter $\sigma^2$, the noise variance. We then estimate $\sigma^2$
using the RSS. **Now, although the formulas for the standard errors do not rely on the linear model being correct, the estimate for $\sigma^2$ does.**
We **earlier** 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$.
Secondly, the standard formulas assume (somewhat unrealistically) that the $x_i$ are fixed, and all the variability comes from the variation in the errors $\epsilon_i$.
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.

Below we compute the bootstrap standard error estimates and the standard linear regression estimates that result from fitting the quadratic model to the data. Since this model provides a good fit to the data (Figure 3.8), there is now a better correspondence between the bootstrap estimates and the standard estimates of ${\rm SE}(\hat{\beta}_0)$, ${\rm SE}(\hat{\beta}_1)$ and ${\rm SE}(\hat{\beta}_2)$.

```{r chunk19}
boot.fn <- function(data, index)
  coef(
      lm(mpg ~ horsepower + I(horsepower^2), 
        data = data, subset = index)
    )
set.seed(1)
boot(Auto, boot.fn, 1000)
q.reg=lm(mpg ~ horsepower + I(horsepower^2), data = Auto)
betas=q.reg$coef
betas
```
How about a picture?
```{r}
plot(horsepower, mpg, 
     main="Dependence of efficiency on engine power",
     pch=20, col="blue")
b.0=betas[[1]]
b.1=betas[[2]]
b.2=betas[[3]]
curve(b.0+b.1*x+b.2*x^2, col="red", lwd=2, add=TRUE)
```
Significance?
```{r}
summary(q.reg)
```