Assignment 3
Farai Muzondo
2026-06-24
The Weekly Data Set
This question should be answered using the Weekly data set, which is part of the ISLR2 package. This data is similar in nature to the Smarket data from this chapter’s lab, except that it contains 1,089 weekly returns for 21 years, from the beginning of 1990 to the end of 2010.
(a) Numerical and Graphical Summaries
Produce some numerical and graphical summaries of the Weekly data. Do there appear to be any patterns?
library(ISLR2)
library(tidyverse)
library(corrplot)
data(Weekly)
head(Weekly)## Year Lag1 Lag2 Lag3 Lag4 Lag5 Volume Today Direction
## 1 1990 0.816 1.572 -3.936 -0.229 -3.484 0.1549760 -0.270 Down
## 2 1990 -0.270 0.816 1.572 -3.936 -0.229 0.1485740 -2.576 Down
## 3 1990 -2.576 -0.270 0.816 1.572 -3.936 0.1598375 3.514 Up
## 4 1990 3.514 -2.576 -0.270 0.816 1.572 0.1616300 0.712 Up
## 5 1990 0.712 3.514 -2.576 -0.270 0.816 0.1537280 1.178 Up
## 6 1990 1.178 0.712 3.514 -2.576 -0.270 0.1544440 -1.372 Downsummary(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
##
##
##
## ggplot(Weekly, aes(x = Today)) +
geom_histogram(bins = 30, fill = "steelblue", color = "white") +
theme_minimal() +
labs(title = "Distribution of Weekly Returns", x = "Today's Return (%)", y = "Count")
ggplot(Weekly, aes(x = Volume)) +
geom_density(fill = "purple", alpha = 0.4) +
theme_minimal() +
labs(title = "Density Plot of Trading Volume", x = "Volume (in billions)", y = "Density")
cor_matrix <- cor(Weekly[, -9])
corrplot(cor_matrix, method = "color", type = "upper", tl.col = "black", tl.srt = 45)
Answer: A few clear patterns emerge from the visual summaries. First, the trading volume exhibits a right-skewed distribution but generally maintains a bell-like shape. Additionally, looking at the correlation matrix, the lag variables show virtually zero correlation with today’s returns, indicating they lack a strong linear relationship.
(b) Logistic Regression
Use the full Weekly data set to fit a logistic regression model with Direction as the response and Lag1–Lag5 along with Volume as predictors. Print the model output using summary(). Based on the regression results, determine whether any predictors show statistically significant relationships with the response. If so, identify which variables are significant.
log_model <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume,
data = Weekly, family = binomial)
summary(log_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: 4Answer: Looking at the p-values in the model
summary, Lag2 is the only variable that demonstrates
statistical significance (p < 0.05).
(c) Confusion Matrix
Compute the confusion matrix and the overall accuracy of the logistic regression model’s predictions. Then interpret the confusion matrix by describing the types of errors the model makes—for example, whether it more frequently predicts “Up” when the true direction is “Down,” or vice versa—and explain what these mistakes reveal about the model’s performance.
library(caret)
train <- Weekly[Weekly$Year < 2009, ]
test <- Weekly[Weekly$Year >= 2009, ]
prob <- predict(log_model, newdata = test, type = "response")
pred_class <- ifelse(prob > 0.5, "Up", "Down")
pred_class <- factor(pred_class, levels = c("Down", "Up"))
actual_class <- factor(test$Direction, levels = c("Down", "Up"))
table(Predicted = pred_class, Actual = actual_class)## Actual
## Predicted Down Up
## Down 17 13
## Up 26 48Answer: The confusion matrix reveals that the logistic regression model makes considerably more false positive errors than false negative ones. Specifically, it incorrectly guessed “Up” on 26 actual “Down” days, while only missing 13 actual “Up” days. Overall, the model is noticeably biased toward predicting positive market days.
(d) Logistic Regression (1990-2008)
Fit a logistic regression model using the training data from 1990–2008, with Lag2 as the sole predictor of Direction. Then evaluate the model on the held‑out observations from 2009 and 2010 by computing the confusion matrix and the overall classification accuracy. Finally, interpret the confusion matrix by describing how well the model predicts the market direction and what types of misclassifications it tends to make.
train <- (Weekly$Year < 2009)
test <- Weekly[Weekly$Year >= 2009, ]
log_model1 <- glm(Direction ~ Lag2, data = Weekly, family = binomial)
prob1 <- predict(log_model1, newdata = test, type = "response")
pred_class <- ifelse(prob1 > 0.5, "Up", "Down")
pred_class <- factor(pred_class, levels = c("Down", "Up"))
actual_class <- factor(test$Direction, levels = c("Down", "Up"))
table(Predicted = pred_class, Actual = actual_class)## Actual
## Predicted Down Up
## Down 9 5
## Up 34 56(e) LDA
Repeat (d) using LDA.
library(MASS)
Direction.test <- Weekly$Direction[!train]
lda_model <- lda(Direction ~ Lag2, data = Weekly, subset = train)
lda_pred <- predict(lda_model, test)
lda.class <- lda_pred$class
table(Predicted = lda.class, Actual = Direction.test)## Actual
## Predicted Down Up
## Down 9 5
## Up 34 56(f) QDA
Repeat (d) using QDA.
qda_model <- qda(Direction ~ Lag2, data = Weekly, subset = train)
qda_pred <- predict(qda_model, test)
qda.class <- qda_pred$class
table(Predicted = qda.class, Actual = Direction.test)## Actual
## Predicted Down Up
## Down 0 0
## Up 43 61(g) KNN (K = 1)
Repeat (d) using KNN with K = 1.
library(class)
train.index <- (Weekly$Year < 2009)
X.train <- Weekly[train.index, "Lag2", drop = FALSE]
X.test <- Weekly[!train.index, "Lag2", drop = FALSE]
Y.train <- Weekly$Direction[train.index]
Direction.test <- Weekly$Direction[!train.index]
set.seed(1)
knn.pred <- knn(train = X.train, test = X.test, cl = Y.train, k = 1)
table(Predicted = knn.pred, Actual = Direction.test)## Actual
## Predicted Down Up
## Down 21 30
## Up 22 31(h) Naive Bayes
Repeat (d) using naive Bayes.
library(e1071)
train <- (Weekly$Year < 2009)
Weekly.test <- Weekly[!train, ]
Direction.test <- Weekly$Direction[!train]
nb.fit <- naiveBayes(Direction ~ Lag2, data = Weekly, subset = train)
nb.class <- predict(nb.fit, Weekly.test)
table(Predicted = nb.class, Actual = Direction.test)## Actual
## Predicted Down Up
## Down 0 0
## Up 43 61Which of these methods appears to provide the best results on this data? Both the logistic regression and the Linear Discriminant Analysis (LDA) models yield the most accurate predictions, tying for the best overall performance.
(i) Experimentation
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.
# Logistic Regression with Lag1 + Lag2
best.logit <- glm(Direction ~ Lag1 + Lag2,
data = Weekly,
family = binomial,
subset = train)
best.probs <- predict(best.logit, Weekly.test, type = "response")
best.pred <- rep("Down", length(best.probs))
best.pred[best.probs > 0.5] <- "Up"
table(Predicted = best.pred, Actual = Direction.test)## Actual
## Predicted Down Up
## Down 7 8
## Up 36 53# KNN with K = 9
X.train <- Weekly[train, "Lag2", drop = FALSE]
X.test <- Weekly[!train, "Lag2", drop = FALSE]
Y.train <- Weekly$Direction[train]
set.seed(1)
knn.pred.9 <- knn(train = X.train, test = X.test, cl = Y.train, k = 9)
table(Predicted = knn.pred.9, Actual = Direction.test)## Actual
## Predicted Down Up
## Down 17 20
## Up 26 41The Auto Data Set
In this problem, you will develop a model to predict whether a given car gets high or low gas mileage based on the Auto data set.
(a) Create Binary Variable mpg01
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.
mpg_median <- median(Auto$mpg)
mpg01 <- as.integer(Auto$mpg > mpg_median)
Auto_new <- data.frame(mpg01, Auto)
head(Auto_new)## mpg01 mpg cylinders displacement horsepower weight acceleration year origin
## 1 0 18 8 307 130 3504 12.0 70 1
## 2 0 15 8 350 165 3693 11.5 70 1
## 3 0 18 8 318 150 3436 11.0 70 1
## 4 0 16 8 304 150 3433 12.0 70 1
## 5 0 17 8 302 140 3449 10.5 70 1
## 6 0 15 8 429 198 4341 10.0 70 1
## name
## 1 chevrolet chevelle malibu
## 2 buick skylark 320
## 3 plymouth satellite
## 4 amc rebel sst
## 5 ford torino
## 6 ford galaxie 500(b) Graphical Exploration
Explore the data graphically in order to investigate the association between mpg01 and the other features.
boxplot(cylinders ~ mpg01, data = Auto_new, main = "Cylinders vs mpg01", xlab = "mpg01", ylab = "Cylinders", col = "lightblue")
boxplot(displacement ~ mpg01, data = Auto_new, main = "Displacement vs mpg01", xlab = "mpg01", ylab = "Displacement", col = "lightgreen")
boxplot(horsepower ~ mpg01, data = Auto_new, main = "Horsepower vs mpg01", xlab = "mpg01", ylab = "Horsepower", col = "lightpink")
boxplot(weight ~ mpg01, data = Auto_new, main = "Weight vs mpg01", xlab = "mpg01", ylab = "Weight", col = "lightyellow")
boxplot(acceleration ~ mpg01, data = Auto_new, main = "Acceleration vs mpg01", xlab = "mpg01", ylab = "Acceleration", col = "lavender")
boxplot(year ~ mpg01, data = Auto_new, main = "Year vs mpg01", xlab = "mpg01", ylab = "Year", col = "wheat")
pairs(Auto_new[, c("mpg01", "cylinders", "displacement", "horsepower", "weight", "acceleration", "year")])
(c) Train/Test Split
Split the data into a training set and a test set.
set.seed(42)
train_size <- floor(0.80 * nrow(Auto_new))
train_indices <- sample(seq_len(nrow(Auto_new)), size = train_size)
train_set <- Auto_new[train_indices, ]
test_set <- Auto_new[-train_indices, ](d) LDA Model
Perform Linear Discriminant Analysis (LDA) on the training data to predict mpg01, using only the variables identified in part (b) as most strongly associated with the response. Evaluate the fitted model on the test set and report the resulting test error rate
lda_automodel <- lda(mpg01 ~ weight + displacement + horsepower + cylinders, data = train_set)
print(lda_automodel)## Call:
## lda(mpg01 ~ weight + displacement + horsepower + cylinders, data = train_set)
##
## Prior probabilities of groups:
## 0 1
## 0.5271565 0.4728435
##
## Group means:
## weight displacement horsepower cylinders
## 0 3597.970 272.5394 129.38788 6.751515
## 1 2345.905 117.0912 78.66892 4.202703
##
## Coefficients of linear discriminants:
## LD1
## weight -0.0009597286
## displacement -0.0016031063
## horsepower 0.0013528155
## cylinders -0.3798340772lda_autopred <- predict(lda_automodel, newdata = test_set)
lda_autoclass <- lda_autopred$class
test_error <- mean(lda_autoclass != test_set$mpg01)
print(paste("LDA Test Error:", test_error))## [1] "LDA Test Error: 0.0632911392405063"Answer: The calculated test error for the LDA model comes out to roughly 6.33%.
(e) QDA Model
Perform QDA on the training data…
qda_automodel <- qda(mpg01 ~ weight + displacement + horsepower + cylinders, data = train_set)
print(qda_automodel)## Call:
## qda(mpg01 ~ weight + displacement + horsepower + cylinders, data = train_set)
##
## Prior probabilities of groups:
## 0 1
## 0.5271565 0.4728435
##
## Group means:
## weight displacement horsepower cylinders
## 0 3597.970 272.5394 129.38788 6.751515
## 1 2345.905 117.0912 78.66892 4.202703qda_autopred <- predict(qda_automodel, newdata = test_set)
qda_autoclass <- qda_autopred$class
test_error_qda <- mean(qda_autoclass != test_set$mpg01)
print(paste("QDA Test Error:", test_error_qda))## [1] "QDA Test Error: 0.0759493670886076"Answer: For the QDA model, the test error increases slightly to approximately 7.60%.
(f) Logistic Regression
Perform logistic regression on the training data…
log_automodel <- glm(mpg01 ~ weight + displacement + horsepower + cylinders,
data = train_set, family = binomial)
summary(log_automodel)##
## Call:
## glm(formula = mpg01 ~ weight + displacement + horsepower + cylinders,
## family = binomial, data = train_set)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 11.4938653 1.8908145 6.079 1.21e-09 ***
## weight -0.0018412 0.0007409 -2.485 0.0130 *
## displacement -0.0109249 0.0088261 -1.238 0.2158
## horsepower -0.0471302 0.0148478 -3.174 0.0015 **
## cylinders -0.0045738 0.3720303 -0.012 0.9902
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 432.99 on 312 degrees of freedom
## Residual deviance: 174.26 on 308 degrees of freedom
## AIC: 184.26
##
## Number of Fisher Scoring iterations: 7prob_autolog <- predict(log_automodel, newdata = test_set, type = "response")
pred_autolog <- ifelse(prob_autolog > 0.5, 1, 0)
test_error_log <- mean(pred_autolog != test_set$mpg01)
print(paste("Logistic Test Error:", test_error_log))## [1] "Logistic Test Error: 0.0632911392405063"(g) Naive Bayes
Perform naive Bayes on the training data…
nb_automodel <- naiveBayes(mpg01 ~ weight + displacement + horsepower + cylinders, data = train_set)
print(nb_automodel)##
## Naive Bayes Classifier for Discrete Predictors
##
## Call:
## naiveBayes.default(x = X, y = Y, laplace = laplace)
##
## A-priori probabilities:
## Y
## 0 1
## 0.5271565 0.4728435
##
## Conditional probabilities:
## weight
## Y [,1] [,2]
## 0 3597.970 670.9549
## 1 2345.905 409.4496
##
## displacement
## Y [,1] [,2]
## 0 272.5394 89.45761
## 1 117.0912 41.22017
##
## horsepower
## Y [,1] [,2]
## 0 129.38788 36.80148
## 1 78.66892 16.11328
##
## cylinders
## Y [,1] [,2]
## 0 6.751515 1.4542036
## 1 4.202703 0.7186482nb_autopred <- predict(nb_automodel, newdata = test_set)
test_error_nb <- mean(nb_autopred != test_set$mpg01)
print(paste("Naive Bayes Test Error:", test_error_nb))## [1] "Naive Bayes Test Error: 0.0632911392405063"(h) KNN Optimization
Perform KNN on the training data, with several values of K…
predictors <- c("weight", "displacement", "horsepower", "cylinders")
scaled_predictors <- scale(Auto_new[, predictors])
train_X <- scaled_predictors[train_indices, ]
test_X <- scaled_predictors[-train_indices, ]
train_y <- Auto_new$mpg01[train_indices]
test_y <- Auto_new$mpg01[-train_indices]
k_values <- c(1, 3, 5, 7, 10, 15, 20, 50)
knn_errors <- numeric(length(k_values))
for (i in seq_along(k_values)) {
set.seed(42)
knn_pred <- knn(train = train_X, test = test_X, cl = train_y, k = k_values[i])
knn_errors[i] <- mean(knn_pred != test_y)
}
knn_results <- data.frame(K = k_values, Test_Error = round(knn_errors, 4))
print(knn_results)## K Test_Error
## 1 1 0.1266
## 2 3 0.0886
## 3 5 0.0759
## 4 7 0.0633
## 5 10 0.0506
## 6 15 0.0633
## 7 20 0.0633
## 8 50 0.0633Answer: Based on the testing loop, K = 10 yields the lowest error rate, making it the optimal choice. The values of K = 15, 20, and 50 also perform comparably well.
The Boston Data Set
Using the Boston data set, fit classification models in order to predict whether a given census tract has a crime rate above or below the median.
crime_median <- median(Boston$crim)
crim01 <- as.integer(Boston$crim > crime_median)
Boston_new <- data.frame(crim01, Boston)
Boston_new$crim <- NULL
cor(Boston_new)[, "crim01"]## crim01 zn indus chas nox rm
## 1.00000000 -0.43615103 0.60326017 0.07009677 0.72323480 -0.15637178
## age dis rad tax ptratio black
## 0.61393992 -0.61634164 0.61978625 0.60874128 0.25356836 -0.35121093
## lstat medv
## 0.45326273 -0.26301673boxplot(rad ~ crim01, data = Boston_new, main = "rad vs crim01", xlab = "crim01", ylab = "rad", col = "lightblue")
boxplot(indus ~ crim01, data = Boston_new, main = "indus vs crim01", xlab = "crim01", ylab = "indus", col = "lightgreen")
boxplot(nox ~ crim01, data = Boston_new, main = "nox vs crim01", xlab = "crim01", ylab = "nox", col = "lightpink")
boxplot(age ~ crim01, data = Boston_new, main = "age vs crim01", xlab = "crim01", ylab = "age", col = "lightyellow")
boxplot(dis ~ crim01, data = Boston_new, main = "dis vs crim01", xlab = "crim01", ylab = "dis", col = "lavender")
boxplot(tax ~ crim01, data = Boston_new, main = "Year vs crim01", xlab = "crim01", ylab = "tax", col = "wheat")
set.seed(123)
train_indices <- sample(1:nrow(Boston_new), 0.8 * nrow(Boston_new))
bostontrain_set <- Boston_new[train_indices, ]
bostontest_set <- Boston_new[-train_indices, ]
# 1. Logistic Regression
boston_log <- glm(crim01 ~ nox + rad + tax + dis, data = bostontrain_set, family = binomial)
summary(boston_log)##
## Call:
## glm(formula = crim01 ~ nox + rad + tax + dis, family = binomial,
## data = bostontrain_set)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.432095 3.566290 -6.010 1.86e-09 ***
## nox 39.659745 6.067602 6.536 6.31e-11 ***
## rad 0.552769 0.123378 4.480 7.45e-06 ***
## tax -0.009133 0.002500 -3.654 0.000258 ***
## dis 0.077851 0.156504 0.497 0.618880
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 560.02 on 403 degrees of freedom
## Residual deviance: 202.12 on 399 degrees of freedom
## AIC: 212.12
##
## Number of Fisher Scoring iterations: 8bostonlog_prob <- predict(boston_log, newdata = bostontest_set, type = "response")
bostonlog_pred <- ifelse(bostonlog_prob > 0.5, 1, 0)
bostonlog_error <- mean(bostonlog_pred != bostontest_set$crim01)
print(paste("Logistic Error:", bostonlog_error))## [1] "Logistic Error: 0.127450980392157"# 2. LDA
boston_lda <- lda(crim01 ~ nox + rad + tax + dis, data = bostontrain_set)
bostonpred_lda <- predict(boston_lda, newdata = bostontest_set)$class
lda_bostonerror <- mean(bostonpred_lda != bostontest_set$crim01)
print(paste("LDA Error:", lda_bostonerror))## [1] "LDA Error: 0.156862745098039"# 3. Naive Bayes
boston_nb <- naiveBayes(crim01 ~ nox + rad + tax + dis, data = bostontrain_set)
boston_nbpred <- predict(boston_nb, newdata = bostontest_set)
boston_nberror <- mean(boston_nbpred != bostontest_set$crim01)
print(paste("Naive Bayes Error:", boston_nberror))## [1] "Naive Bayes Error: 0.156862745098039"# 4. KNN
bostonscaled_features <- scale(Boston_new[, -1])
train_boston <- bostonscaled_features[train_indices, c("nox", "rad", "tax", "dis")]
test_boston <- bostonscaled_features[-train_indices, c("nox", "rad", "tax", "dis")]
train_y <- Boston_new$crim01[train_indices]
test_y <- Boston_new$crim01[-train_indices]
set.seed(123)
knn_pred <- knn(train_boston, test_boston, train_y, k = 5)
err_knn <- mean(knn_pred != test_y)
print(paste("KNN Error:", err_knn))## [1] "KNN Error: 0.0686274509803922"