Lab5

13

library(MASS)
library(class)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::select() masks MASS::select()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(corrplot)

## Warning: package 'corrplot' was built under R version 4.4.3

## corrplot 0.95 loaded

library(ISLR2)

## Warning: package 'ISLR2' was built under R version 4.4.3

## 
## Attaching package: 'ISLR2'
## 
## The following object is masked from 'package:MASS':
## 
##     Boston

library(e1071)

## Warning: package 'e1071' was built under R version 4.4.3

data(Weekly)

summary(Weekly)

##       Year           Lag1               Lag2               Lag3         
##  Min.   :1990   Min.   :-18.1950   Min.   :-18.1950   Min.   :-18.1950  
##  1st Qu.:1995   1st Qu.: -1.1540   1st Qu.: -1.1540   1st Qu.: -1.1580  
##  Median :2000   Median :  0.2410   Median :  0.2410   Median :  0.2410  
##  Mean   :2000   Mean   :  0.1506   Mean   :  0.1511   Mean   :  0.1472  
##  3rd Qu.:2005   3rd Qu.:  1.4050   3rd Qu.:  1.4090   3rd Qu.:  1.4090  
##  Max.   :2010   Max.   : 12.0260   Max.   : 12.0260   Max.   : 12.0260  
##       Lag4               Lag5              Volume            Today         
##  Min.   :-18.1950   Min.   :-18.1950   Min.   :0.08747   Min.   :-18.1950  
##  1st Qu.: -1.1580   1st Qu.: -1.1660   1st Qu.:0.33202   1st Qu.: -1.1540  
##  Median :  0.2380   Median :  0.2340   Median :1.00268   Median :  0.2410  
##  Mean   :  0.1458   Mean   :  0.1399   Mean   :1.57462   Mean   :  0.1499  
##  3rd Qu.:  1.4090   3rd Qu.:  1.4050   3rd Qu.:2.05373   3rd Qu.:  1.4050  
##  Max.   : 12.0260   Max.   : 12.0260   Max.   :9.32821   Max.   : 12.0260  
##  Direction 
##  Down:484  
##  Up  :605  
##            
##            
##            
##

str(Weekly)

## 'data.frame':    1089 obs. of  9 variables:
##  $ Year     : num  1990 1990 1990 1990 1990 1990 1990 1990 1990 1990 ...
##  $ Lag1     : num  0.816 -0.27 -2.576 3.514 0.712 ...
##  $ Lag2     : num  1.572 0.816 -0.27 -2.576 3.514 ...
##  $ Lag3     : num  -3.936 1.572 0.816 -0.27 -2.576 ...
##  $ Lag4     : num  -0.229 -3.936 1.572 0.816 -0.27 ...
##  $ Lag5     : num  -3.484 -0.229 -3.936 1.572 0.816 ...
##  $ Volume   : num  0.155 0.149 0.16 0.162 0.154 ...
##  $ Today    : num  -0.27 -2.576 3.514 0.712 1.178 ...
##  $ Direction: Factor w/ 2 levels "Down","Up": 1 1 2 2 2 1 2 2 2 1 ...

a. Produce some numerical and graphical summaries of the Weekly data. Do there appear to be any patterns?

pairs(Weekly)  # Scatterplot matrix

cor(Weekly[, -9])

##               Year         Lag1        Lag2        Lag3         Lag4
## Year    1.00000000 -0.032289274 -0.03339001 -0.03000649 -0.031127923
## Lag1   -0.03228927  1.000000000 -0.07485305  0.05863568 -0.071273876
## Lag2   -0.03339001 -0.074853051  1.00000000 -0.07572091  0.058381535
## Lag3   -0.03000649  0.058635682 -0.07572091  1.00000000 -0.075395865
## Lag4   -0.03112792 -0.071273876  0.05838153 -0.07539587  1.000000000
## Lag5   -0.03051910 -0.008183096 -0.07249948  0.06065717 -0.075675027
## Volume  0.84194162 -0.064951313 -0.08551314 -0.06928771 -0.061074617
## Today  -0.03245989 -0.075031842  0.05916672 -0.07124364 -0.007825873
##                Lag5      Volume        Today
## Year   -0.030519101  0.84194162 -0.032459894
## Lag1   -0.008183096 -0.06495131 -0.075031842
## Lag2   -0.072499482 -0.08551314  0.059166717
## Lag3    0.060657175 -0.06928771 -0.071243639
## Lag4   -0.075675027 -0.06107462 -0.007825873
## Lag5    1.000000000 -0.05851741  0.011012698
## Volume -0.058517414  1.00000000 -0.033077783
## Today   0.011012698 -0.03307778  1.000000000

