---
title: "Principal component regression (PCR)"
author: "Trevor Hastie and Robert Tibshirani"
output:
  pdf_document: default
---

**Here, I am adapting part of the lab associated with Chapter 6 of the textbook.** 

# Data Prep

We wish to predict a baseball player's `Salary` on the basis of various statistics associated with performance in the previous year.

First of all, we note that the `Salary` variable is missing for some of the players.  The `is.na()` function can be used to identify the missing observations. It returns a vector of the same length as the input vector, with a `TRUE` for any elements that are missing, and a `FALSE` for non-missing elements.
 The `sum()` function can then be used to count all of the missing elements.

```{r chunk1}
library(ISLR2)
data(Hitters)
names(Hitters)
#View(Hitters)
dim(Hitters)
sum(is.na(Hitters$Salary))
```

Hence we see that `Salary` is missing for $59$ players. The `na.omit()` function removes all of the rows that have missing values in any variable.

```{r chunk2}
Hitters <- na.omit(Hitters)
dim(Hitters)
sum(is.na(Hitters))
```

All is well now, so we `attach` the data for ease of use. 

```{r}
attach(Hitters)
```

# Principal Components Regression

Principal components regression (PCR) can be performed using the `pcr()` function, which is part of the `pls` library. We now apply PCR to the `Hitters` data, in order to predict `Salary`. 

```{r chunk42}
#install.packages("pls")
library(pls)
set.seed(1)
pcr.fit <- pcr(Salary ~ ., data = Hitters, scale = TRUE,
    validation = "CV")
```

The syntax for the `pcr()` function is similar to that for `lm()`, with a few additional
options. Setting `scale = TRUE` has the effect of *standardizing* each
predictor, prior to generating the principal
components, so that the scale on which each variable is measured will not have an effect.
 Setting `validation = "CV"` causes
`pcr()` to compute the ten-fold cross-validation error for each possible
value of $M$, the number of principal components used. The resulting fit can be examined using `summary()`.

```{r chunk43}
summary(pcr.fit)
```

The CV score is provided for each possible number of components, ranging
from $M=0$ onwards. Note that  `pcr()` reports the *root mean squared error*; in order to obtain the usual MSE, we must square this quantity. For instance, a root mean squared error of $352.8$ corresponds to an MSE of
$352.8^2=124{,}468$.

One can also plot the cross-validation scores using the
`validationplot()` function. Using `val.type = "MSEP"`
will cause the cross-validation MSE to be plotted.

```{r chunk44}
validationplot(pcr.fit, val.type = "MSEP")
```

We see that the smallest cross-validation error occurs when $M=18$ components are used. This is barely fewer than $M=19$, which amounts to simply performing least squares, because when all of the components are used in PCR no dimension reduction occurs. However, from the plot we also see that the cross-validation error is roughly the same when only one component is included in the model. This suggests that a model that uses just a small number of components might suffice.

The `summary()` function also provides the *percentage of variance explained* in the predictors and in the response using different numbers of components. Briefly, we can think of this as
the amount of information about the predictors or the
response that is captured using $M$ principal components. For example, setting $M=1$ only captures $38.31 \%$ of all the variance, or information, in the predictors. In contrast, using $M=5$ increases the value to $84.29 \%$. If we
were to use all $M=p=19$ components, this would increase to $100 \%$.

We now perform PCR on the training data and evaluate its test set performance.

```{r chunk45}
set.seed(1)
#splitting the data set
train <- sample(c(TRUE, FALSE), nrow(Hitters),
    replace = TRUE)
test <- (!train)

pcr.fit <- pcr(Salary ~ ., data = Hitters, subset = train,
    scale = TRUE, validation = "CV")
validationplot(pcr.fit, val.type = "MSEP", legendpos = t(c(12,200000)))
summary(pcr.fit)
```

Now we find that the lowest cross-validation error occurs when $M=7$ components are used.
We compute the test MSE as follows.

```{r chunk46}
x <- model.matrix(Salary ~ ., Hitters)[, -1]
#x
y <- Hitters$Salary
y.test <- y[test]

pcr.pred <- predict(pcr.fit, x[test, ], ncomp = 7)
mean((pcr.pred - y.test)^2)
```


Finally, we fit PCR on the full data set, using $M=7$, the number of components identified by cross-validation.

```{r chunk47}
pcr.fit <- pcr(y ~ x, scale = TRUE, ncomp = 7)
summary(pcr.fit)
```