Assignment 3
Alexis Chinnappan (UTSA ID: sez126)
6/27/2021
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?
names(Weekly)## [1] "Year" "Lag1" "Lag2" "Lag3" "Lag4" "Lag5"
## [7] "Volume" "Today" "Direction"dim(Weekly)## [1] 1089 9summary(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
##
##
##
## 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.000000000Since the Direction variable is qualitative, we remove for the cor() function, and it produces the pairwise correlations among the Weekly data set predictors. The correlations between the lag variables and Today are close to zero, this implies that there is very little correlation between today’s returns and previous day’s returns. Only Volume and Year has substantial correlation with each other.
attach(Weekly)
plot(Volume)
From the above plot, we are able to see that Volume of trades increase over time (that’s is, average number of shares traded weekly increased from 1990 to 2010).
(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.fits <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, data=Weekly, family = binomial )
summary(glm.fits)##
## 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: 4Since the p-value for the Lag2 (0.0296) is less than the significant value (0.05), the Lag2 predictor is statistically significant and all the other predictors’ p-values are greater than the significant value (0.05). So, those ARE NOT statistically significant.
(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.fits,type="response")
glm.pred <- rep("Down", length(glm.probs))
glm.pred[glm.probs > 0.5] <- "Up"
table(glm.pred, Direction)## Direction
## glm.pred Down Up
## Down 54 48
## Up 430 557The diagonal elements of the above confusion matrix indicate correct predictions, while the off-diagonals represent incorrect predictions. Hence the model correctly predicted that the market would go up on 557 weeks and that it would go down on 54 weeks, for a total of 557 + 54 = 611 correct predictions. The mean() function can be used to compute the fraction of days for which the prediction was correct. In this case, logistic regression correctly predicted the movement of the market 56.1% (611/1089) of the time. However, this result is misleading because we trained and tested the model on the same set of 1089 observations. In other words, 100− 56.1 = 43.9% is the training error rate, which is often overly optimistic.
Also, we can say that the model predicts correctly that for the weeks that the markets goes up is (557/(557+48)) 92.1% right. But, when the market goes down for the weeks, the prediction is only mere (54/(430+54)) 11.2% right. This is not correct 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 <- (Year < 2009)
Weekly.2009 <- Weekly[!train,]
Direction.2009 <- Direction[!train]
glm.train <- glm(Direction ~ Lag2, data = Weekly, family = binomial, subset = train)
summary(glm.train)##
## 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: 4glm.probs2 <- predict(glm.train, Weekly.2009, type = "response")
glm.pred2 <- rep("Down", length(glm.probs2))
glm.pred2[glm.probs2 > .5] <- "Up"
table(glm.pred2, Direction.2009)## Direction.2009
## glm.pred2 Down Up
## Down 9 5
## Up 34 56From the above confusion matrix, the correct predictions would be (56 + 9)=65, that is, (65/104) => 62.5% of the time the model predicts correctly and the test error rate would be 37.5%. For the weeks the market go up, the predictions is 91.8% (56/61) right and for the weeks the market goes down, the prediction is 20.9% (9/43) right. Still the down prediction is low in percent but this is better than the previous model with all the predictors.
(e) Repeat (d) using LDA.
lda.fit=lda(Direction ~ Lag2, data = Weekly, subset = train)
lda.fit## 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.4414162The LDA output indicates that ˆπ1 = 0.4477157 and ˆπ2 = 0.5522843; in other words, 44.8% of the training observations correspond to weeks during which the market went down and 55.2% of the training observations correspond to weeks during which market went up.
lda.pred <- predict(lda.fit, Weekly.2009)
table(lda.pred$class, Direction.2009)## Direction.2009
## Down Up
## Down 9 5
## Up 34 56From the above confusion matrix of the lda also reflects the similar percentages from the logistic regression model in (d).
(f) Repeat (d) using QDA.
qda.fit <- qda(Direction ~ Lag2, data = Weekly, subset = train)
qda.fit## 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.fit, Weekly.2009)
table(qda.pred$class, Direction.2009)## Direction.2009
## Down Up
## Down 0 0
## Up 43 61From the above confusion matrix, the correct predictions would be (61 + 0)=61, that is, (61/104) => 58.7% of the time the model predicts correctly and the test error rate would be 41.3%. For the weeks the market go up, the predictions is 100% (61/61) right and for the weeks the market goes down, the prediction is 0% (0/43) right. This qda model always opts for Up all the time even though it correctly predicts 58.7%.
(g) Repeat (d) using KNN with K = 1.
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.2009)## Direction.2009
## knn.pred Down Up
## Down 21 30
## Up 22 31From the above confusion matrix for knn model, the correct predictions would be (31 + 21)=52, that is, (52/104) => 50% of the time the model predicts correctly and the test error rate would be 50%. For the weeks the market go up, the predictions is 50.8% (31/61) right and for the weeks the market goes down, the prediction is 48.8% (21/43) right.
(h) Which of these methods appears to provide the best results on this data?
From the test error rates between the different methods, logistic regression and LDA have minimum test error rates compare to the others, that is logistic regression and LDA provides best results on this data. After this, QDA method provides decent error rate and followed by KNN.
(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.
Logistic regression with Lag2:Lag1
glm.fitlg12 <- glm(Direction ~ Lag2:Lag1, data = Weekly, family = binomial, subset = train)
glm.probslg12 <- predict(glm.fitlg12, Weekly.2009, type = "response")
glm.predlg12 <- rep("Down", length(glm.probslg12))
glm.predlg12[glm.probslg12 > 0.5] = "Up"
table(glm.predlg12, Direction.2009)## Direction.2009
## glm.predlg12 Down Up
## Down 1 1
## Up 42 60mean(glm.predlg12 == Direction.2009)## [1] 0.5865385LDA with Lag2:Lag5 interaction
lda.fitlg25 <- lda(Direction ~ Lag2:Lag5, data = Weekly, subset = train)
lda.predlg25 <- predict(lda.fitlg25, Weekly.2009)
table(lda.predlg25$class, Direction.2009)## Direction.2009
## Down Up
## Down 0 0
## Up 43 61mean(lda.predlg25$class == Direction.2009)## [1] 0.5865385QDA with Lag2:Lag3 interaction
qda.fitlg23 <- qda(Direction ~ Lag2:Lag3, data = Weekly, subset = train)
qda.predlg23 <- predict(qda.fitlg23, Weekly.2009)
table(qda.predlg23$class, Direction.2009)## Direction.2009
## Down Up
## Down 6 8
## Up 37 53mean(qda.predlg23$class == Direction.2009)## [1] 0.5673077QDA with sqrt(abs(Lag2))
qda.fitlg2 <- qda(Direction ~ Lag2 + sqrt(abs(Lag2)), data = Weekly, subset = train)
qda.predlg2 <- predict(qda.fitlg2, Weekly.2009)
table(qda.predlg2$class, Direction.2009)## Direction.2009
## Down Up
## Down 12 13
## Up 31 48mean(qda.predlg2$class == Direction.2009)## [1] 0.5769231KNN with k=10
knn.pred10 <- knn(train.X, test.X, train.Direction, k=10)
table(knn.pred10, Direction.2009)## Direction.2009
## knn.pred10 Down Up
## Down 17 18
## Up 26 43mean(knn.pred10 == Direction.2009)## [1] 0.5769231KNN with k=100
knn.pred100 <- knn(train.X, test.X, train.Direction, k=100)
table(knn.pred100, Direction.2009)## Direction.2009
## knn.pred100 Down Up
## Down 9 12
## Up 34 49mean(knn.pred100 == Direction.2009)## [1] 0.5576923When compare these model with the possibilities/combinations to the original Logistic regression model(d) and LDA (e), the logistic regression model and the LDA performs better than these with respect to the test error rates.
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.
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) :365attach(Auto)
mpg01 <- rep(0, length(mpg))
mpg01[mpg > median(mpg)] <- 1
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.
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.0000000pairs(Auto)
From the above correlation and pair-wise scatterplots, we may conclude that mpg01 might have association with cylinders, displacement, horsepower and weight. Since mpg01 is derived from mpg, mpg01 will have association with mpg.
boxplot(cylinders ~ mpg01, data = Auto, main = "cylinders vs mpg01")
boxplot(displacement ~ mpg01, data = Auto, main = "displacement vs mpg01")
boxplot(horsepower ~ mpg01, data = Auto, main = "horsepower vs mpg01")
boxplot(weight ~ mpg01, data = Auto, main = "weight vs mpg01")
boxplot(mpg ~ mpg01, data = Auto, main = "mpg vs mpg01")
(c) Split the data into a training set and a test set.
Splitting the Auto data into training set and a test set using the year, that is, even year as training set and odd year as test set.
year.train <- (year %% 2 == 0)
auto.train <- Auto[year.train, ]
auto.test <- Auto[!year.train, ]
mpg01.test <- mpg01[!year.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?
From (b), we identified cylinders, displacement, horsepower and weight has some association with mpg01.
fit.auto.train.lda <- lda(mpg01 ~ cylinders + displacement + horsepower + weight, data = Auto, subset = year.train)
fit.auto.train.lda## Call:
## lda(mpg01 ~ cylinders + displacement + horsepower + weight, data = Auto,
## subset = year.train)
##
## Prior probabilities of groups:
## 0 1
## 0.4571429 0.5428571
##
## Group means:
## cylinders displacement horsepower weight
## 0 6.812500 271.7396 133.14583 3604.823
## 1 4.070175 111.6623 77.92105 2314.763
##
## Coefficients of linear discriminants:
## LD1
## cylinders -0.6741402638
## displacement 0.0004481325
## horsepower 0.0059035377
## weight -0.0011465750lda.auto.test.pred <- predict(fit.auto.train.lda, auto.test)
table(lda.auto.test.pred$class, mpg01.test)## mpg01.test
## 0 1
## 0 86 9
## 1 14 73From the above confusion matrix, the correct predictions would be (86 + 73)=159, that is, (159/182) => 87.4% of the time the model predicts correctly and the test error rate would be 12.6%.
mean(lda.auto.test.pred$class != mpg01.test)## [1] 0.1263736From the above mean also we can find out the test error rate, which is approximately 12.6%.
(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?
fit.auto.train.qda <- qda(mpg01 ~ cylinders + displacement + horsepower + weight, data = Auto, subset = year.train)
fit.auto.train.qda## Call:
## qda(mpg01 ~ cylinders + displacement + horsepower + weight, data = Auto,
## subset = year.train)
##
## Prior probabilities of groups:
## 0 1
## 0.4571429 0.5428571
##
## Group means:
## cylinders displacement horsepower weight
## 0 6.812500 271.7396 133.14583 3604.823
## 1 4.070175 111.6623 77.92105 2314.763qda.auto.test.pred <- predict(fit.auto.train.qda, auto.test)
table(qda.auto.test.pred$class, mpg01.test)## mpg01.test
## 0 1
## 0 89 13
## 1 11 69From the above confusion matrix, the correct predictions would be (89 + 69)=158, that is, (158/182) => 86.8% of the time the model predicts correctly and the test error rate would be 13.2%.
mean(qda.auto.test.pred$class != mpg01.test)## [1] 0.1318681From the above mean also we can find out the test error rate, which is approximately 13.2%.
(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?
fit.auto.train.glm <- glm(mpg01 ~ cylinders + displacement + horsepower + weight, data = Auto, family = binomial, subset = year.train)
summary(fit.auto.train.glm)##
## Call:
## glm(formula = mpg01 ~ cylinders + displacement + horsepower +
## weight, family = binomial, data = Auto, subset = year.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
## displacement 0.002462 0.015030 0.164 0.8699
## horsepower -0.050611 0.025209 -2.008 0.0447 *
## weight -0.002922 0.001137 -2.569 0.0102 *
## ---
## 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: 7glm.auto.test.probs <- predict(fit.auto.train.glm, auto.test, type = "response")
glm.auto.test.pred <- rep(0, length(glm.auto.test.probs))
glm.auto.test.pred[glm.auto.test.probs > 0.5] <- 1
table(glm.auto.test.pred, mpg01.test)## mpg01.test
## glm.auto.test.pred 0 1
## 0 89 11
## 1 11 71From the above confusion matrix, the correct predictions would be (89 + 71)=160, that is, (160/182) => 87.9% of the time the model predicts correctly and the test error rate would be 12.1%.
mean(glm.auto.test.pred != mpg01.test)## [1] 0.1208791From the above mean also we can find out the test error rate, which is approximately 12.1%.
(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?
auto.train.X <- cbind(cylinders, displacement, horsepower, weight)[year.train, ]
auto.test.X <- cbind(cylinders, displacement, horsepower, weight)[!year.train, ]
auto.train.mpg01 <- mpg01[year.train]
set.seed(1)
knn.pred1 <- knn(auto.train.X, auto.test.X, auto.train.mpg01, k = 1)
table(knn.pred1, mpg01.test)## mpg01.test
## knn.pred1 0 1
## 0 83 11
## 1 17 71mean(knn.pred1 != mpg01.test)## [1] 0.1538462From the above, we can find out the test error rate, which is approximately 15.4% for K=1.
knn.pred10 <- knn(auto.train.X, auto.test.X, auto.train.mpg01, k = 10)
table(knn.pred10, mpg01.test)## mpg01.test
## knn.pred10 0 1
## 0 77 7
## 1 23 75mean(knn.pred10 != mpg01.test)## [1] 0.1648352From the above, we can find out the test error rate, which is approximately 16.5% for K=10.
knn.pred20 <- knn(auto.train.X, auto.test.X, auto.train.mpg01, k = 20)
table(knn.pred20, mpg01.test)## mpg01.test
## knn.pred20 0 1
## 0 81 7
## 1 19 75mean(knn.pred20 != mpg01.test)## [1] 0.1428571From the above, we can find out the test error rate, which is approximately 13.2% for K=20.
knn.pred30 <- knn(auto.train.X, auto.test.X, auto.train.mpg01, k = 30)
table(knn.pred30, mpg01.test)## mpg01.test
## knn.pred30 0 1
## 0 83 8
## 1 17 74mean(knn.pred30 != mpg01.test)## [1] 0.1373626From the above, we can find out the test error rate, which is approximately 13.7% for K=30.
knn.pred100 <- knn(auto.train.X, auto.test.X, auto.train.mpg01, k = 100)
table(knn.pred100, mpg01.test)## mpg01.test
## knn.pred100 0 1
## 0 81 7
## 1 19 75mean(knn.pred100 != mpg01.test)## [1] 0.1428571From the above, we can find out the test error rate, which is approximately 14.3% for K=100.
From above all K values (that’s 1, 10, 20, 30, 100), K = 20 with 13.2% error rate, performs the best among all the other K values.
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.
Create a new predictor for crime rate crim01 with 0 if the crime rate is below the median and 1 if the crime rate is above the median.
attach(Boston)
crim01 <- rep(0, length(crim))
crim01[crim > median(crim)] <- 1
Boston <- data.frame(Boston, crim01)
summary(Boston)## crim zn indus chas
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 Min. :0.00000
## 1st Qu.: 0.08205 1st Qu.: 0.00 1st Qu.: 5.19 1st Qu.:0.00000
## Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.00000
## Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.06917
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.00000
## Max. :88.97620 Max. :100.00 Max. :27.74 Max. :1.00000
## nox rm age dis
## Min. :0.3850 Min. :3.561 Min. : 2.90 Min. : 1.130
## 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02 1st Qu.: 2.100
## Median :0.5380 Median :6.208 Median : 77.50 Median : 3.207
## Mean :0.5547 Mean :6.285 Mean : 68.57 Mean : 3.795
## 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08 3rd Qu.: 5.188
## Max. :0.8710 Max. :8.780 Max. :100.00 Max. :12.127
## rad tax ptratio black
## Min. : 1.000 Min. :187.0 Min. :12.60 Min. : 0.32
## 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 1st Qu.:375.38
## Median : 5.000 Median :330.0 Median :19.05 Median :391.44
## Mean : 9.549 Mean :408.2 Mean :18.46 Mean :356.67
## 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.23
## Max. :24.000 Max. :711.0 Max. :22.00 Max. :396.90
## lstat medv crim01
## Min. : 1.73 Min. : 5.00 Min. :0.0
## 1st Qu.: 6.95 1st Qu.:17.02 1st Qu.:0.0
## Median :11.36 Median :21.20 Median :0.5
## Mean :12.65 Mean :22.53 Mean :0.5
## 3rd Qu.:16.95 3rd Qu.:25.00 3rd Qu.:1.0
## Max. :37.97 Max. :50.00 Max. :1.0Split the Boston data set into 2 (70% train and 30% test).
train <- sort(sample(nrow(Boston), nrow(Boston)*.7))
boston.train <- Boston[train, ]
boston.test <- Boston[-train, ]
crim01.test <- crim01[-train]Run the logistic regression on the predictors except crim and predict the model.
glm.boston.fit <- glm (crim01 ~ . - crim, data = Boston, family = binomial, subset = train)
glm.boston.probs <- predict(glm.boston.fit, boston.test, type = "response")
glm.boston.pred <- rep(0, length(glm.boston.probs))
glm.boston.pred[glm.boston.probs > 0.5] <- 1
table(glm.boston.pred, crim01.test)## crim01.test
## glm.boston.pred 0 1
## 0 65 8
## 1 6 73mean(glm.boston.pred != crim01.test)## [1] 0.09210526For the above logistic regression, the test error rate is approximately 14.5%.
Let’s run the LDA with similar parameters of logistic regression.
lda.boston.fit <- lda(crim01 ~ . - crim, data = Boston, subset = train)
lda.boston.pred <- predict(lda.boston.fit, boston.test)
table(lda.boston.pred$class, crim01.test)## crim01.test
## 0 1
## 0 68 19
## 1 3 62mean(lda.boston.pred$class != crim01.test)## [1] 0.1447368For the above lda model, the test error rate is approximately 16.4%.
Let’s remove some of the predictors from the above lda model (that’s zn, rm and dis) along with crim.
lda.boston.fit2 <- lda(crim01 ~ . - crim - zn - rm - dis, data = Boston, subset = train)
lda.boston.pred2 <- predict(lda.boston.fit2, boston.test)
table(lda.boston.pred2$class, crim01.test)## crim01.test
## 0 1
## 0 70 18
## 1 1 63mean(lda.boston.pred2$class != crim01.test)## [1] 0.125For the above lda model, the test error rate is approximately 15.9%.
Now, we include all the parameters except crim in the train and test for KNN and apply KNN classifications for various K values starting with K = 1.
train.X <- cbind(zn, indus, chas, nox, rm, age, dis, rad, tax, ptratio, black, lstat, medv)[train, ]
test.X <- cbind(zn, indus, chas, nox, rm, age, dis, rad, tax, ptratio, black, lstat, medv)[-train, ]
train.crim01 <- crim01[train]
set.seed(1)
knn.pred.k1 <- knn(train.X, test.X, train.crim01, k = 1)
table(knn.pred.k1, crim01.test)## crim01.test
## knn.pred.k1 0 1
## 0 63 8
## 1 8 73mean(knn.pred.k1 != crim01.test)## [1] 0.1052632For the above KNN with K = 1, the test error rate is approximately 9.9%.
knn.pred.k5 <- knn(train.X, test.X, train.crim01, k = 5)
table(knn.pred.k5, crim01.test)## crim01.test
## knn.pred.k5 0 1
## 0 64 10
## 1 7 71mean(knn.pred.k5 != crim01.test)## [1] 0.1118421For the above KNN with K = 5, the test error rate is approximately 8.6%.
knn.pred.k10 <- knn(train.X, test.X, train.crim01, k = 10)
table(knn.pred.k10, crim01.test)## crim01.test
## knn.pred.k10 0 1
## 0 63 10
## 1 8 71mean(knn.pred.k10 != crim01.test)## [1] 0.1184211For the above KNN with K = 10, the test error rate is approximately 9.9%.
knn.pred.k20 <- knn(train.X, test.X, train.crim01, k = 20)
table(knn.pred.k20, crim01.test)## crim01.test
## knn.pred.k20 0 1
## 0 61 13
## 1 10 68mean(knn.pred.k20 != crim01.test)## [1] 0.1513158For the above KNN with K = 20, the test error rate is approximately 13.2%.
knn.pred.k100 <- knn(train.X, test.X, train.crim01, k = 100)
table(knn.pred.k100, crim01.test)## crim01.test
## knn.pred.k100 0 1
## 0 65 22
## 1 6 59mean(knn.pred.k100 != crim01.test)## [1] 0.1842105For the above KNN with K = 100, the test error rate is approximately 18.4%.
From above all K values (that’s 1, 5, 10, 20, 100), K = 5 with 8.6% error rate, performs the best among all the other K values.
From all the above classification models, the KNN with K = 5 model performs the best among all the others with lowest test error rate of 8.6%.