Volume is strongly positively correlated with Year. Other correlations are week, but Lag1 is negatively correlated with Lag2 but positively correlated with Lag3.

b. Use the full data set to perform a logistic regression with Direction as the response and the five lag variables plus Volume as predictors. Use the summary function to print the results. Do any of the predictors appear to be statistically significant? If so, which ones?

logistic_model <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, 
                      data = Weekly, family = binomial)
summary(logistic_model)

## 
## Call:
## glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + 
##     Volume, family = binomial, data = Weekly)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.26686    0.08593   3.106   0.0019 **
## Lag1        -0.04127    0.02641  -1.563   0.1181   
## Lag2         0.05844    0.02686   2.175   0.0296 * 
## Lag3        -0.01606    0.02666  -0.602   0.5469   
## Lag4        -0.02779    0.02646  -1.050   0.2937   
## Lag5        -0.01447    0.02638  -0.549   0.5833   
## Volume      -0.02274    0.03690  -0.616   0.5377   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1496.2  on 1088  degrees of freedom
## Residual deviance: 1486.4  on 1082  degrees of freedom
## AIC: 1500.4
## 
## Number of Fisher Scoring iterations: 4

Intercept: The intercept is statistically significant (p=0.0019p = 0.0019p=0.0019), meaning that even when all predictors are zero, the model has a significant baseline effect.

Lag1, Lag3, Lag4, Lag5, and Volume: These variables are not statistically significant as their p-values are greater than 0.05.

Lag2: This is the only statistically significant predictor (p=0.0296p = 0.0296p=0.0296), meaning it has a meaningful impact on predicting the Direction

c. Compute the confusion matrix and overall fraction of correct predictions. Explain what the confusion matrix is telling you about the types of mistakes made by logistic regression.

logistic_pred <- ifelse(predict(logistic_model, type = "response") > 0.5, "Up", "Down")

conf_matrix_logistic <- table(Predicted = logistic_pred, Actual = Weekly$Direction)

accuracy_logistic <- sum(diag(conf_matrix_logistic)) / sum(conf_matrix_logistic)
conf_matrix_logistic

##          Actual
## Predicted Down  Up
##      Down   54  48
##      Up    430 557

accuracy_logistic

## [1] 0.5610652

Observations:

The model predicts “Up” very frequently, as seen by the large number of false negatives (430), meaning it struggles to identify “Down” correctly.
Accuracy is only slightly better than random guessing (50%), suggesting the model might not be very effective.
Since Lag2 was the only statistically significant variable, a model with just Lag2 might perform similarly with fewer parameters.
The overall fraction of correct predictions is 0.56. Although logistic regression correctly predicts upwards movements well, it incorrectly predicts most downwards movements as up.

d. Now fit the logistic regression model using a training data period from 1990 to 2008, with Lag2 as the only predictor. Compute the confusion matrix and the overall fraction of correct predictions for the held out data (that is, the data from 2009 and 2010).

train <- Weekly$Year < 2009
test <- !train

logistic_model2 <- glm(Direction ~ Lag2, data = Weekly, subset = train, family = binomial)
logistic_pred2 <- ifelse(predict(logistic_model2, Weekly[test,], type = "response") > 0.5, "Up", "Down")

conf_matrix_logistic2 <- table(Predicted = logistic_pred2, Actual = Weekly$Direction[test])
accuracy_logistic2 <- sum(diag(conf_matrix_logistic2)) / sum(conf_matrix_logistic2)

conf_matrix_logistic2

##          Actual
## Predicted Down Up
##      Down    9  5
##      Up     34 56

accuracy_logistic2

## [1] 0.625

Observations:

Accuracy: 62.5%: The model correctly predicts the direction 62.5% of the time.
Compared to the full dataset model (56.1% accuracy), this model is slightly better at predicting new data.
The model still struggles with “Down” days, as shown by the relatively high number of false negatives (34).
Since we are only using Lag2 as the predictor, other features or transformations might improve performance.

