Assignment #3

Question 13

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.

library(ISLR2)
data(Weekly)

(a)

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

summary(Weekly)
cor(Weekly[, -9])
table(Weekly$Direction)

pairs(Weekly)

plot(Weekly$Volume, type = "l", xlab = "Week (1990-2010)", ylab = "Volume")

Overall, thee is an upward trend in the trading volume over the years, which is not indicative of a linear relationship. This is explained as there is a positive correlation between volume and year of 0.84.

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

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

summary(glm.fit)

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

The only significant predictor is Lag2, with a p value is 0.0296 < 0.05. There is a slight higher chance that there is an upward trend if there was a larger return two weeks ago.

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

glm.probs <- predict(glm.fit, type = "response")
glm.pred  <- rep("Down", nrow(Weekly))
glm.pred[glm.probs > 0.5] <- "Up"

table(glm.pred, Weekly$Direction)

        
glm.pred Down  Up
    Down   54  48
    Up    430 557

mean(glm.pred == Weekly$Direction) 

[1] 0.5610652

The confusion matrix shows that the model predicts Up for 987 of the 1089 weeks. However, the model has an accuracy of 56.1% since the logictic regression correctly classifies 611 of the 1089 weeks where (54+557)/1089 = 611 as there are fall Up predictions.

(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 <= 2008
Weekly.test <- Weekly[!train, ]
Direction.test <- Weekly$Direction[!train]

# lag2 as only predictor
glm.fit <- glm(Direction ~ Lag2, data = Weekly, family = binomial,
               subset = train)

# data from 2009 and 2010
glm.probs <- predict(glm.fit, Weekly.test, type = "response")
glm.pred  <- rep("Down", nrow(Weekly.test))
glm.pred[glm.probs > 0.5] <- "Up"

table(glm.pred, Direction.test)

        Direction.test
glm.pred Down Up
    Down    9  5
    Up     34 56

mean(glm.pred == Direction.test)

[1] 0.625

(e)

Repeat (d) using LDA.

library(MASS)

lda.fit <- lda(Direction ~ Lag2, data = Weekly, subset = train)

lda.pred  <- predict(lda.fit, Weekly.test)
lda.class <- lda.pred$class

table(lda.class, Direction.test)

         Direction.test
lda.class Down Up
     Down    9  5
     Up     34 56

mean(lda.class == Direction.test)

[1] 0.625

(f)

Repeat (d) using QDA.

qda.fit <- qda(Direction ~ Lag2, data = Weekly, subset = train)

qda.class <- predict(qda.fit, Weekly.test)$class

table(qda.class, Direction.test)

         Direction.test
qda.class Down Up
     Down    0  0
     Up     43 61

mean(qda.class == Direction.test)

[1] 0.5865385

(g)

Repeat (d) using KNN with K = 1.

library(class)
train.X <- as.matrix(Weekly$Lag2[train])
test.X  <- as.matrix(Weekly$Lag2[!train])
train.Direction <- Weekly$Direction[train]

set.seed(240)
knn.pred <- knn(train.X, test.X, train.Direction, k = 1)

table(knn.pred, Direction.test)

        Direction.test
knn.pred Down Up
    Down   21 30
    Up     22 31

mean(knn.pred == Direction.test)

[1] 0.5

(h)

Repeat (d) using naive Bayes.

library(e1071)

nb.fit <- naiveBayes(Direction ~ Lag2, data = Weekly, subset = train)

nb.class <- predict(nb.fit, Weekly.test)

table(nb.class, Direction.test)

        Direction.test
nb.class Down Up
    Down    0  0
    Up     43 61

mean(nb.class == Direction.test)

[1] 0.5865385

(i)

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

Out of the five methods, logistic regression and LDA have the best results, both at 62.5% since Lag2 as a predictor is only slightly significant and contributes to bias-variance.

(j)

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.

train <- Weekly$Year <= 2008

X <- cbind(Lag2 = Weekly$Lag2, Lag1Lag2 = Weekly$Lag1 * Weekly$Lag2)
X <- scale(X)

train.X <- X[train, ]
test.X  <- X[!train, ]
train.Direction <- Weekly$Direction[train]
Direction.test  <- Weekly$Direction[!train]

set.seed(240)
knn.pred <- knn(train.X, test.X, train.Direction, k = 100)
table(knn.pred, Direction.test)

        Direction.test
knn.pred Down Up
    Down   16 11
    Up     27 50

mean(knn.pred == Direction.test)

[1] 0.6346154

By using Lag2 and the interaction between Lag1 and Lag2, there is very minimal difference between 62.5% (logistic regression using only Lag2) and 63.5% with only an additional week when using KNN. The predictor Lag2 is too weak to vastly improve the model.

Question 14

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 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.

#library(ISLR2)
data(Auto)

# prints median
med_mpg <- median(Auto$mpg)
print(med_mpg)

[1] 22.75

# assigns 1 or 0 if over or under median
mpg01 <- ifelse(Auto$mpg > median(Auto$mpg), 1, 0)
Auto2 <- data.frame(mpg01 = factor(mpg01), Auto)

head(Auto2)

table(mpg01)

mpg01
  0   1 
196 196 

The result is balanced where there are 196 variables that are over the median of 22.75 mpg and 196 that are under 22.75 mpg.

(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, 3))                       # 2x3 grid of plots
boxplot(displacement ~ mpg01, data = Auto2, xlab = "mpg01", ylab = "displacement")
boxplot(horsepower   ~ mpg01, data = Auto2, xlab = "mpg01", ylab = "horsepower")

