```{r setup, echo = FALSE}
knitr::opts_chunk$set(error = TRUE)
```

# Linear Classification

This lab is a minor modification of the lab from `An Introduction to Statistical Learning with Applications in R', Chapter 4.


## The Stock Market Data

We will begin by examining some numerical and graphical summaries of
the `Smarket` data, which is part of the `ISLR2` library. This
data set consists of percentage returns for the S\&P 500 stock index
over $1,250$ days, from the beginning of 2001 until the end of
2005. For each date, we have recorded the percentage returns for each
of the five previous trading days, `lagone` through `lagfive`. We
have also recorded `volume` (the number of shares traded on the previous day, in billions), `Today` (the percentage return on the date in question)  and `direction` (whether the market was `Up` or `Down` on this date). Our goal is to predict `direction` (a qualitative response) using the other features.

```{r chunk1}
library(ISLR2)
names(Smarket)
dim(Smarket)
summary(Smarket)
pairs(Smarket)
```

The `cor()` function produces a matrix that contains all of the pairwise correlations among the predictors in a data set. The first command below gives an error message because the `direction` variable is qualitative.

```{r chunk2, error=TRUE}
cor(Smarket)
cor(Smarket[, -9])
```

As one would expect, the correlations between the lag variables and today's returns are close to zero. In other words, there appears to be little correlation between today's returns and previous days' returns. The only substantial correlation is between `Year` and `volume`. By plotting the data, which is ordered chronologically, we see that `volume` is increasing over time. In other words, the average number of shares traded daily increased from 2001 to 2005.

```{r chunk3}
attach(Smarket)
plot(Volume)
```


## Logistic Regression

Next, we will fit a logistic regression model in order to predict `direction` using `lagone` through `lagfive` and `volume`. The `glm()` function  can be used to fit many types of generalized linear models, including logistic regression.
 The syntax of the `glm()` function is similar to that of  `lm()`,  except that we must pass in the argument `family = binomial` in order to tell `R` to run a logistic regression rather than some other type of generalized linear model.


```{r chunk4}
glm.fits <- glm(
    Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume,
    data = Smarket, family = binomial
  )
summary(glm.fits)
```

The smallest $p$-value here is associated with `lagone`. The negative coefficient for this predictor suggests that if the market had a positive return yesterday, then it is less likely to go up today. However, at a value of $0.15$, the $p$-value is still relatively large, and so there is no clear evidence of a real association between `lagone` and `direction`. 

 We use the `coef()` function in order to access just the coefficients for this fitted model. We can also use the `summary()` function to access  particular aspects of the fitted model, such as the $p$-values for the coefficients.

```{r chunk5}
coef(glm.fits)
summary(glm.fits)$coef
summary(glm.fits)$coef[, 4]
```



The `predict()` function can be used to predict the probability that the market will go up, given values of the predictors. The `type = "response"` option tells `R` to output probabilities of the form $P(Y=1|X)$, as opposed to other information such as the logit. If no data set is supplied to the `predict()` function,
then  the probabilities are computed for the training data that was used to fit the logistic regression model. Here we have printed only the first ten probabilities. We know that these values correspond to the probability of the market going up, rather than down, because the `contrasts()` function indicates that `R` has created a dummy variable with a 1 for `Up`.

```{r chunk6}
glm.probs <- predict(glm.fits, type = "response")
glm.probs[1:10]
contrasts(Direction)
```

In order to make a prediction as to whether the market will go up or down on a particular day, we must convert these predicted probabilities into class labels, `Up` or `Down`.
The following two commands create a vector of class predictions based on whether the predicted probability of a market increase is greater than or less than $0.5$.

```{r chunk7}
glm.pred <- rep("Down", 1250)
glm.pred[glm.probs > .5] = "Up"
```

The first command creates a vector of 1,250  `Down` elements. The second line transforms to `Up` all of the elements for which the predicted probability of a market increase exceeds $0.5$. Given these predictions, the `table()` function can be used to produce a confusion matrix  in order to determine how many observations were correctly or incorrectly classified.

```{r chunk8}
table(glm.pred, Direction)
(507 + 145) / 1250
mean(glm.pred == Direction)
```

The diagonal elements of the confusion matrix indicate correct predictions, while the off-diagonals represent incorrect predictions. Hence our model correctly predicted that the market would go up on $507$ days and that it would go down on $145$ days, for a total of $507+145 = 652$ correct predictions. The `mean()` function can be used to compute the fraction of days for which the prediction was correct. In this case, logistic regression correctly predicted the movement of the market  $52.2$ \% of the time.

At first glance, it appears that the logistic regression model is working a little better than random guessing. However, this result is misleading because we trained and tested the  model on the same set of $1,250$ observations. In other words,
 $100\%-52.2\%=47.8\%$, is the *training* error  rate. As we have seen previously, the training error rate is often overly optimistic---it tends to underestimate the  test error rate.  In order to better assess the accuracy of the logistic regression model in this setting, we can fit the model using part of the data, and then examine how well it predicts the *held out* data.
This will yield a more realistic error rate, in the sense that in practice we will be interested in our model's performance not on the data that we used to fit the model, but rather on days in the future for which the market's movements are unknown.

To implement this strategy, we will first create a vector corresponding to the observations from 2001 through 2004. We will then use this vector to create a held out data set of observations from 2005.

```{r chunk9}
train <- (Year < 2005)
Smarket.2005 <- Smarket[!train, ]
dim(Smarket.2005)
Direction.2005 <- Direction[!train]
```

The object `train` is a vector of $1{,}250$ elements, corresponding to the observations in our data set. The elements of the vector that correspond to observations that occurred before 2005 are set to `TRUE`, whereas those that correspond to observations  in 2005 are set to `FALSE`.
The object `train` is a  *Boolean* vector, since its elements are `TRUE` and `FALSE`.
Boolean vectors can be used to obtain a subset of the rows or columns of a matrix. For instance, the command `Smarket[train, ]` would pick out a submatrix of the stock market data set, corresponding only to the dates before 2005, since
those are the ones for which the elements of `train` are `TRUE`.
The `!` symbol can be used to reverse all of the elements of  a Boolean vector. That is, `!train` is a vector similar to `train`, except that  the elements that are `TRUE` in `train` get swapped to `FALSE` in `!train`, and the elements that are `FALSE`
in `train` get swapped to `TRUE` in `!train`. Therefore, `Smarket[!train, ]` yields a submatrix of the stock market data containing only
the observations for which `train` is `FALSE`---that is, the observations with dates in 2005. The output above indicates that there are 252 such observations.

We now fit a logistic regression model using only the subset of the observations that correspond to dates before 2005, using the `subset` argument. We then obtain predicted probabilities of the stock market going up for each of the days in our test set---that is, for the days in 2005.

```{r chunk10}
glm.fits <- glm(
    Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume,
    data = Smarket, family = binomial, subset = train
  )
