```{r, message=FALSE, warning=FALSE, echo=FALSE}
require(knitr)
opts_chunk$set(eval=FALSE)
```

This lab on Linear Regression in R comes from p. 109-119 of "Introduction to Statistical Learning with Applications in R" by Gareth James, Daniela Witten, Trevor Hastie and Robert Tibshirani. It was re-implemented in Fall 2016 in `tidyverse` format by Amelia McNamara and R. Jordan Crouser at Smith College.

#  3.6.1 Libraries

The `library()` function is used to load libraries, or groups of functions and data sets that are not included in the base R distribution. Basic functions that perform least squares linear regression and other simple analyses come standard with the base distribution, but more exotic functions require additional libraries. Here we load the `MASS` package, which is a very large collection of data sets and functions. We also load the `ISLR` package, which includes the data sets associated with this book.

```{r}
library(MASS)
library(ISLR)
```

# 3.6.2 Simple Linear Regression

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

We can fit linear models using the `lm()` function. For example, we might want to fit a model with `medv` as the response and `lstat` as the predictor, which has the syntax `medv~lstat`:

```{r}
model_SL = lm(medv~lstat)
```

Unfortunately, the above line will break because `R` doesn't know where to look for the data. We'll need to tell it explicitly that `data = Boston` within this function call. Now let's fit a simple linear model (we'll call it `model_SL`):

```{r}
model_SL = lm(medv~lstat, data=Boston)
```

If we type `model_SL`, some basic information about the model is output. For more detailed information, we use `summary(model_SL)`

```{r}
model_SL
summary(model_SL)
```

We can use the `names()` function in order to find out what other pieces of information are stored in `model_SL`. Although we can extract these quantities by name—e.g. `model_SL$coefficients` — it is safer to use extractor functions like `coef()` to access them.

```{r}
names(model_SL)
coef(model_SL)
```

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

```{r}
confint(model_SL)
```

The `predict()` function can be used to produce both confidence intervals and prediction intervals for the prediction of `medv` for a given value of `lstat`. First, we'll make a data frame with some new values for `lstat`

```{r}
new_values = data.frame(lstat=c(5,10,15))
```

Now, we'll call the `predict()` function to see what our model predicts for the corresponding `medv` value, as well as the `confidence` or `prediction` intervals:

```{r}
predict(model_SL,new_values, interval="confidence")
predict(model_SL,new_values, interval="prediction")
```

We will now plot `medv` and `lstat` along with the least squares regression line using `ggplot()` and `geom_smooth()`:

```{r}
library(ggplot2)
ggplot(Boston, aes(x=lstat, y=medv)) +
  geom_smooth(method = "lm") +
  geom_point()
```

Below we experiment with some additional settings for plotting lines and points. The `lwd=3` command causes the width of the regression line to be increased by a factor of 3. We can also use the `shape` and `size` options to create different plotting symbols.

```{r}
ggplot(Boston, aes(x=lstat, y=medv)) +
  geom_smooth(method = "lm", lwd=3, se = FALSE) +
  geom_point(shape = "+", size = 5)
```

We can use the `ggfortify` package to produce diagnostic plots. These plots are automatically produced by applying the `autoplot()` function directly to the model.

```{r}
library(ggfortify)
autoplot(model_SL)
```

Alternatively, 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. To draw these plots side by side, we can use the `grid.arrange()` function from the `gridExtra` library. We'll also use the `labs()` function to add labels to our plots.

```{r}
library(gridExtra)

plot1 = ggplot() +
  geom_point(aes(predict(model_SL), residuals(model_SL))) +
  labs(x="Predicted Value", y="Residual")

plot2 = ggplot() +
  geom_point(aes(predict(model_SL), rstudent(model_SL))) +
  labs(x="Predicted Value", y="Studentized Residual")

grid.arrange(plot1, plot2, ncol=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}
hatvalues(model_SL)
```

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.

```{r}
which.max(hatvalues(model_SL))
```

#  3.6.3 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}
model_ML = lm(medv~lstat+age, data=Boston)
summary(model_ML)
```

The Boston data set contains 13 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}
model_ML = lm(medv~., data=Boston)
summary(model_ML)
```

The `vif()` function, part of the `car` package, can be used to compute variance inflation factors. Most VIFs are low to moderate for this data. The car package is not part of the base R installation so it must be downloaded the first time you use it via the `install.packages` option in R.