e. Repeat (d) using LDA.

lda_model <- lda(Direction ~ Lag2, data = Weekly, subset = train)
lda_pred <- predict(lda_model, Weekly[test,])$class
conf_matrix_lda <- table(Predicted = lda_pred, Actual = Weekly$Direction[test])
accuracy_lda <- sum(diag(conf_matrix_lda)) / sum(conf_matrix_lda)
conf_matrix_lda

##          Actual
## Predicted Down Up
##      Down    9  5
##      Up     34 56

accuracy_lda

## [1] 0.625

f. Repeat (d) using QDA.

qda_model <- qda(Direction ~ Lag2, data = Weekly, subset = train)
qda_pred <- predict(qda_model, Weekly[test,])$class
conf_matrix_qda <- table(Predicted = qda_pred, Actual = Weekly$Direction[test])
accuracy_qda <- sum(diag(conf_matrix_qda)) / sum(conf_matrix_qda)
conf_matrix_qda

##          Actual
## Predicted Down Up
##      Down    0  0
##      Up     43 61

accuracy_qda

## [1] 0.5865385

g. Repeat (d) using KNN with K=1.

x_train <- Weekly[train, "Lag2", drop = FALSE]
x_test <- Weekly[test, "Lag2", drop = FALSE]
y_train <- Weekly$Direction[train]
y_test <- Weekly$Direction[test]

knn_pred <- knn(x_train, x_test, y_train, k = 1)
conf_matrix_knn <- table(Predicted = knn_pred, Actual = y_test)
accuracy_knn <- sum(diag(conf_matrix_knn)) / sum(conf_matrix_knn)
conf_matrix_knn

##          Actual
## Predicted Down Up
##      Down   21 29
##      Up     22 32

accuracy_knn

## [1] 0.5096154

h. Repeat (d) using naive Bayes.

nb_model <- naiveBayes(Direction ~ Lag2, data = Weekly, subset = train)
nb_pred <- predict(nb_model, Weekly[test,])
conf_matrix_nb <- table(Predicted = nb_pred, Actual = Weekly$Direction[test])
accuracy_nb <- sum(diag(conf_matrix_nb)) / sum(conf_matrix_nb)
conf_matrix_nb

##          Actual
## Predicted Down Up
##      Down    0  0
##      Up     43 61

accuracy_nb

## [1] 0.5865385

i. Which of these methods appears to provide the best results on this data?

accuracies <- c(Logistic = accuracy_logistic2, LDA = accuracy_lda, 
                QDA = accuracy_qda, KNN = accuracy_knn, NaiveBayes = accuracy_nb)
accuracies

##   Logistic        LDA        QDA        KNN NaiveBayes 
##  0.6250000  0.6250000  0.5865385  0.5096154  0.5865385

Logistic Regression and LDA performed the best (62.5%), meaning they are making similar predictions.

QDA and Naïve Bayes (both ~58.65%) performed slightly worse, indicating that a quadratic decision boundary or probabilistic approach might not be capturing the pattern well.

KNN (K=1) performed the worst (50%), suggesting that a simple nearest-neighbor approach is not effective here. Trying different values of K (e.g., 3, 5, 10) might help.

k. Experiment with different combinations of predictors, including possible transformations and interactions, for each of the methods. Report the variables, method, and associated confusion matrix that appears to provide the best results on the held out data. Note that you should also experiment with values for K in the KNN classifier.

logistic_model3 <- glm(Direction ~ Lag2 + Lag2^2, data = Weekly, subset = train, family = binomial)
logistic_pred3 <- ifelse(predict(logistic_model3, Weekly[test,], type = "response") > 0.5, "Up", "Down")
conf_matrix_logistic3 <- table(Predicted = logistic_pred3, Actual = Weekly$Direction[test])
accuracy_logistic3 <- sum(diag(conf_matrix_logistic3)) / sum(conf_matrix_logistic3)

# Experiment with different K values for KNN
knn_pred5 <- knn(x_train, x_test, y_train, k = 5)
conf_matrix_knn5 <- table(Predicted = knn_pred5, Actual = y_test)
accuracy_knn5 <- sum(diag(conf_matrix_knn5)) / sum(conf_matrix_knn5)
conf_matrix_knn5