glm.probs <- predict(glm.fits, Smarket.2005,
    type = "response")
```

Notice that we have trained and tested our model on two completely separate data sets: training was performed using only the dates before 2005, and testing was performed using only the dates in 2005.
Finally, we compute the predictions for 2005 and compare them to the actual movements of the market over that time period.

```{r chunk11}
glm.pred <- rep("Down", 252)
glm.pred[glm.probs > .5] <- "Up"
table(glm.pred, Direction.2005)
mean(glm.pred == Direction.2005)
mean(glm.pred != Direction.2005)
```

The `!=` notation means *not equal to*, and  so the last command computes the  test set error rate. The results are rather disappointing: the test error rate is $52$ \%, which is worse than random guessing! Of course this result is not all that surprising, given that one would not generally expect to be able to  use previous days' returns to predict future market performance. (After all, if it were possible to do so, then the authors of this book would be out striking it rich rather than writing a statistics textbook.)

We recall that the logistic regression model had very underwhelming $p$-values associated with all of the predictors, and that the smallest $p$-value, though not very small, corresponded to `lagone`. Perhaps by removing the variables that appear not  to be helpful in
predicting
`direction`, we can obtain a more effective model. After all, using predictors that have no relationship with the response tends to cause a deterioration in the test error rate (since such predictors cause an increase in variance without a corresponding decrease in bias), and so removing such predictors may in turn yield an improvement.
 Below we have refit the logistic regression using just `lagone` and `lagtwo`, which seemed to have the highest predictive power in the original logistic regression model.

```{r chunk12}
glm.fits <- glm(Direction ~ Lag1 + Lag2, data = Smarket,
    family = binomial, subset = train)
glm.probs <- predict(glm.fits, Smarket.2005,
    type = "response")
