Assignment 3
Raegan Gutierrez
2/26/2021
Question 10
This question should be answered using the Weekly data set, which is part of the ISLR 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. Produce some numerical and graphical summaries of the Weekly data. Do there appear to be any patterns?
library(ISLR)
attach(Weekly)
pairs(Weekly)
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.000000000It appears the only relationship is between Volume and Year. At 0.84194, there is a strong positive relationship between the two variables. This correlation can even be seen on the pairs scatter plot between the two variables. All other variables seem to have little to no correlation based both on their scatter plots and correlation values given.
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.direct= glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, data = Weekly, family = "binomial")
summary(glm.direct)##
## Call:
## glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 +
## Volume, family = "binomial", data = Weekly)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6949 -1.2565 0.9913 1.0849 1.4579
##
## 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: 4The only predictor that appears to be statistically significant is Lag 2.
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.
probs= predict(glm.direct, type = "response")
glm.pred = rep("Down", 1089)
glm.pred[probs > 0.5] = "Up"
table(glm.pred, Direction)## Direction
## glm.pred Down Up
## Down 54 48
## Up 430 557(54+557)/1089## [1] 0.5610652557/(48+557)## [1] 0.920661254/(54+430)## [1] 0.1115702Based on the confusion matrix, the model’s predictions were correct about 56.11% of the time, leaving a test error rate of 43.89%. The model accurately guessed the stock market going up about 92.07% of the time, but only accurately predicted the market going down for the week at 11.16% of the time.
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 = (Year < 2009)
Weekly.pred = Weekly[!train, ]
Direction.pred = Direction[!train]
glm.2 = glm(Direction ~ Lag2, data = Weekly, family = binomial, subset = train)
summary(glm.2)##
## Call:
## glm(formula = Direction ~ Lag2, family = binomial, data = Weekly,
## subset = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.536 -1.264 1.021 1.091 1.368
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.20326 0.06428 3.162 0.00157 **
## Lag2 0.05810 0.02870 2.024 0.04298 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1354.7 on 984 degrees of freedom
## Residual deviance: 1350.5 on 983 degrees of freedom
## AIC: 1354.5
##
## Number of Fisher Scoring iterations: 4probs2 = predict(glm.2, Weekly.pred, type = "response")
glm.2 = rep("Down", length(probs2))
glm.2[probs2 > 0.5] = "Up"
table(glm.2, Direction.pred)## Direction.pred
## glm.2 Down Up
## Down 9 5
## Up 34 56(9+56)/104## [1] 0.62556/(56+5)## [1] 0.91803289/(9+34)## [1] 0.2093023This model was more accurate, making correct predictions 62.5% of the time, leaving an error rate of 37.5%. The model accurately predicted the stock market going up for the week 91.80% of the time. The model still struggles to predict when the market is down for the week, with accurate estimates at 20.93% of the time.
e. Repeat (d) using LDA.
library(MASS)
lda.direct = lda(Direction ~ Lag2, data = Weekly, subset = train)
lda.direct## Call:
## lda(Direction ~ Lag2, data = Weekly, subset = train)
##
## Prior probabilities of groups:
## Down Up
## 0.4477157 0.5522843
##
## Group means:
## Lag2
## Down -0.03568254
## Up 0.26036581
##
## Coefficients of linear discriminants:
## LD1
## Lag2 0.4414162lda.pred = predict(lda.direct, Weekly.pred)
table(lda.pred$class, Direction.pred)## Direction.pred
## Down Up
## Down 9 5
## Up 34 56(9+56)/104## [1] 0.62556/(56+5)## [1] 0.91803289/(9+34)## [1] 0.2093023In this case, the LDA model got the same results from the GLM model, with an accuracy of 62.5% and a test error rate of 37.5%. Again, the model predicted the stock market going up and down for the week at 91.80% and 20.93%, respectively.
f. Repeat (d) using QDA.
qda.direct = qda(Direction ~ Lag2, data = Weekly, subset = train)
qda.direct## Call:
## qda(Direction ~ Lag2, data = Weekly, subset = train)
##
## Prior probabilities of groups:
## Down Up
## 0.4477157 0.5522843
##
## Group means:
## Lag2
## Down -0.03568254
## Up 0.26036581qda.pred = predict(qda.direct, Weekly.pred)
table(qda.pred$class, Direction.pred)## Direction.pred
## Down Up
## Down 0 0
## Up 43 6161/104## [1] 0.5865385In this case, the model accurately predicts the stock market’s movements for the week 58.65% of the time, leaving an error rate of 41.35%. However, the model only guesses “up” every time. Therefore, it guesses when the stock market increases for the week at 100% accuracy, but never guesses when the stock market decreases accurately. Even with close to 60% accuracy, it is still not a trustworthy model.
g. Repeat (d) using KNN with \(K=1\).
library(class)
train.X = as.matrix(Lag2[train])
test.X = as.matrix(Lag2[!train])
train.Direction = Direction[train]
set.seed(1)
knn.pred = knn(train.X, test.X, train.Direction, k = 1)
table(knn.pred, Direction.pred)## Direction.pred
## knn.pred Down Up
## Down 21 30
## Up 22 31(21+31)/104## [1] 0.531/61## [1] 0.508196721/43## [1] 0.4883721The KNN model predicted the stock market’s movements with 50% accuracy, leaving a test error rate of 50%. The model accurately predicted when the stock market would increase for the week 50.82% of the time, and accurately predicted when the market would go down for the week 48.84% of the time.
h. Which of these methods appears to provide the best results on this data?
The LDA and the second GLM model(part d) have the best results, both providing accurate predictions 62.5% of the time, leaving the lowest test error rate of 37.5%. A notable second was the QDA model, or better known as the “up” model, with an accuracy of 58.65% by just predicting “up” every time. Lastly, we have the KNN model at 50% accuracy, which is as accurate as randomly guessing every week.
i. 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. linear model
glm.3 = glm(Direction ~ Lag1+Lag2, data = Weekly, family = binomial, subset = train)
probs3 = predict(glm.3, Weekly.pred, type = "response")
glm3.pred = rep("Down", length(probs3))
glm3.pred[probs3 > 0.5] = "Up"
table(glm3.pred, Direction.pred)## Direction.pred
## glm3.pred Down Up
## Down 7 8
## Up 36 53mean(glm3.pred == Direction.pred)## [1] 0.5769231LDA
lda.direct2 = lda(Direction ~ Lag2+Lag1, data = Weekly, subset = train)
lda.pred2 = predict(lda.direct2, Weekly.pred)
mean(lda.pred2$class == Direction.pred)## [1] 0.5769231QDA
qda.direct2 = qda(Direction ~ Lag2 + sqrt(abs(Lag2)), data = Weekly, subset = train)
qda.pred2 = predict(qda.direct2, Weekly.pred)
table(qda.pred2$class, Direction.pred)## Direction.pred
## Down Up
## Down 12 13
## Up 31 48mean(qda.pred2 == Direction.pred)## [1] 0KKN
knn.pred2 = knn(train.X, test.X, train.Direction, k = 15)
table(knn.pred2, Direction.pred)## Direction.pred
## knn.pred2 Down Up
## Down 20 20
## Up 23 41mean(knn.pred2 == Direction.pred)## [1] 0.5865385knn.pred3 = knn(train.X, test.X, train.Direction, k = 50)
table(knn.pred3, Direction.pred)## Direction.pred
## knn.pred3 Down Up
## Down 21 20
## Up 22 41mean(knn.pred3 == Direction.pred)## [1] 0.5961538Even after all of the experimental models, the original linear regression from part d and original LDA model are the most accurate predictors at 62.5% accuracy.
Question 11
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.
attach(Auto)
mpg01 = rep(1, length(mpg))
mpg01[mpg < median(mpg)] = 0
Auto = data.frame(Auto, mpg01)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.
pairs(Auto)
boxplot(weight ~ mpg01, data = Auto, main = "Weight vs mpg01")
boxplot(displacement ~ mpg01, data = Auto, main = "Displacement vs mpg01")
boxplot(horsepower ~ mpg01, data = Auto, main = "Horsepower vs mpg01")
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.0000000The features that seem most likely to predict mpg01 based on graphical data are weight, displacement, and horsepower. All of these variables show a strong, negative correlation with mpg01 in the scatter plots, which is verified by the box plots. Looking at the correlations of the data, cylinders is also a good predictor of mpg01.
c. Split the data into a training set and a test set.
train = (year %% 2 == 0)
Auto.train = Auto[train, ]
Auto.test = Auto[!train, ]
mpg01.test = mpg01[!train]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?
lda.mpg01 = lda(mpg01 ~ cylinders + weight + displacement + horsepower, data = Auto, subset = train)
lda.mpg01## Call:
## lda(mpg01 ~ cylinders + weight + displacement + horsepower, data = Auto,
## subset = train)
##
## Prior probabilities of groups:
## 0 1
## 0.4571429 0.5428571
##
## Group means:
## cylinders weight displacement horsepower
## 0 6.812500 3604.823 271.7396 133.14583
## 1 4.070175 2314.763 111.6623 77.92105
##
## Coefficients of linear discriminants:
## LD1
## cylinders -0.6741402638
## weight -0.0011465750
## displacement 0.0004481325
## horsepower 0.0059035377lda.pred =predict(lda.mpg01, Auto.test)
table(lda.pred$class, mpg01.test)## mpg01.test
## 0 1
## 0 86 9
## 1 14 73mean(lda.pred$class != mpg01.test)## [1] 0.1263736The test error rate of the LDA model is about 12.64%.
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.mpg01 = qda(mpg01 ~ cylinders + weight + displacement + horsepower, data = Auto, subset = train)
qda.mpg01## Call:
## qda(mpg01 ~ cylinders + weight + displacement + horsepower, data = Auto,
## subset = train)
##
## Prior probabilities of groups:
## 0 1
## 0.4571429 0.5428571
##
## Group means:
## cylinders weight displacement horsepower
## 0 6.812500 3604.823 271.7396 133.14583
## 1 4.070175 2314.763 111.6623 77.92105qda.pred =predict(qda.mpg01, Auto.test)
table(qda.pred$class, mpg01.test)## mpg01.test
## 0 1
## 0 89 13
## 1 11 69mean(qda.pred$class != mpg01.test)## [1] 0.1318681The test error rate of the QDA model is about 13.19%.
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.mpg01 = glm(mpg01 ~ cylinders + weight + displacement + horsepower, data = Auto, family= binomial, subset = train)
summary(glm.mpg01)##
## Call:
## glm(formula = mpg01 ~ cylinders + weight + displacement + horsepower,
## family = binomial, data = Auto, subset = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.48027 -0.03413 0.10583 0.29634 2.57584
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 17.658730 3.409012 5.180 2.22e-07 ***
## cylinders -1.028032 0.653607 -1.573 0.1158
## weight -0.002922 0.001137 -2.569 0.0102 *
## displacement 0.002462 0.015030 0.164 0.8699
## horsepower -0.050611 0.025209 -2.008 0.0447 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 289.58 on 209 degrees of freedom
## Residual deviance: 83.24 on 205 degrees of freedom
## AIC: 93.24
##
## Number of Fisher Scoring iterations: 7probs = predict(glm.mpg01, Auto.test, type = "response")
glm.pred = rep(0, length(probs))
glm.pred[probs > 0.5] = 1
table(glm.pred, mpg01.test)## mpg01.test
## glm.pred 0 1
## 0 89 11
## 1 11 71mean(glm.pred != mpg01.test)## [1] 0.1208791The test error rate of the linear regression model is about 12.09%.
g. 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?
\(K=1\)
train.X = cbind(cylinders, weight, displacement, horsepower)[train, ]
test.X = cbind(cylinders, weight, displacement, horsepower)[!train, ]
train.mpg01 = mpg01[train]
set.seed(1)
knn.pred = knn(train.X, test.X, train.mpg01, k = 1)
table(knn.pred, mpg01.test)## mpg01.test
## knn.pred 0 1
## 0 83 11
## 1 17 71mean(knn.pred != mpg01.test)## [1] 0.1538462\(K=10\)
knn.pred2 = knn(train.X, test.X, train.mpg01, k = 10)
table(knn.pred2, mpg01.test)## mpg01.test
## knn.pred2 0 1
## 0 77 7
## 1 23 75mean(knn.pred2 != mpg01.test)## [1] 0.1648352\(K=50\)
knn.pred3 = knn(train.X, test.X, train.mpg01, k = 50)
table(knn.pred3, mpg01.test)## mpg01.test
## knn.pred3 0 1
## 0 80 7
## 1 20 75mean(knn.pred3 != mpg01.test)## [1] 0.1483516After testing K values of 1, 10, and 50, I obtained error rates of about 15.38%, 16.48%, and 11.86%, respectively. This means that the K value of 50 performed the best on this data set.
Question 13
Using the Boston data set, fit classification models in order to predict whether a given suburb has a crime rate above or below the median.Explore logistic regression, LDA, and KNN models using various subsets of the predictors. Describe your findings.
library(MASS)
attach(Boston)
crim01 = rep(0, length(crim))
crim01[crim > median(crim)] = 1
Boston = data.frame(Boston, crim01)
train = 1:(length(crim) / 2)
test = (length(crim) / 2 + 1):length(crim)
Boston.train = Boston[train, ]
Boston.test = Boston[test, ]
crim01.test = crim01[test]Linear Regression: all variables
glm.crim01 = glm(crim01 ~ . - crim01 - crim, data = Boston, family = binomial, subset = train)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredprobs = predict(glm.crim01, Boston.test, type = "response")
glm.pred = rep(0, length(probs))
glm.pred[probs > 0.5] = 1
table(glm.pred, crim01.test)## crim01.test
## glm.pred 0 1
## 0 68 24
## 1 22 139mean(glm.pred != crim01.test)## [1] 0.1818182The first linear regression including all variables gives a test error rate of about 18.18%.
Linear Regression: selective variables
glm.crim01 = glm(crim01 ~ lstat+medv+tax+rad+age, data = Boston, family = binomial, subset = train)
probs = predict(glm.crim01, Boston.test, type = "response")
glm.pred = rep(0, length(probs))
glm.pred[probs > 0.5] = 1
table(glm.pred, crim01.test)## crim01.test
## glm.pred 0 1
## 0 83 21
## 1 7 142mean(glm.pred != crim01.test)## [1] 0.1106719The second linear regression, including the variables lstat, medv, tax, rad, and age, provides a test error rate of about 11.07%, making it more accurate than the first regression model.
LDA
lda.crim01 = lda(crim01 ~ dis+rad+nox, data = Boston, subset = train)
lda.pred = predict(lda.crim01, Boston.test)
table(lda.pred$class, crim01.test)## crim01.test
## 0 1
## 0 80 16
## 1 10 147mean(lda.pred$class != crim01.test)## [1] 0.1027668The LDA model, including the variables dis, nox, and rad, provides a test error rate of about 10.28%. This is the most accurate model so far, with almost 90% accuracy on estimating Boston crime rates.
KNN
\(K=1\)
train.X <- cbind(nox, rm, age, dis, tax, medv)[train, ]
test.X <- cbind(nox, rm, age, dis, tax, medv)[test, ]
train.crim01 <- crim01[train]
set.seed(1)
knn.pred <- knn(train.X, test.X, train.crim01, k = 1)
table(knn.pred, crim01.test)## crim01.test
## knn.pred 0 1
## 0 83 153
## 1 7 10mean(knn.pred != crim01.test)## [1] 0.6324111\(K=10\)
knn.pred2 <- knn(train.X, test.X, train.crim01, k = 10)
table(knn.pred2, crim01.test)## crim01.test
## knn.pred2 0 1
## 0 83 23
## 1 7 140mean(knn.pred2 != crim01.test)## [1] 0.1185771\(K=50\)
knn.pred3 <- knn(train.X, test.X, train.crim01, k = 50)
table(knn.pred3, crim01.test)## crim01.test
## knn.pred3 0 1
## 0 79 22
## 1 11 141mean(knn.pred3 != crim01.test)## [1] 0.1304348The KNN models, which contain the variables nox, rm, age, dis, tax, and medv, had K values of 1, 10, and 50, respectively. K=1 produced the worst model of all, with a test error rate of 63.24%, meaning the predictions are wrong more often than they are right. The models with a K value of 50 performed much better, with a test error rate of 13.04%, but the best fit for this data set was a K value of 10, giving us the smallest error rate of 11.86%.
Overall, the best model was the LDA model with selective variables, producing a test error rate of only 10.28%. This means that 89.72% of the time, this model will accurately predict if a Boston suburb falls above or below the median crime rate. This is much better than two of the KNN models, and a significant improvement from the original linear regression model containing all of the variables. Other notable, successful models are the second linear regression model with selective variables and the KNN model with a K value of 10. Both models produce predictions with 88.93% and 88.14% accuracy, respectively.