## Chapter 4.6 Lab: LDA, QDA, and KNN

   
### Linear Discriminant Analysis method:

We aim to use discriminant analysis to predict the stock price change:

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

attach(Smarket)

train <- (Year < 2005) # split data into training and testing data
test <- Smarket[!train, ]  # test data
dim(test) # 252 * 9 for the testing data
Direction.test <- Direction[!train] # test data directions

lda.fit <- lda(Direction ~ Lag1 + Lag2, data = Smarket, subset = train)
lda.fit  #if −0.642 × Lag1 − 0.514 × Lag2 is large, we predict up.
lda.pred <- predict(lda.fit, test)
names(lda.pred)
lda.class <- lda.pred$class
table(lda.class, Direction.test)
paste('Correct prediction rate is',mean(lda.class == Direction.test))  #similar results to logistic regression
```

Graphically what we are trying to achieve, this should also explain why it is not working well, the factors are too clustered together.

```{r}
library(ggplot2)

# Assuming Smarket data is already loaded and 'train' is defined
lda.fit <- lda(Direction ~ Lag1 + Lag2, data = Smarket, subset = train)

# Creating a grid to plot decision boundary
xrange <- range(Smarket$Lag1)
yrange <- range(Smarket$Lag2)
grid <- expand.grid(Lag1 = seq(from = xrange[1], to = xrange[2], length.out = 200),
                    Lag2 = seq(from = yrange[1], to = yrange[2], length.out = 200))

# Predicting class for each point in the grid
lda.pred <- predict(lda.fit, newdata = grid)
class.pred <- lda.pred$class

# Convert to data frame for ggplot
grid['class'] <- as.factor(class.pred)
data_to_plot <- Smarket[train, c('Lag1', 'Lag2', 'Direction')]

# Plot
ggplot(data_to_plot, aes(x = Lag1, y = Lag2, color = Direction)) + 
  geom_point(alpha = 0.8) + 
  scale_color_manual(values = c('blue', 'red')) +
  geom_contour(data = grid, aes(z = as.numeric(class)), breaks = 1.5, color = 'black') +
  ggtitle('LDA Decision Boundary for Smarket Data')
```

## Quadratic Discriminant Analysis

similar work on a QDA format:

```{r}
###   Quadratic Discriminant Analysis

qda.fit <- qda(Direction ~ Lag1 + Lag2, data = Smarket, subset = train)
qda.fit
qda.class <- predict(qda.fit, test)$class
table(qda.class, Direction.test)
paste('Correct prediction rate is',mean(qda.class == Direction.test))
```

it is getting better! let's see what the prediction curve is doing graphically we are trying to:

```{r}
# Assuming Smarket data is already loaded and 'train' is defined
qda.fit <- qda(Direction ~ Lag1 + Lag2, data = Smarket, subset = train)

# Creating a grid to plot decision boundary
xrange <- range(Smarket$Lag1)
yrange <- range(Smarket$Lag2)
grid <- expand.grid(Lag1 = seq(from = xrange[1], to = xrange[2], length.out = 200),
                    Lag2 = seq(from = yrange[1], to = yrange[2], length.out = 200))

# Predicting class for each point in the grid
qda.pred <- predict(qda.fit, newdata = grid)
class.pred <- qda.pred$class

# Convert to data frame for ggplot
grid['class'] <- as.factor(class.pred)
data_to_plot <- Smarket[train, c('Lag1', 'Lag2', 'Direction')]

# Plot
ggplot(data_to_plot, aes(x = Lag1, y = Lag2, color = Direction)) + 
  geom_point(alpha = 0.8) + 
  scale_color_manual(values = c('blue', 'red')) +
  geom_contour(data = grid, aes(z = as.numeric(class)), breaks = 1.5, color = 'black') +
  ggtitle('QDA Decision Boundary for Smarket Data')
```

## K nearest Neighbors

### For K = 1:

```{r}
library(class)
train.X <- cbind(Lag1,Lag2)[train, ] #X for training data
test.X <- cbind(Lag1,Lag2)[!train, ] #X for testing data
train.Direction <- Direction[train] #Y for training data

