homework03

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.3.2

library(ISLR)
library(corrplot)

## corrplot 0.92 loaded

library(MASS)
library(class)
library(e1071)

Question 13

head(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  
##            
##            
##            
## 

# I only used numerics for better results

cor_matrix <- cor(Weekly[, sapply(Weekly, is.numeric)])

corrplot(cor_matrix, method = "circle")

## Question 13(a)

According to my results, I found out that Volume and Year have a little strong correlation with each other which might suggest that the volume of trading has increased over the years.

Question 13(b)

# fitting the logistic regression model 

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

# printing the summary of the model 

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

According to my results, I saw that variable Lag2 is statistically significant at the 0.05 level. Which means, when Lag2 increases, the log-odds of the market going up vesus down increase. Meaning that there is a positive relationship between Lag2 and the likelihood of positive market direction

Question 13(c)

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

# confusion matrix
confusion_matrix <- table(Predicted = fitted_results, Actual = Weekly$Direction)

print(confusion_matrix)

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

# overall fraction of correct predictions
overall_accuracy <- sum(diag(confusion_matrix)) / sum(confusion_matrix)

# overall accuracy
print(overall_accuracy)

## [1] 0.5610652

According to my results, my overall accuracy is around 56.1%, which means my model predicts the model a little bit more than half of the time:

1. my model predictics more up than down.

Question 13(d)

# Subseting the Weekly data into training and test sets
train <- subset(Weekly, Year < 2009)
test <- subset(Weekly, Year > 2008)

# Fitting the logistic regression model on the training set
train_model <- glm(Direction ~ Lag2, data=train, family=binomial)

# Predicting the Direction on the test set
test$predicted_direction <- predict(train_model, newdata=test, type="response")
test$predicted_class <- ifelse(test$predicted_direction > 0.5, "Up", "Down")

# confusion matrix
confusion_matrix_test <- table(Predicted = test$predicted_class, Actual = test$Direction)

print(confusion_matrix_test)

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

# overall fraction of correct predictions
overall_accuracy_test <- sum(diag(confusion_matrix_test)) / sum(confusion_matrix_test)

# overall accuracy
print(overall_accuracy_test)

## [1] 0.625

According to my results, using Lag2 as my only predictor has an overall accuracy of 62.5%. Which means it can predict the market direction 62.5% of the time for the years 2009 and 2010.

Question 13(e)

train <- subset(Weekly, Year < 2009)
test <- subset(Weekly, Year > 2008)

# Fit the LDA model on the training set
lda_model <- lda(Direction ~ Lag2, data=train)


test$predicted_class <- predict(lda_model, newdata=test)$class

# confusion matrix
confusion_matrix_test <- table(Predicted = test$predicted_class, Actual = test$Direction)

print(confusion_matrix_test)

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

# overall fraction of correct predictions
overall_accuracy_test <- sum(diag(confusion_matrix_test)) / sum(confusion_matrix_test)

# overall accuracy
print(overall_accuracy_test)

## [1] 0.625

According to my results, using the Linear Discriminant Analysis, I see that I also got 62.5% on the test data(years 2009 and 2010) using Lag 2 as the only predictor

Question 13(d)

# Fitting the QDA model on the training data
qda_model <- qda(Direction ~ Lag2, data = train)


qda_predictions <- predict(qda_model, newdata = test)
test$qda_predicted_class <- qda_predictions$class

# confusion matrix for the QDA predictions
confusion_matrix_qda <- table(Predicted = test$qda_predicted_class, Actual = test$Direction)


print(confusion_matrix_qda)

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

# overall fraction of correct predictions for the QDA model
overall_accuracy_qda <- sum(diag(confusion_matrix_qda)) / sum(confusion_matrix_qda)

#overall accuracy for the QDA model
print(overall_accuracy_qda)

## [1] 0.5865385

According to my results, I saw that the accuracy rate dropped to 58.65 which is lower than the logistic regression and LDA models. The model did not predict any down market weeks.

Question 13(d)

train_matrix <- as.matrix(train$Lag2)
test_matrix <- as.matrix(test$Lag2)

# factor
train_direction <- train$Direction

# KNN with K = 1
knn_pred <- knn(train = train_matrix, test = test_matrix, cl = train_direction, k = 1)

# confusion matrix
confusion_matrix_knn <- table(Predicted = knn_pred, Actual = test$Direction)

print(confusion_matrix_knn)

##          Actual
## Predicted Down Up
##      Down   21 30
##      Up     22 31

# overall fraction of correct predictions
overall_accuracy_knn <- sum(diag(confusion_matrix_knn)) / sum(confusion_matrix_knn)


print(overall_accuracy_knn)

## [1] 0.5

Accoring to my results, I got 50% for accuracy rate which is basically a random chance.

Question 13(h)

# Fitting the Naive Bayes model on the training data
nb_model <- naiveBayes(Direction ~ Lag2, data = train)


nb_predictions <- predict(nb_model, newdata = test)

# Creating the confusion matrix for the Naive Bayes predictions
confusion_matrix_nb <- table(Predicted = nb_predictions, Actual = test$Direction)


print(confusion_matrix_nb)

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

# overall fraction of correct predictions for the Naive Bayes model
overall_accuracy_nb <- sum(diag(confusion_matrix_nb)) / sum(confusion_matrix_nb)

# overall accuracy for the Naive Bayes model
print(overall_accuracy_nb)

## [1] 0.5865385

According to my results, I saw that the model predictied around 58.65%, which was higher than the KNN with K = 1

Question 13(i)

According to my results, the logistic regression model and LDA gave me the best results with an accuracy rate of 62.5% and I saw that KNN with K = 1 was the worst with 50% accuracy rate.

Question 13(j)

results <- data.frame(K=integer(), Accuracy=double())

# Assuming 'train_matrix', 'test_matrix', and 'train_direction' are already defined
for (k in 1:20) {
  knn_pred <- knn(train = train_matrix, test = test_matrix, cl = train_direction, k = k)
  confusion_matrix_knn <- table(Predicted = knn_pred, Actual = test$Direction)
  accuracy <- sum(diag(confusion_matrix_knn)) / sum(confusion_matrix_knn)
  results <- rbind(results, data.frame(K=k, Accuracy=accuracy))
}

# Find the K with the highest accuracy
best_model <- results[which.max(results$Accuracy), ]
print(best_model)

##   K  Accuracy
## 2 2 0.6057692

According to my results, among the range of K values that I tested (1-20),

K = 17 gave me the best accuracy.

Question 14

head(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

Question 14(a)

median_mpg <- median(Auto$mpg)

# Create the mpg01 variable based on the median
Auto$mpg01 <- ifelse(Auto$mpg > median_mpg, 1, 0)



str(Auto)

## 'data.frame':    392 obs. of  10 variables:
##  $ mpg         : num  18 15 18 16 17 15 14 14 14 15 ...
##  $ cylinders   : num  8 8 8 8 8 8 8 8 8 8 ...
##  $ displacement: num  307 350 318 304 302 429 454 440 455 390 ...
##  $ horsepower  : num  130 165 150 150 140 198 220 215 225 190 ...
##  $ weight      : num  3504 3693 3436 3433 3449 ...
##  $ acceleration: num  12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
##  $ year        : num  70 70 70 70 70 70 70 70 70 70 ...
##  $ origin      : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ name        : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...
##  $ mpg01       : num  0 0 0 0 0 0 0 0 0 0 ...

# deleting name variable 

cor(Auto[,-9])

##                     mpg  cylinders displacement horsepower     weight
## mpg           1.0000000 -0.7776175   -0.8051269 -0.7784268 -0.8322442
## cylinders    -0.7776175  1.0000000    0.9508233  0.8429834  0.8975273
## displacement -0.8051269  0.9508233    1.0000000  0.8972570  0.9329944
## horsepower   -0.7784268  0.8429834    0.8972570  1.0000000  0.8645377
## weight       -0.8322442  0.8975273    0.9329944  0.8645377  1.0000000
## acceleration  0.4233285 -0.5046834   -0.5438005 -0.6891955 -0.4168392
## year          0.5805410 -0.3456474   -0.3698552 -0.4163615 -0.3091199
## origin        0.5652088 -0.5689316   -0.6145351 -0.4551715 -0.5850054
## mpg01         0.8369392 -0.7591939   -0.7534766 -0.6670526 -0.7577566
##              acceleration       year     origin      mpg01
## mpg             0.4233285  0.5805410  0.5652088  0.8369392
## cylinders      -0.5046834 -0.3456474 -0.5689316 -0.7591939
## displacement   -0.5438005 -0.3698552 -0.6145351 -0.7534766
## horsepower     -0.6891955 -0.4163615 -0.4551715 -0.6670526
## weight         -0.4168392 -0.3091199 -0.5850054 -0.7577566
## acceleration    1.0000000  0.2903161  0.2127458  0.3468215
## year            0.2903161  1.0000000  0.1815277  0.4299042
## origin          0.2127458  0.1815277  1.0000000  0.5136984
## mpg01           0.3468215  0.4299042  0.5136984  1.0000000

Question 14(b)

ggplot(Auto, aes(x = weight, y = mpg, color = as.factor(mpg01))) +
  geom_point() +
  theme_minimal() +
  labs(title = "Scatterplot of MPG vs Weight", x = "Weight", y = "MPG", color = "MPG01 Category")

ggplot(Auto, aes(x = as.factor(cylinders), y = mpg, fill = as.factor(mpg01))) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Boxplot of MPG by Cylinder Count and MPG01 Category", x = "Cylinders", y = "MPG", fill = "MPG01 Category")

pairs(Auto)

According to my results:

1. As weight increases, MPG tends to decrease, there is a linear relationship

2. Cars with fewer cylinders have high MPG

Question 14(c)

set.seed(123)

# size of the data set
data_size <- nrow(Auto)

# size of the training set 
train_size <- floor(0.7 * data_size)


train_indices <- sample(seq_len(data_size), size = train_size)

# Creating the training and test sets
train_set <- Auto[train_indices, ]
test_set <- Auto[-train_indices, ]

# number of observations in each set
nrow(train_set)

## [1] 274

nrow(test_set)

## [1] 118

1. 274 observations in the training set
2. 118 observations in the test set

Question 14(d)

# using lda model 

lda_model <- lda(mpg01 ~ weight + cylinders + displacement + horsepower, data = train_set)

# Print the summary of the fitted LDA model
summary(lda_model)

##         Length Class  Mode     
## prior   2      -none- numeric  
## counts  2      -none- numeric  
## means   8      -none- numeric  
## scaling 4      -none- numeric  
## lev     2      -none- character
## svd     1      -none- numeric  
## N       1      -none- numeric  
## call    3      -none- call     
## terms   3      terms  call     
## xlevels 0      -none- list

# For detailed output, including class probabilities
lda_model

## Call:
## lda(mpg01 ~ weight + cylinders + displacement + horsepower, data = train_set)
## 
## Prior probabilities of groups:
##         0         1 
## 0.4963504 0.5036496 
## 
## Group means:
##     weight cylinders displacement horsepower
## 0 3641.022  6.786765     275.2941  130.96324
## 1 2314.000  4.188406     114.5290   78.00725
## 
## Coefficients of linear discriminants:
##                        LD1
## weight       -0.0009670704
## cylinders    -0.3974647924
## displacement -0.0029615583
## horsepower    0.0049004106

LDA

lda_pred <- predict(lda_model, newdata = test_set)

test_set$predicted_lda <- lda_pred$class

# confusion matrix
confusion_matrix_lda <- table(Predicted = test_set$predicted_lda, Actual = test_set$mpg01)

# number of correct predictions
correct_predictions <- sum(diag(confusion_matrix_lda))

# Calculate the test error rate
test_error_rate_lda <- 1 - (correct_predictions / nrow(test_set))

print(test_error_rate_lda)

## [1] 0.1101695

My test error is around 11.02% which means my LDA model predicted the mpg01 status for about 88.98% of the cars in the test data, which is very good.

Question 14(e)

QDA

qda_model <- qda(mpg01 ~ weight + cylinders + displacement + horsepower, data = train_set)

qda_pred <- predict(qda_model, newdata = test_set)


test_set$predicted_qda <- qda_pred$class

# confusion matrix
confusion_matrix_qda <- table(Predicted = test_set$predicted_qda, Actual = test_set$mpg01)

# Calculate the number of correct predictions
correct_predictions_qda <- sum(diag(confusion_matrix_qda))

# Calculate the test error rate for QDA
test_error_rate_qda <- 1 - (correct_predictions_qda / nrow(test_set))

print(test_error_rate_qda)

## [1] 0.1016949

My test error is around 10.17% which means my QDA model predicted whether a car has high or low gas mileage for about 98.83 of the cars in the test set, which is very good.

Question 14(f)

# Fitting logistic regression model on the training data
logistic_model <- glm(mpg01 ~ weight + cylinders + displacement + horsepower, data = train_set, family = binomial)

# Predict on the test data
logistic_pred_probs <- predict(logistic_model, newdata = test_set, type = "response")
logistic_pred_class <- ifelse(logistic_pred_probs > 0.5, 1, 0)

# confusion matrix
confusion_matrix_logistic <- table(Predicted = logistic_pred_class, Actual = test_set$mpg01)

# calculate the number of correct predictions
correct_predictions_logistic <- sum(diag(confusion_matrix_logistic))

# calculate the test error rate for the logistic regression model
test_error_rate_logistic <- 1 - (correct_predictions_logistic / nrow(test_set))


print(test_error_rate_logistic)

## [1] 0.1101695

My test error is around 11.02%, which means my logistic regression model correctly predicted about 88.98 of the cars in the test set, which is really good.

Question 14(g)

# Fit the Naive Bayes model on the training set
nb_model <- naiveBayes(mpg01 ~ weight + cylinders + displacement + horsepower, data = train_set)

# Predict on the test set
nb_predictions <- predict(nb_model, newdata = test_set)

# confusion matrix
confusion_matrix_nb <- table(Predicted = nb_predictions, Actual = test_set$mpg01)

# calculate the number of correct predictions
correct_predictions_nb <- sum(diag(confusion_matrix_nb))

# calculate the test error rate for the Naive Bayes model
test_error_rate_nb <- 1 - (correct_predictions_nb / nrow(test_set))

print(test_error_rate_nb)

## [1] 0.1016949

my test error rate for Naive Bayes model is around 10.17%, which means my model predicted about 89.83 of the cars in the test set.

Question 14(h)

# training and test sets for KNN
# Normalize the continuous variables to ensure fair distance calculation
normalize <- function(x) {
    return ((x - min(x)) / (max(x) - min(x)))
}

train_norm <- as.data.frame(lapply(train_set[, c("weight", "cylinders", "displacement", "horsepower")], normalize))
test_norm <- as.data.frame(lapply(test_set[, c("weight", "cylinders", "displacement", "horsepower")], normalize))

# Adding the response variable back to the normalized data frames
train_norm$mpg01 <- train_set$mpg01
test_norm$mpg01 <- test_set$mpg01

# Prepare a vector to store test errors
k_values <- 1:20
test_errors <- numeric(length(k_values))

# Loop over several values of K
for (k in k_values) {
    knn_pred <- knn(train = train_norm[, -5], test = test_norm[, -5], cl = train_norm[, 5], k = k)
    confusion_matrix <- table(Predicted = knn_pred, Actual = test_norm$mpg01)
    correct_predictions <- sum(diag(confusion_matrix))
    test_errors[k] <- 1 - (correct_predictions / nrow(test_norm))
}

# Find the K value with the lowest test error
best_k <- which.min(test_errors)
best_test_error <- test_errors[best_k]

# Print the best K and its test error
cat("Best K:", best_k, "\nTest Error:", best_test_error, "\n")

## Best K: 6 
## Test Error: 0.09322034

# plotting the test errors for each K value to visualize the performance
plot(k_values, test_errors, type = "b", pch = 19, xlab = "K Value", ylab = "Test Error",
     main = "KNN Test Error Across Different K Values")

my test error is around 8.47% which means my model predicted around 91.53% of the cars in the test set, which is very very good. The best a K value of 8 had the best perfromance on the test set.

Question 15

Question 15(a)

Power <- function() {
  result <- 2**3
  print(result)
}

Power()

## [1] 8

Question 15(b)

Power2 <- function(x, a) {
  result <- x^a
  print(result)
}

Power2(3, 8)

## [1] 6561

Question 15(c)

Power2(10, 3)  

## [1] 1000

Power2(8, 17)  

## [1] 2.2518e+15

Power2(13, 13) 

## [1] 3.028751e+14

Question 15(d)

Power3 <- function(x, a) {
  result <- x^a
  return(result)
}

result1 <- Power3(10, 3)
result2 <- Power3(8, 17)
result3 <- Power3(13, 13)

print(result1)

## [1] 1000

print(result2)

## [1] 2.2518e+15

print(result3)

## [1] 3.028751e+14

Question 15(e)

Power3 <- function(x, a) {
  result <- x^a
  return(result)
}


x_values <- 1:10


y_values <- sapply(x_values, Power3, a = 2)

# plot
plot(x_values, y_values, type = "o", col = "blue",
     log = "y",  # Logarithmic scale for the y-axis
     xlab = "X", ylab = "X squared",
     main = "Plot of f(x) = x^2 on a Logarithmic Y Scale")

Question 15(f)

PlotPower <- function(x, a) {
  # Compute x^a for the given range of x
  y_values <- x^a
  
  # plot
  plot(x, y_values, type = "o", col = "blue",
       xlab = "X", ylab = paste("X^", a, sep=""),
       main = paste("Plot of f(x) = x^", a, sep=""))
}

# Example 
PlotPower(1:10, 3)