---
title: "Multiple linear regression"
author: "Trevor Hastie and Robert Tibshirani"
output:
  pdf_document: default
---

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

We load the `ISLR2` package, which includes the data sets associated with this book.

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

The `ISLR2` library contains the `Boston`  data set, which
records `medv` (median house value) for $506$ census tracts in Boston. We will seek to predict `medv` using $12$ predictors such as `rmvar` (average number of  rooms per house), `age` (proportion of owner-occupied units built prior to 1940) and `lstat` (percent of households with low socioeconomic status).

```{r chunk2}
head(Boston)
```

We should also investigate the variables therein and `attach` so that we don't have to use the `$` notation.

```{r}
names(Boston)
attach(Boston)
```

To find out more about the data set, we can type `?Boston`.

We will start by fitting a simple  linear regression model, with `medv` as the response and `lstat`  as the predictor. Here is the scatterplot.

```{r}
plot(lstat, medv, 
     main="Median value vs. lower status",
     xlab="Lower status in the population (%)",
     ylab="Median home value (in Ks)",
     pch=20, col="blue")
```

For the implementation of the linear model, the basic syntax is `lm(y ~ x, data)`, where `y` is the response, `x` is the predictor, and `data` is the data set in which these two variables are kept.


```{r chunk3, error=TRUE}
lm.fit <- lm(medv ~ lstat)
summary(lm.fit)

plot(lstat, medv, 
     main="Median value vs. lower status",
     xlab="Lower status in the population (%)",
     ylab="Median home value (in Ks)",
     pch=20, col="blue")
abline(lm.fit, col="red", lwd=2)
```
We can use the `names()` function in order to find out what other pieces of information  are stored in `lm.fit`.
Although we can extract these quantities by name---e.g. `lm.fit$coefficients`---it is safer to use the extractor functions like `coef()` to access them.

```{r chunk6}
names(lm.fit)
coef(lm.fit)
```

In order to obtain a confidence interval for the coefficient estimates, we can use the `confint()` command.

```{r chunk7}
confint(lm.fit)
```

The `predict()` function can be used to produce confidence intervals and prediction intervals for the prediction of `medv` for a given value of `lstat`.

```{r chunk8}
predict(lm.fit, data.frame(lstat = (c(5, 10, 15))),
    interval = "confidence")
predict(lm.fit, data.frame(lstat = (c(5, 10, 15))),
    interval = "prediction")
```

How about the full array of confidence and prediction intervals?

```{r}
newdata<-data.frame(lstat=seq(min(lstat),max(lstat),0.1))

conf<-predict(lm.fit,newdata,interval="confidence")
#conf
pred<-predict(lm.fit,newdata,interval="prediction")
#pred
all.ints<-cbind(conf,pred[,-1])
#all.ints

plot(lstat, medv, 
     main="Median value vs. lower status",
     xlab="Lower status in the population (%)",
     ylab="Median home value (in Ks)",
     pch=20, col="blue")
matplot(newdata, all.ints, type="l", lty=1, lwd=2, col=c("red", "purple", "purple","green","green"), add=TRUE)
```

We can compute the residuals from a linear regression
fit using the `residuals()` function. The function
`rstudent()` will return the studentized residuals, and we
can use this function to plot the residuals against the fitted values.

```{r chunk12}
plot(predict(lm.fit), residuals(lm.fit),
     main="Do residuals depend on the explanatory?",
     xlab="Explanatory",
     ylab="Residuals",
     pch=20, col="blue")
abline(0,0, col="red", lwd=2)
plot(predict(lm.fit), rstudent(lm.fit),
     main="Do residuals depend on the explanatory?",
     xlab="Explanatory",
     ylab="Residuals in t-units",
     pch=20, col="blue")
abline(0,0, col="red", lwd=2)
```

On the basis of the residual plots, there is some evidence of non-linearity.

Leverage statistics can be computed for any number of predictors using the `hatvalues()` function.

```{r chunk13}
plot(hatvalues(lm.fit), 
     pch=20, col="blue")
which.max(hatvalues(lm.fit))
```

The `which.max()` function identifies the index of the largest element of a vector. In this case, it tells us which observation has the largest leverage statistic.

## Multiple Linear Regression

In order to fit a multiple linear regression model using least squares, we again use the `lm()` function. The syntax 
`lm(y ~ x1 + x2 + x3)` is used to fit a model with three predictors, `x1`, `x2`, and `x3`.
The `summary()` function now outputs the regression coefficients for all the predictors.

```{r chunk14}
lm.fit <- lm(medv ~ lstat + age, data = Boston)
summary(lm.fit)
```

The `Boston` data set contains 12 variables, and so it would be cumbersome to have to type all of these in order to perform a regression using all of the predictors.
Instead, we can use the following short-hand:

```{r chunk15}
lm.fit <- lm(medv ~ ., data = Boston)
summary(lm.fit)
```
What if we would like to perform a regression using all of the variables but one?  For example, in the above regression output,  `age` has a high $p$-value. So we may wish to run a regression excluding this predictor. The following syntax results in a regression using all predictors except `age`.

```{r chunk17}
lm.fit1 <- lm(medv ~ . - age, data = Boston)
summary(lm.fit1)
```

Alternatively, the `update()` function can be used.

```{r chunk18}
lm.fit1 <- update(lm.fit, ~ . - age)
```

*Take away the next explanatory with a sizeable $p-$value, i.e., `indus`. *

```{r chunk17a}
lm.fit2 <- lm(medv ~ . - age - indus, data = Boston)
summary(lm.fit2)
```

*The above approach is called* **backward selection**. *How would you design* **forward selection?**

## Interaction Terms

It is easy to include interaction terms in a linear model using the `lm()` function. The syntax `lstat:age` tells `R` to include an interaction term between `lstat` and `age`.

```{r chunk19z}
summary(lm(medv ~ lstat:age, data = Boston))
```
```{r chunk19y}
summary(lm(medv ~ lstat + age + lstat:age, data = Boston))
```
```{r chunk19x}
summary(lm(medv ~ lstat + lstat:age, data = Boston))
```

The syntax `lstat * age` simultaneously includes `lstat`, `age`, and the interaction term `lstat`$\times$`age` as predictors; it is a shorthand for `lstat + age + lstat:age`.
We can also pass in transformed versions of the predictors.

```{r chunk19}
summary(lm(medv ~ lstat * age, data = Boston))
```