glm.pred <- rep("Down", 252)
glm.pred[glm.probs > .5] <- "Up"
table(glm.pred, Direction.2005)
mean(glm.pred == Direction.2005)
106 / (106 + 76)
```

Now the results appear to be a little better: $56\%$ of the daily movements have been correctly predicted. It is worth noting that in this case, a much simpler strategy of predicting that the market will increase every day will also be correct $56\%$ of the time! Hence, in terms of overall error rate, the logistic regression method is no better than the naive approach. However, the confusion matrix shows that on days when logistic regression predicts an increase in the market, it has a $58\%$ accuracy rate. This suggests a possible trading strategy of buying on days when the model predicts an increasing market, and avoiding trades on days when a decrease is predicted. Of course one would need to investigate more carefully whether this small improvement was real or just due to random chance.


Suppose that we want to predict the returns associated with particular values of `lagone` and `lagtwo`. In particular, we want to predict `direction` on a day when `lagone` and `lagtwo` equal 1.2 and~1.1, respectively, and  on a day when they equal  1.5 and $-$0.8.
We do this using the `predict()` function.

```{r chunk13}
predict(glm.fits,
    newdata =
      data.frame(Lag1 = c(1.2, 1.5),  Lag2 = c(1.1, -0.8)),
    type = "response"
  )
```

Now, we plot the decision boundary of the logistic classifier overlayed with data. We do this by set a grid and predict classification values on the grid. We first set up the grid.
```{r grid}
color <- c('#e66101','#5e3c99','#fdb863','#b2abd2')
grid_n <- 400
grid_L1 <- seq(from = min(Smarket[["Lag1"]]), to = max(Smarket[["Lag1"]]),
               length.out = grid_n)
grid_L2 <- seq(from = min(Smarket[["Lag2"]]), to = max(Smarket[["Lag2"]]),
               length.out = grid_n)
grid_L <- expand.grid(Lag1 = grid_L1, Lag2 = grid_L2)
```

Then, we predict classification values on the grid, which will determine the decision boundary. Then we plot the decision boundary using the `.filled.contour()` function and overlay data on it.
```{r plot_glm}
p_grid_glm <- matrix(
  predict(glm.fits, newdata = grid_L, type = "response") > 0.5,
  nrow = grid_n, ncol = grid_n)

plot(NA, main = "Logistic Classifier", xlab = "Lag1", ylab = "Lag2",
      xlim = range(grid_L1), ylim = range(grid_L2))
.filled.contour(x = grid_L1, y = grid_L2, z = p_grid_glm, levels = c(0, 0.5, 1),
      col = color[3:4])
points(Smarket[["Lag1"]], Smarket[["Lag2"]],
      col = color[Smarket[["Direction"]]], pch = 16)
legend("topright", legend = c("Down", "Up"), col = color[1:2], pch = 16)
```

You can see that the decision boundary is linear.

## Linear Discriminant Analysis

Now we will perform LDA on the `Smarket` data. In `R`, we fit an LDA model using the  `lda()` function, which is part of the `MASS` library. Notice that the syntax for the `lda()` function is identical to that of `lm()`, and to that of `glm()` except for the absence of the `family` option. We fit the model using only the observations before 2005.

```{r chunk14}
library(MASS)
lda.fit <- lda(Direction ~ Lag1 + Lag2, data = Smarket,
    subset = train)
lda.fit
plot(lda.fit)
```

The LDA output indicates that $\hat\pi_1=0.492$ and $\hat\pi_2=0.508$; in other words, $49.2$ \% of the training observations correspond to days during which the market went down.
 It also provides the group means; these are the average of each predictor within each class, and are used by LDA as estimates of $\mu_k$.
 These suggest that there is a tendency for the previous 2~days' returns to be negative on days when the market increases, and a tendency for the previous days' returns to be positive on days when the market declines.
  The *coefficients of linear discriminants* output provides the linear combination of `lagone` and `lagtwo` that are used to form the LDA decision rule. In other words, these are the multipliers of the elements of $X=x$ in (4.24).
  If $-0.642 \times \text{lagone} - 0.514 \times \text{lagtwo}$ is large, then the LDA classifier will predict a market increase, and if it is small, then the LDA classifier will predict a market decline.

  The `plot()` function produces plots of the *linear discriminants*, obtained by computing $-0.642 \times \text{lagone} - 0.514 \times \text{lagtwo}$ for each of the training observations. The `Up` and `Down` observations are displayed separately.




The `predict()` function returns a list with three elements. The first element,  `class`,  contains LDA's predictions about the movement of the market. The second element, `posterior`, is a matrix whose $k$th column contains the posterior probability that the corresponding observation belongs to the $k$th class, computed from (4.15). Finally, `x` contains the linear discriminants, described earlier.

```{r chunk15}
lda.pred <- predict(lda.fit, Smarket.2005)
names(lda.pred)
```

As we observed in Section 4.5, the LDA and logistic regression predictions are almost identical.

```{r chunk16}
lda.class <- lda.pred$class
table(lda.class, Direction.2005)
mean(lda.class == Direction.2005)
```

Applying a $50$ \% threshold to the posterior probabilities allows us to recreate the predictions contained in `lda.pred$class`.

```{r chunk17}
sum(lda.pred$posterior[, 1] >= .5)
sum(lda.pred$posterior[, 1] < .5)
```

Notice that the posterior probability output by the model corresponds to the probability that the market will *decrease*:

```{r chunk18}
lda.pred$posterior[1:20, 1]
lda.class[1:20]
```

If we wanted to use a posterior probability threshold other than $50$ \% in order to make predictions, then we could easily do so. For instance, suppose that we wish to predict a market decrease only if we are very certain that the market will indeed decrease on that day---say, if the posterior probability
is at least $90$ \%.

```{r chunk19}
sum(lda.pred$posterior[, 1] > .9)
```

No days in 2005 meet that threshold! In fact, the greatest posterior probability of decrease in all of 2005 was $52.02$ \%.

We again plot the decision boundary of the logistic classifier overlayed with data. We do this by set a grid and predict classification values on the grid. We use the grid computed before. Then, we predict classification values on the grid, which will determine the decision boundary. Then we plot the decision boundary using the `.filled.contour()` function and overlay data on it.
```{r plot_lda}
p_grid_lda <- matrix(as.numeric(predict(lda.fit, newdata = grid_L)[["class"]]),
                     nrow = grid_n, ncol = grid_n)