boxplot(weight       ~ mpg01, data = Auto2, xlab = "mpg01", ylab = "weight")
boxplot(acceleration ~ mpg01, data = Auto2, xlab = "mpg01", ylab = "acceleration")

boxplot(year         ~ mpg01, data = Auto2, xlab = "mpg01", ylab = "year")
boxplot(cylinders    ~ mpg01, data = Auto2, xlab = "mpg01", ylab = "cylinders")

par(mfrow = c(1, 1))    

pairs(Auto2[, c("displacement","horsepower","weight","acceleration")], 
      col = ifelse(mpg01 == 1, "blue", "red"))

cor(Auto2[, c("mpg01","cylinders","displacement","horsepower","weight", "acceleration","year","origin")])[,"mpg01"]

       mpg01    cylinders displacement   horsepower       weight acceleration 
   1.0000000   -0.7591939   -0.7534766   -0.6670526   -0.7577566    0.3468215 
        year       origin 
   0.4299042    0.5136984 

Other factors, such as cylinders, displacement, horsepower, and weight have a relationship with mpg01. There is a negative correlation from -0.67 and -0.76, indicating that cars that are lighter, have fewer cylinders, and horsepower, have higher mileage.

(c)

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

mpg01 <- ifelse(Auto$mpg > median(Auto$mpg), 1, 0)
Auto2 <- data.frame(mpg01 = factor(mpg01), Auto)

set.seed(240)
n <- nrow(Auto2)
train <- sample(n, size = 0.7 * n)

Auto.train <- Auto2[train, ]
Auto.test  <- Auto2[-train, ]
mpg01.test <- Auto2$mpg01[-train]

dim(Auto.train)

[1] 274  10

dim(Auto.test)

[1] 118  10

274 + 118 = 392, which is the full size of the Auto dataset.

274/392 = 0.699 so this is a 70/30 split.

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

#library(MASS)

# predictors from mpg01 in part (b)
lda.fit <- lda(mpg01 ~ cylinders + displacement + horsepower + weight,
               data = Auto.train)

lda.pred  <- predict(lda.fit, Auto.test)
lda.class <- lda.pred$class

table(lda.class, mpg01.test)

         mpg01.test
lda.class  0  1
        0 44  4
        1 10 60

mean(lda.class != mpg01.test)

[1] 0.1186441

There is a test error of 11% using the predictors of cylinders, displacement, horsepower, and weight, meaning that there is around 90% accuracy when applying a fitted model. Engine and size variables offer a clean separation of the high and low mileage cars.

(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 ~ cylinders + displacement + horsepower + weight, data = Auto.train)

qda.class <- predict(qda.fit, Auto.test)$class

table(qda.class, mpg01.test)

         mpg01.test
qda.class  0  1
        0 46  6
        1  8 58

mean(qda.class != mpg01.test)

[1] 0.1186441

QDA performed the same as LDA with an error of 11%, meaning that there is a near 90% accuracy.

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

glm.fit <- glm(mpg01 ~ cylinders + displacement + horsepower + weight, data = Auto.train, family = binomial)

glm.probs <- predict(glm.fit, Auto.test, type = "response")
glm.pred  <- rep(0, nrow(Auto.test))
glm.pred[glm.probs > 0.5] <- 1