##          Actual
## Predicted Down Up
##      Down   16 20
##      Up     27 41

accuracy_knn5

## [1] 0.5480769

set.seed(1)
res <- sapply(1:30, function(k) {
  fit <- knn(
    Weekly[train, 2:4, drop = FALSE],
    Weekly[!train, 2:4, drop = FALSE],
    Weekly$Direction[train],
    k = k
  )
  mean(fit == Weekly[!train, ]$Direction)
})
plot(1:30, res, type = "o", xlab = "k", ylab = "Fraction correct")

# Final accuracy comparison
accuracies_extended <- c(Logistic = accuracy_logistic3, LDA = accuracy_lda, 
                         QDA = accuracy_qda, KNN_1 = accuracy_knn, KNN_5 = accuracy_knn5, NaiveBayes = accuracy_nb)
accuracies_extended

##   Logistic        LDA        QDA      KNN_1      KNN_5 NaiveBayes 
##  0.6250000  0.6250000  0.5865385  0.5096154  0.5480769  0.5865385

Observations:

Logistic Regression and LDA remain the best models (62.5%).
QDA and Naïve Bayes (58.65%) are slightly worse.
KNN (K=1) performed very poorly (50%), but KNN (K=5) improved slightly (53.85%).

-   Since KNN is very sensitive to the choice of K, testing values like K=10, K=15, k=20 might be worthwhile.

-   KNN using the first 3 Lag variables performs marginally better than logistic regression with Lag2 if we tune K to be K=26.

14.

data(Auto)

summary(Auto)

##       mpg          cylinders      displacement     horsepower        weight    
##  Min.   : 9.00   Min.   :3.000   Min.   : 68.0   Min.   : 46.0   Min.   :1613  
##  1st Qu.:17.00   1st Qu.:4.000   1st Qu.:105.0   1st Qu.: 75.0   1st Qu.:2225  
##  Median :22.75   Median :4.000   Median :151.0   Median : 93.5   Median :2804  
##  Mean   :23.45   Mean   :5.472   Mean   :194.4   Mean   :104.5   Mean   :2978  
##  3rd Qu.:29.00   3rd Qu.:8.000   3rd Qu.:275.8   3rd Qu.:126.0   3rd Qu.:3615  
##  Max.   :46.60   Max.   :8.000   Max.   :455.0   Max.   :230.0   Max.   :5140  
##                                                                                
##   acceleration        year           origin                      name    
##  Min.   : 8.00   Min.   :70.00   Min.   :1.000   amc matador       :  5  
##  1st Qu.:13.78   1st Qu.:73.00   1st Qu.:1.000   ford pinto        :  5  
##  Median :15.50   Median :76.00   Median :1.000   toyota corolla    :  5  
##  Mean   :15.54   Mean   :75.98   Mean   :1.577   amc gremlin       :  4  
##  3rd Qu.:17.02   3rd Qu.:79.00   3rd Qu.:2.000   amc hornet        :  4  
##  Max.   :24.80   Max.   :82.00   Max.   :3.000   chevrolet chevette:  4  
##                                                  (Other)           :365

a. Create a binary variable, mpg01, that contains a 1 if mpg contains a value above its median, and a 0 if mpg contains a value below its median. You can compute the median using the median() function. Note you may find it helpful to use the data.frame() function to create a single data set containing both mpg01 and the other Auto variables.

Auto$mpg01 <- ifelse(Auto$mpg > median(Auto$mpg), 1, 0)

b. Explore the data graphically in order to investigate the association between mpg01 and the other features. Which of the other features seem most likely to be useful in predicting mpg01? Scatterplots and boxplots may be useful tools to answer this question. Describe your findings.

par(mfrow = c(2, 2))
boxplot(Auto$horsepower ~ Auto$mpg01, main = "Horsepower vs mpg01", col = "lightblue")
boxplot(Auto$weight ~ Auto$mpg01, main = "Weight vs mpg01", col = "lightgreen")
boxplot(Auto$acceleration ~ Auto$mpg01, main = "Acceleration vs mpg01", col = "lightcoral")
boxplot(Auto$displacement ~ Auto$mpg01, main = "Displacement vs mpg01", col = "lightyellow")

Observations:

