Baseball.Rmd

---
title: "Baseball"
output: html_document
---

Now we will aply the best subset selection approach to the Hitters dataset.  Recall, this involves evaluating all possible combinations of parameters for model selection.  

```{r}
library(ISLR)
fix(Hitters)
names(Hitters)
nrow(Hitters)
```
Check the dimensions of the data.  Next, find the number of players that are missing salary data:
```{r}
dim(Hitters)
sum(is.na(Hitters$Salary))
```
We see there are 322 players (n) with 20 data features (p).  Of the 322 players, 59 are missing salary data, so we will omit them from our analysis using the "na.omit" function.  Also, we recheck the dimension of the dataframe after applying na.omit.
```{r}
Hitters=na.omit(Hitters)
dim(Hitters)
```

In order to perform the best subset selection using the regsubsets() function.
```{r}
library(leaps)
regfit.full =regsubsets(Salary~.,data=Hitters)
summary(regfit.full)
```
Above is a summary of the top selected models, from p=1 to p=8.  Asterisks indicate the inclusion of a predictor in a given model (each row is a model).

Note that, by default, the regsubsets() function considers from p=1 up to p=8.  However, we can set the "nvmax" parameter in order to increase the upper limit for p.  Below, we set it to 19.
```{r}
regfit.full=regsubsets(Salary~.,data=Hitters, nvmax=19)
reg.summary = summary(regfit.full)
```
We can examine the residual sum of squares and the R-squared from the summary:
```{r}
names(reg.summary)
```
```{r}
reg.summary$rss
```
Note that the residual sum of squares in monotonically decreasing as we increase to 19 predictors.
```{r}
reg.summary$rsq
```

Note that the R-squared value is monotonically increasing as we increase to 19 predictors.

We can plot, RSS, the adjusted R-squared, C-p, and BIC for all 19 models to help us pick the optimal model.

```{r}
par(mfrow=c(2,2))
plot(reg.summary$rss, xlab="Number of Variables", ylab="RSS", type='l')
plot(reg.summary$adjr2, xlab="Number of Variables", ylab="Adjusted R-squared", type='l')
plot(reg.summary$cp, xlab="Number of Variables", ylab="cp", type='l')
plot(reg.summary$bic, xlab="Number of Variables", ylab="BIC", type='l')
```
We can found the extremum values indicating the top performing model according to each metric.
```{r}
which.min(reg.summary$rss)
which.max(reg.summary$adjr2)
which.min(reg.summary$cp)
which.min(reg.summary$bic)
```
We can highlight these values on the graphs by doing the following:
```{r}
par(mfrow=c(2,2))

plot(reg.summary$rss, xlab="Number of Variables", ylab="RSS", type='l')


plot(reg.summary$adjr2, xlab="Number of Variables", ylab="Adjusted R-squared", type='l')
points(11, reg.summary$adjr2[11], col="red", cex=2, pch=20)

plot(reg.summary$cp, xlab="Number of Variables", ylab="cp", type='l')
points(10, reg.summary$cp[10], col="red", cex=2, pch=20)

plot(reg.summary$bic, xlab="Number of Variables", ylab="BIC", type='l')
points(6, reg.summary$bic[6], col="red", cex=2, pch=20)
```

We can also use the plot() command that is built into the regsubsets() function.
```{r}
plot(regfit.full, scale="r2")
plot(regfit.full, scale="adjr2")
plot(regfit.full, scale="Cp")
plot(regfit.full, scale="bic")
```
We can see from the above that several models share a BIC around -150, but the lowest BIC belongs to a six predictor model.  We can find out which predictors those are by calling the coef() functin on the regfit.full model.  Note we also have to specify the number of predictors (6).
```{r}
coef(regfit.full, 6)
```


We can also use forward or backward stepwise selection in order to generate our model.  We still use the regsubsets() function:
```{r}
regfit.fwd=regsubsets(Salary~., data=Hitters, nvmax=19, method="forward")
summary(regfit.fwd)

regfit.bwd=regsubsets(Salary~., data=Hitters, nvmax=19, method="backward")
summary(regfit.bwd)
```
Interestingly, the models found by the forward stepwise, backward stepwise, and best subset are different for 7 predictors.
```{r}
coef(regfit.full, 7)
coef(regfit.fwd, 7)
coef(regfit.bwd, 7)
```