```{r}
library(car)
vif(model_ML)
```

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}
model_ML_no_age = lm(medv~.-age, data=Boston)
summary(model_ML_no_age)
```

Alternatively, we can use the `update()` function to return an updated version of our previous `model_ML`.

```{r}
model_ML_no_age = update(model_ML, ~.-age)
```

# 3.6.4 Interaction Terms

It is easy to include interaction terms in a linear model using the `lm()` function. The syntax `lstat:black` tells R to include an interaction term between `lstat` and `black`. The syntax `lstat*age` simultaneously includes `lstat`, `age`, and the interaction term `lstat×age` as predictors; it is a shorthand for `lstat+age+lstat:age`.

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

# 3.6.5 Non-linear Transformations of the Predictors

The `lm()` function can also accommodate non-linear transformations of the predictors. For instance, given a predictor `X`, we can create a predictor `X2` using `I(X^2)`. The function `I()` is needed since the ^ has a special meaning in a formula; wrapping as we do allows the standard usage in R, which is to raise `X` to the power 2. We now perform a regression of `medv` onto `lstat` and `lstat2`.

```{r}
model_ML_quadratic=lm(medv~lstat+I(lstat^2), data=Boston)
summary(model_ML_quadratic)
```

The near-zero p-value associated with the quadratic term suggests that it leads to an improved model. We use the `anova()` function to further quantify the extent to which the quadratic fit is superior to the linear fit.

```{r}
anova(model_SL, model_ML_quadratic)
```

Here Model 1 represents the linear submodel containing only one predictor, `lstat`, while Model 2 corresponds to the larger quadraticmodel that has two predictors, `lstat` and `lstat2`. The `anova()` function performs a hypothesis test comparing the two models. The null hypothesis is that the two models fit the data equally well, and the alternative hypothesis is that the full model is superior. 

The F-statistic is 135, and the associated p-value is virtually zero. This provides very clear evidence that the model containing the predictors `lstat` and `lstat2` is far superior to the model that only contains the predictor `lstat`. This is not surprising, since earlier we saw evidence for non-linearity in the relationship between `medv` and `lstat`. 

If we type:

```{r}
autoplot(model_ML_quadratic)
```

then we see that when the `lstat2` term is included in the model, there is little discernible pattern in the residuals.

In order to create a cubic fit, we can include a predictor of the form `I(X^3)`. However, this approach can start to get cumbersome for higher order polynomials. A better approach involves using the `poly()` function to create the polynomial within `lm()`. For example, the following command produces a fifth-order polynomial fit:

```{r}
model_ML_5th_order_poly = lm(medv~poly(lstat, 5, raw=TRUE), data=Boston)
summary(model_ML_5th_order_poly)
```

This suggests that including additional polynomial terms, up to fifth order, leads to an improvement in the model fit! However, further investigation of the data reveals that no polynomial terms beyond fifth order have significant p-values in a regression fit.

Of course, we are in no way restricted to using polynomial transformations of the predictors. Here we try a log transformation.

```{r}
summary(lm(medv~log(rm), data=Boston))
```

# 3.6.6 Qualitative Predictors

We will now examine the `Carseats` data, which is part of the `ISLR` library. We will attempt to predict `Sales` (child car seat sales) in 400 locations based on a number of predictors.

```{r}
names(Carseats)
```

The `Carseats` data includes qualitative predictors such as `Shelveloc`, an indicator of the quality of the shelving location—that is, the space within a store in which the car seat is displayed—at each location. The predictor `Shelveloc` takes on three possible values, `Bad`, `Medium`, and `Good`.

Given a qualitative variable such as `Shelveloc`, R generates dummy variables automatically. Below we fit a multiple regression model that includes some interaction terms.

```{r}
lm.fit=lm(Sales~.+Income:Advertising+Price:Age, data=Carseats)
summary(lm.fit)
```

The `contrasts()` function returns the coding that R uses for the dummy variables. Use `?contrasts` to learn about other contrasts, and how to set them.

```{r}
contrasts(Carseats$ShelveLoc)
```

# Getting Credit
To get credit for this lab, please reply to the prompt posted to the [#lab2](https://sds293.slack.com/messages/C76KRHW14) Slack channel.