plot(NA, main = "Linear Discriminant Analysis", xlab = "Lag1", ylab = "Lag2",
      xlim = range(grid_L1), ylim = range(grid_L2))
.filled.contour(x = grid_L1, y = grid_L2, z = p_grid_lda, levels = c(1, 1.5, 2),
      col = color[3:4])
points(Smarket[["Lag1"]], Smarket[["Lag2"]],
      col = color[Smarket[["Direction"]]], pch = 16)
points(lda.fit[["means"]], pch = "+", cex = 3, col = color[3:4])
legend("topright", legend = c("Down", "Up"), col = color[1:2], pch = 16)
```

You can see that the decision boundary is linear. Two crosses are the estimated centers for Gaussian distributions for each class. Note that the decision boundary need not pass between the centers.

## Quadratic Discriminant Analysis

We will now fit a QDA model to the `Smarket` data. QDA is  implemented in `R` using the `qda()` function, which is also part of the `MASS` library. The syntax is identical to that of `lda()`.

```{r chunk20}
qda.fit <- qda(Direction ~ Lag1 + Lag2, data = Smarket,
    subset = train)
qda.fit
```

The output  contains the group means. But it does not contain the coefficients of the linear discriminants, because the QDA classifier involves  a quadratic, rather than a linear, function of the predictors. The `predict()` function works in exactly the same fashion as for LDA.

```{r chunk21}
qda.class <- predict(qda.fit, Smarket.2005)$class
table(qda.class, Direction.2005)
mean(qda.class == Direction.2005)
```

Interestingly, the QDA predictions are accurate almost $60$ \% of the time, even though the 2005 data was not used  to fit the model. This level of accuracy is quite impressive for stock market data, which is known to be quite hard to model accurately.
This  suggests that the quadratic form assumed by QDA may capture the true relationship more accurately than the linear forms assumed by LDA and logistic regression.
 However, we recommend evaluating this method's performance on a larger test set  before betting  that this approach will consistently beat the market!
 
We again plot the decision boundary of the quadratic discriminant analysis overlayed with data. We do this by set a grid and predict classification values on the grid. We use the grid computed before. Then, we predict classification values on the grid, which will determine the decision boundary. Then we plot the decision boundary using the `.filled.contour()` function and overlay data on it.
```{r plot_qda}
p_grid_qda <- matrix(as.numeric(predict(qda.fit, newdata = grid_L)[["class"]]),
                     nrow = grid_n, ncol = grid_n)

plot(NA, main = "Quadratic Discriminant Analysis", xlab = "Lag1", ylab = "Lag2",
     xlim = range(grid_L1), ylim = range(grid_L2))
.filled.contour(x = grid_L1, y = grid_L2, z = p_grid_qda, levels = c(1, 1.5, 2),
                col = color[3:4])
points(Smarket[["Lag1"]], Smarket[["Lag2"]],
       col = color[Smarket[["Direction"]]], pch = 16)
points(qda.fit[["means"]], pch = "+", cex = 3, col = color[3:4])
legend("topright", legend = c("Down", "Up"), col = color[1:2], pch = 16)
```

You can see that the decision boundary is quadratic. Two crosses are the estimated centers for Gaussian distributions for each class. Note that the decision boundary need not pass between the centers.