set.seed(1) # to ensure the reproducibility, set a random seed beforehand
knn.pred <- knn(train.X, test.X, train.Direction, k = 1)
table(knn.pred, Direction.test)
paste('Correct prediction rate is', mean(knn.pred == Direction.test))
```

it is just as good as random guessing, let's take a look at why:

```{r}
train_data <- data.frame(Lag1 = train.X[,1], Lag2 = train.X[,2], Direction = train.Direction)
test_data <- data.frame(Lag1 = test.X[,1], Lag2 = test.X[,2], Prediction = knn.pred)
# Plot
ggplot() +
  geom_point(data = train_data, aes(x = Lag1, y = Lag2, color = Direction), shape = 1) +
  geom_point(data = test_data, aes(x = Lag1, y = Lag2, color = Prediction), shape = 4) +
  scale_color_manual(values = c('blue', 'red')) +
  labs(color = 'Class', title = 'KNN Classification with k = 1') +
  theme_minimal()
```

for the exterior part it was doing with less confidence, and the interior part is dealing with too much complexity, let's see if increasing K will get better:

### For K = 3

```{r}
knn.pred <- knn(train.X, test.X, train.Direction, k = 3)
table(knn.pred, Direction.test)
paste('Correct prediction rate is',mean(knn.pred == Direction.test))
```

getting slightly better

## Results Summary:

| Algorithms                      | Correct Prediction Rate |
|---------------------------------|-------------------------|
| GLM with every variable         | 0.480                   |
| GLM with Lag1 and Lag2          | 0.560                   |
| LDA with Lag1 and Lag2          | 0.560                   |
| QDA with Lag1 and Lag2          | 0.600                   |
| KNN (K = 1 ) with Lag1 and Lag2 | 0.500                   |
| KNN (K = 3) with Lag1 and Lag2  | 0.536                   |

: Summary of Prediction Accuracy for Different Methods

Notice this still largely depends on the parameter choice, students could feel free to test other combinations.

## An Application on Other data --- Caravan Insurance Data

```{r, message = FALSE}
### An Application to Caravan Insurance Data
?Caravan
dim(Caravan)
attach(Caravan)
summary(Purchase)
paste('yes ratio of the data is', sum(Purchase == 'Yes')/ length(Purchase))  # the ratio of Yes is very small, i.e. the two classes are highly imbalanced

```

Let's perform a train test split, with scaling the train matrix at the same variance, we will use the first 1000 data as testing dataset

```{r}
standardized.X <- scale(Caravan[,-86])  # the 86th column is the Purchase outcome
# st
var(Caravan[,1]) # the variance is too high..
var(Caravan[,2]) # and this one is a bit too low...
var(standardized.X[,1]) # now this is better! 
var(standardized.X[,2])

# perform a train test split
test_ind <- 1:1000 # we will 
train.X <- standardized.X[-test_ind,]
train.Y <- Purchase[-test_ind]

test.X <- standardized.X[test_ind,]
test.Y <- Purchase[test_ind]
```

### KNN method with K = 1

```{r}

set.seed(1)
knn.pred <- knn(train.X, test.X, train.Y, k = 1)
table(knn.pred, test.Y)
mean(test.Y != knn.pred)
mean(test.Y != "No")  # Overall, KNN is no better than blindly classifying to "No"

9 / (9 + 68)  # However, it improves the chance of Purchase to 11.7% comparing to 5.8% if we look at the predicted 'Yes' group. (think of a compaignian, this forecasted group might be good to target compare to target a general population)
```

### KNN with K = 3

```{r}
knn.pred <- knn(train.X, test.X, train.Y, k = 3)
table(knn.pred, test.Y)
5 / (5 + 21) # increasing, but with less and less samples in the group


```

### KNN with k = 5

```{r}
knn.pred <- knn(train.X, test.X, train.Y, k = 5)
table(knn.pred,test.Y)
4 / (4 + 11) # even fewer samples


```

### GLM with 0.5 as cutting threshold

```{r}
glm.fit <- glm(Purchase~., data = Caravan, family = binomial, subset = -test_ind)
glm.probs <- predict(glm.fit, Caravan[test_ind,], type="response")
glm.pred <- rep("No", 1000)
glm.pred[glm.probs > 0.5] <- "Yes"
table(glm.pred,test.Y) # not working well as there is no predicted Yes value correctly. This might indicate us to lower the threshold to make the decision


```

### GLM with 0.25 as cutting threshold

```{r}
glm.pred <- rep("No",1000)
glm.pred[glm.probs > 0.25] <- "Yes"
table(glm.pred,test.Y) 
11/(22+11)  # slightly better but still very high FALSE POSITIVE rate
```