table(glm.pred, mpg01.test)

        mpg01.test
glm.pred  0  1
       0 45  6
       1  9 58

mean(glm.pred != mpg01.test)

[1] 0.1271186

The test error is now at 12% compared to the 11% obtained from LDA and QDA, so there is now 89% accuracy. But overall, all three methods are almost identical and indicates there is a clear boundary separating high and low mpg for the Auto dataset.

(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 ~ cylinders + displacement + horsepower + weight, data = Auto.train)

nb.class <- predict(nb.fit, Auto.test)


table(nb.class, mpg01.test)

        mpg01.test
nb.class  0  1
       0 45  6
       1  9 58

mean(nb.class != mpg01.test)

[1] 0.1271186

The Naive Bayes model is the same result as the logistic regression of 12% test error with 89% accuracy.

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

vars <- c("cylinders", "displacement", "horsepower", "weight")
X <- scale(Auto2[, vars])

train.X <- X[train, ]
test.X  <- X[-train, ]
train.mpg01 <- Auto2$mpg01[train]

set.seed(1)
for (k in c(1, 3, 5, 7, 10, 15, 20, 50, 100)) {
  knn.pred <- knn(train.X, test.X, train.mpg01, k = k)
  cat("K =", k, " test error =", round(mean(knn.pred != mpg01.test), 4), "\n")}

K = 1  test error = 0.1441 
K = 3  test error = 0.1102 
K = 5  test error = 0.1186 
K = 7  test error = 0.1186 
K = 10  test error = 0.1186 
K = 15  test error = 0.1186 
K = 20  test error = 0.1186 
K = 50  test error = 0.1186 
K = 100  test error = 0.1271 

K = 3 performs the best as there is a test error of 11.02% and around 89% accuracy. Using a very small K of 1 over fits, while using a very large K of 100 under fits. K of 5 - 50 have the same test error, but K = 3 is more accurate.

Question 16

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. Explore logistic regression, LDA, naive Bayes, and KNN models using various subsets of the predictors. Describe your findings. Hint: You will have to create the response variable yourself, using the variables that are contained in the Boston data set.

data(Boston)
crim01 <- ifelse(Boston$crim > median(Boston$crim), 1, 0)
Boston2 <- data.frame(crim01 = factor(crim01), Boston)

# Which predictors are most associated with crim01?
cor(data.frame(crim01, Boston[, -1]))[, "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.26301673 

# train and test split
set.seed(240)
n <- nrow(Boston2)
train <- sample(n, 0.7 * n)
Boston.train <- Boston2[train, ]
Boston.test  <- Boston2[-train, ]
crim01.test  <- Boston2$crim01[-train]

# logistic regression
glm.fit  <- glm(crim01 ~ nox + rad + dis + age + tax + indus, data = Boston.train, family = binomial)
glm.prob <- predict(glm.fit, Boston.test, type = "response")
glm.pred <- ifelse(glm.prob > 0.5, 1, 0)
mean(glm.pred != crim01.test)

[1] 0.1118421

# LDA
lda.fit <- lda(crim01 ~ nox + rad + dis + age + tax + indus, data = Boston.train)
mean(predict(lda.fit, Boston.test)$class != crim01.test)

[1] 0.1578947

# Naive Bayes
nb.fit <- naiveBayes(crim01 ~ nox + rad + dis + age + tax + indus, data = Boston.train)
mean(predict(nb.fit, Boston.test) != crim01.test)

[1] 0.1776316

## KNN
vars <- c("nox","rad","dis")
X <- scale(Boston2[, vars])
set.seed(240)
for (k in c(1,3,5,10)) {
  knn.pred <- knn(X[train,], X[-train,], Boston2$crim01[train], k = k)
  cat("KNN K=",k," error=", round(mean(knn.pred != crim01.test),4), "\n")}

KNN K= 1  error= 0.0658 
KNN K= 3  error= 0.0526 
KNN K= 5  error= 0.0592 
KNN K= 10  error= 0.0658 

KNN model is a better fit as there is significant lower test error by separating the high and low crime tracts among all K values. Additionally, using fewer but stronger predictors makes a drastic difference.

Logistic regression is a better fit of LDA with a test error of 11% compared to about 16%.

Naive Bayes has the worst test error of about 18% since it assumes the predictors are independent as they are actually very correlated with one another.