Horsepower vs mpg01:
- Cars with lower fuel efficiency (mpg01 = 0) tend to have higher horsepower.
- Cars with higher fuel efficiency (mpg01 = 1) show significantly lower horsepower.
Weight vs mpg01:
- Heavier cars are associated with lower fuel efficiency (mpg01 = 0).
- Lighter cars are more fuel-efficient (mpg01 = 1).
Acceleration vs mpg01:
- There is no clear distinction in acceleration between the two categories of mpg01. Both groups seem to have overlapping distributions.
Displacement vs mpg01:
- Cars with lower fuel efficiency (mpg01 = 0) tend to have larger engine displacement.
- Cars with higher fuel efficiency (mpg01 = 1) generally have smaller engine displacement.

c. Split the data into a training set and a test set.

set.seed(1)
train_index <- sample(1:nrow(Auto), size = 0.7 * nrow(Auto))
train_data <- Auto[train_index, ]
test_data <- Auto[-train_index, ]

d. Perform LDA on the training data in order to predict mpg01 using the variables that seemed most associated with mpg01 in (b). What is the test error of the model obtained?

predictors <- c("horsepower", "weight", "displacement", "cylinders")

# (d) Perform LDA
lda_fit <- lda(mpg01 ~ horsepower + weight + displacement + cylinders, data = train_data)
lda_pred <- predict(lda_fit, test_data)$class
lda_error <- mean(lda_pred != test_data$mpg01)
print(paste("LDA Test Error:", round(lda_error, 4)))

## [1] "LDA Test Error: 0.1186"

e. Perform QDA on the training data in order to predict mpg01 using the variables that seemed most associated with mpg01 in (b). What is the test error of the model obtained?

qda_fit <- qda(mpg01 ~ horsepower + weight + displacement + cylinders, data = train_data)
qda_pred <- predict(qda_fit, test_data)$class
qda_error <- mean(qda_pred != test_data$mpg01)
print(paste("QDA Test Error:", round(qda_error, 4)))

## [1] "QDA Test Error: 0.1186"

f. Perform logistic regression on the training data in order to predict mpg01 using the variables that seemed most associated with mpg01 in (b). What is the test error of the model obtained?

logit_fit <- glm(mpg01 ~ horsepower + weight + displacement + cylinders, data = train_data, family = binomial)
logit_probs <- predict(logit_fit, test_data, type = "response")
logit_pred <- ifelse(logit_probs > 0.5, 1, 0)
logit_error <- mean(logit_pred != test_data$mpg01)
print(paste("Logistic Regression Test Error:", round(logit_error, 4)))

## [1] "Logistic Regression Test Error: 0.0932"

g. Perform naive Bayes on the training data in order to predict mpg01 using the variables that seemed most associated with mpg01 in (b). What is the test error of the model obtained?

nb_fit <- naiveBayes(mpg01 ~ horsepower + weight + displacement + cylinders, data = train_data)
nb_pred <- predict(nb_fit, test_data)
nb_error <- mean(nb_pred != test_data$mpg01)
print(paste("Naïve Bayes Test Error:", round(nb_error, 4)))

## [1] "Naïve Bayes Test Error: 0.1102"

h. Perform KNN on the training data, with several values of K, in order to predict mpg01. Use only the variables that seemed most associated with mpg01 in (b). What test errors do you obtain? Which value of K seems to perform the best on this data set?

set.seed(1)
train_X <- scale(train_data[, predictors])  # Standardizing the predictors
test_X <- scale(test_data[, predictors])
train_Y <- train_data$mpg01
test_Y <- test_data$mpg01

knn_errors <- sapply(1:20, function(k) {
  knn_pred <- knn(train_X, test_X, train_Y, k = k)
  mean(knn_pred != test_Y)
})

# Find best K
best_k <- which.min(knn_errors)
print(paste("Best K:", best_k, "with Test Error:", round(knn_errors[best_k], 4)))

## [1] "Best K: 12 with Test Error: 0.1102"

#plot KNN test errors
plot(1:20, knn_errors, type = "o", col = "blue", xlab = "K", ylab = "Test Error", main = "KNN Test Error for Different K")

For the models tested here, K=12 appears to perform best.

Lab5

Kiran

2025-03-11

13

Observations:

Observations:

Observations:

14.

Observations: