Assignment #3
Yuan Torres
2023-02-27
PROBLEM 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.
(a)Produce some numerical and graphical summaries of the weekly data, do there appear to be any patterns?
library(ISLR)
library(tidyverse)
attach(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
##
##
##
## 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.000000000plot(Volume)
There is a strong correlation between ‘year’ and ‘volume’ with a number
0.8419. Volume is increasing overtime, average number of shares that are
traded on a daily basis have been increasing 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.weekly <- glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume, data=Weekly,family=binomial )
summary(glm.weekly)##
## 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: 4all the predictors` p-value are bigger than 0.05,except for “Lag2” with a P-value 0.0296, that means only the “Lag2” is a statistical significance predictor.
(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.weekly,type="response")
contrasts(Direction)## Up
## Down 0
## Up 1glm.pred <- rep("Up" ,1089)
glm.pred[glm.probs <=0.5]="Down"
table(glm.pred, Direction)## Direction
## glm.pred Down Up
## Down 54 48
## Up 430 557mean(glm.pred==Direction )## [1] 0.5610652This illustrates that the model predicted the weekly market trend correctly 56.11% of the time.Observations From the confusion matrix shows that total of 484 Down and the model has predicted 54 correctly and 430 wrongly. Similarly, total of 605 Up, the model has predicted 557 correctly and 48 wrongly. Percentage of correct predicted for Up is (557,605) * 100 = 92%, which tells us the model is doing well at predicting for Up value. Calculating for Down is (54/484) * 100 = 11%, which tells us the model is not good at predicting for Down value.
(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.0910 <- Weekly[!train,]
dim(Weekly.0910)## [1] 104 9glm.fit=glm(Direction~Lag2, data = Weekly,family = binomial,subset = train)
glm.probs=predict(glm.fit, Weekly.0910, type="response")
glm.pred=rep("Down", 104)
glm.pred[glm.probs>0.5]="Up"
Direction.0910=Direction[!train]
table(glm.pred, Direction.0910)## Direction.0910
## glm.pred Down Up
## Down 9 5
## Up 34 56mean(glm.pred==Direction.0910)## [1] 0.625When spliting up the whole Weekly dataset into a training and test dataset, the model correctly predicted weekly trends at rate of 62.5%, which is a moderate improvement from the model that utilized the whole dataset. The confusion matrix shows that true positive rate for “up” is (56/61)100=91.8%. rate for “down” is (9/43)100=21%. This model was able to improve significantly on correctly predicting overall,upwards and downward trends thank the previous one.
(e) Repeat (d) using LDA.
library(MASS)
lda.fit=lda(Direction~Lag2,data = Weekly,subset = train)
lda.pred=predict(lda.fit,Weekly.0910)
names(lda.pred)## [1] "class" "posterior" "x"lda.class <- lda.pred$class
table(lda.class, Direction.0910)## Direction.0910
## lda.class Down Up
## Down 9 5
## Up 34 56mean(lda.class==Direction.0910)## [1] 0.625Using Linear Discriminant Analysis to develop a classifying model yielded similar results as the logistic regression model created in part D.
(f) Repeat (d) using QDA.
qda.fit=qda(Direction~Lag2 ,data=Weekly ,subset=train)
qda.class=predict(qda.fit , Weekly.0910)$class
table(qda.class ,Direction.0910)## Direction.0910
## qda.class Down Up
## Down 0 0
## Up 43 61mean(qda.class==Direction.0910)## [1] 0.5865385Quadratic Linear Analysis created a model with an accuracy of 58.65%, In other words, 41.35% is the test error rate. which is lower than the previous LDA methods. For weeks when the market goes down, the model is right about 0% of the time.
(g) Repeat (d) using KNN with K = 1.
library(class)
weekly.train=cbind(Lag2[train ])
weekly.test=cbind(Lag2[!train ])
train.Direction =Direction [train]
set.seed(1)
knn.pred=knn(weekly.train,weekly.test,train.Direction ,k=1)
table(knn.pred ,Direction.0910)## Direction.0910
## knn.pred Down Up
## Down 21 30
## Up 22 31mean(knn.pred == Direction.0910)## [1] 0.5The K-Nearest neighbors resulted in a classifying model with an accuracy rate of 50% which is equal to random chance.
(h) Repeat (d) using naive Bayes.
library (e1071)
nb.fit=naiveBayes(Direction~Lag2 ,data=Weekly ,subset=train)
nb.class=predict(nb.fit ,Weekly.0910)
table(nb.class ,Direction.0910)## Direction.0910
## nb.class Down Up
## Down 0 0
## Up 43 61mean (nb.class == Direction.0910)## [1] 0.5865385(i) Which of these methods appears to provide the best results on this data?
The methods that have the highest accuracy rates are the Logistic Regression and Linear Discriminant Analysis; both having rates of 62.5%.
(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.*
glm #2 full use Lag1-Lag3
glm.weekly <- glm(Direction~Lag1+Lag2+Lag3+Volume, data=Weekly,family=binomial )
glm.probs <- predict(glm.weekly,type="response")
glm.pred <- rep("Down" ,1089)
glm.pred[glm.probs>0.5]="Up"
table(glm.pred, Direction)## Direction
## glm.pred Down Up
## Down 46 44
## Up 438 561mean(glm.pred==Direction )## [1] 0.5573921compare to the first glm model with full dataset (0.5610652), this model (glm#2 with predictors fromLag1-Lag3) correct rate went down slightly.
full dataset , training data, testing data separated
train=(Year<2009)
Weekly.0910 <- Weekly[!train,]
dim(Weekly.0910)## [1] 104 9glm.fit=glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume, data = Weekly,family = binomial,subset = train)
glm.probs=predict(glm.fit, Weekly.0910, type="response")
glm.pred=rep("Down", 104)
glm.pred[glm.probs>0.5]="Up"
Direction.0910=Direction[!train]
table(glm.pred, Direction.0910)## Direction.0910
## glm.pred Down Up
## Down 31 44
## Up 12 17mean(glm.pred==Direction.0910)## [1] 0.4615385compare to (d)model, which used a separeted trainning data and test data with only one predictor, this model used all predictors but come out the correct rate is lower than (d)model (0.625)
lda with predictor lag2 squared
lda.fit=lda(Direction~Lag2^2,data = Weekly,subset = train)
lda.pred=predict(lda.fit,Weekly.0910)
names(lda.pred)## [1] "class" "posterior" "x"lda.class <- lda.pred$class
table(lda.class, Direction.0910)## Direction.0910
## lda.class Down Up
## Down 9 5
## Up 34 56mean(lda.class==Direction.0910)## [1] 0.625It comes out the same result as the LDA regular model(e).
**qda with predictor Lag2*Lag3**
qda.fit=qda(Direction~Lag2*Lag3 ,data=Weekly ,subset=train)
qda.class=predict(qda.fit , Weekly.0910)$class
table(qda.class ,Direction.0910)## Direction.0910
## qda.class Down Up
## Down 9 10
## Up 34 51mean(qda.class==Direction.0910)## [1] 0.5769231correct rate went down.
KNN model with K=15
weekly.train=cbind(Lag2[train ])
weekly.test=cbind(Lag2[!train ])
train.Direction =Direction [train]
set.seed(1)
knn.pred=knn(weekly.train,weekly.test,train.Direction ,k=15)
table(knn.pred ,Direction.0910)## Direction.0910
## knn.pred Down Up
## Down 20 20
## Up 23 41mean(knn.pred == Direction.0910)## [1] 0.5865385correct rate goes up by 8.6%, it is a big improvement. the number of K has a important influence in KNN model.
##PROBLEM 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(ISLR)
library(tidyverse)
attach(Auto)
dim(Auto)## [1] 392 9summary(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) :365head(Auto)## mpg cylinders displacement horsepower weight acceleration year origin
## 1 18 8 307 130 3504 12.0 70 1
## 2 15 8 350 165 3693 11.5 70 1
## 3 18 8 318 150 3436 11.0 70 1
## 4 16 8 304 150 3433 12.0 70 1
## 5 17 8 302 140 3449 10.5 70 1
## 6 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 500median(mpg)## [1] 22.75mpg01 <- ifelse( mpg > median(mpg), yes = 1, no = 0)
Auto01 <- data.frame(Auto, mpg01)
head(Auto01)## mpg cylinders displacement horsepower weight acceleration year origin
## 1 18 8 307 130 3504 12.0 70 1
## 2 15 8 350 165 3693 11.5 70 1
## 3 18 8 318 150 3436 11.0 70 1
## 4 16 8 304 150 3433 12.0 70 1
## 5 17 8 302 140 3449 10.5 70 1
## 6 15 8 429 198 4341 10.0 70 1
## name mpg01
## 1 chevrolet chevelle malibu 0
## 2 buick skylark 320 0
## 3 plymouth satellite 0
## 4 amc rebel sst 0
## 5 ford torino 0
## 6 ford galaxie 500 0dim(Auto01)## [1] 392 10(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(Auto01[,-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.0000000from this correlation matrix we can see, “cylinders, displancement, horsepower and weight” associates with “mpg01” more than others.
plot(mpg01, cylinders)
plot(mpg01,horsepower)
plot(mpg01, displacement)
plot(mpg01,weight)
plot(mpg01,origin )
Since mpg01 is a categorical field, it’s effects are not so easily
observed on a scatter plot.
ggplot() + geom_point(aes(horsepower, displacement, colour = as.factor(mpg01)))
ggplot() + geom_point(aes(horsepower,weight, colour = as.factor(mpg01)))
ggplot() + geom_point(aes(displacement,weight, colour = as.factor(mpg01)))
ggplot() + geom_point(aes(cylinders,displacement, colour = as.factor(mpg01)))
The plots show strong negative correlation between mpg01 and
cylinders, displacement, horsepower and weight.
boxplot(cylinders ~ mpg01, data = Auto01, main = "Cylinders vs mpg01")
boxplot(displacement ~ mpg01, data = Auto01, main = "Displacement vs mpg01")
boxplot(horsepower ~ mpg01, data = Auto01, main = "Horsepower vs mpg01")
boxplot(weight ~ mpg01, data = Auto01, main = "Weight vs mpg01")
boxplot(year ~ mpg01, data = Auto01, main = "Year vs mpg01")
(c) Split the data into a training set and a test set.
set.seed(1)
train <- sample(nrow(Auto01), size = nrow(Auto01) * 0.7)
auto.train<-Auto01[train, ]
auto.test<-Auto01[-train, ]
mpg01.test=mpg01[-train]
dim(Auto01)## [1] 392 10dim(auto.train)## [1] 274 10dim(auto.test)## [1] 118 10(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)
lda.new=lda(mpg01 ~ cylinders + weight + displacement + horsepower,
data = Auto01, subset = train)
lda.pred=predict(lda.new,auto.test)
names(lda.pred)## [1] "class" "posterior" "x"lda.class = lda.pred$class
table(lda.class, mpg01.test)## mpg01.test
## lda.class 0 1
## 0 50 3
## 1 11 54mean(lda.class!=mpg01.test)## [1] 0.1186441The LDA model does well. The prediction test error rate of 11.86%
(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.new=qda(mpg01 ~ cylinders + weight + displacement + horsepower,
data = Auto01, subset = train)
qda.class=predict(qda.new , auto.test)$class
table(qda.class ,mpg01.test)## mpg01.test
## qda.class 0 1
## 0 52 5
## 1 9 52mean(qda.class!=mpg01.test)## [1] 0.1186441The QDA model does the same as LDA, The prediction test error rate of 11.86%
(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.auto <- glm(mpg01 ~ cylinders + weight + displacement + horsepower,
data = Auto01,family=binomial,subset = train )
glm.probs <- predict(glm.auto,auto.test,type="response")
glm.pred <- rep("0" ,118)
glm.pred[glm.probs >0.5]="1"
table(glm.pred, mpg01.test)## mpg01.test
## glm.pred 0 1
## 0 53 3
## 1 8 54mean(glm.pred!=mpg01.test )## [1] 0.09322034The logistic regression model does better than the LDA and QDA model, The prediction test error rate is 9.32%
(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?
library (e1071)
nb.fit=naiveBayes(mpg01 ~ cylinders + weight + displacement + horsepower,
data = Auto01,subset = train,)
nb.class=predict(nb.fit ,auto.test)
table(nb.class ,mpg01.test)## mpg01.test
## nb.class 0 1
## 0 52 4
## 1 9 53mean (nb.class !=mpg01.test)## [1] 0.1101695The naive Bayes model did a little better than LDA and QDA, but not as good as logistic regression model. The prediction test error rate is 11.02%
(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?
library(class)
train <- sample(nrow(Auto01), size = nrow(Auto01) * 0.7)
train.k=cbind(cylinders,weight,displacement,horsepower)[train,]
test.k=cbind(cylinders,weight,displacement,horsepower)[-train,]
train.mpg01 =mpg01 [train]
mpg01.test=mpg01 [-train]
dim(train.k)## [1] 274 4dim(test.k)## [1] 118 4set.seed(1)
knn.pred=knn(train.k,test.k,train.mpg01 ,k=1)
table(knn.pred ,mpg01.test)## mpg01.test
## knn.pred 0 1
## 0 54 2
## 1 12 50mean(knn.pred != mpg01.test)## [1] 0.1186441When K=1, the predition test error rate is 13.56%, which is the highest among all of the previous models.
train.k=cbind(cylinders,weight,displacement,horsepower)[train,]
test.k=cbind(cylinders,weight,displacement,horsepower)[-train,]
train.mpg01 =mpg01 [train]
mpg01.test=mpg01 [-train]
dim(train.k)## [1] 274 4dim(test.k)## [1] 118 4set.seed(1)
knn.pred=knn(train.k,test.k,train.mpg01 ,k=2)
table(knn.pred ,mpg01.test)## mpg01.test
## knn.pred 0 1
## 0 53 5
## 1 13 47mean(knn.pred != mpg01.test)## [1] 0.1525424When K=2, the prediction test error rate is 11.86%, better than K=1
train.k=cbind(cylinders,weight,displacement,horsepower)[train,]
test.k=cbind(cylinders,weight,displacement,horsepower)[-train,]
train.mpg01 =mpg01 [train]
mpg01.test=mpg01 [-train]
dim(train.k)## [1] 274 4dim(test.k)## [1] 118 4set.seed(1)
knn.pred=knn(train.k,test.k,train.mpg01 ,k=3)
table(knn.pred ,mpg01.test)## mpg01.test
## knn.pred 0 1
## 0 58 4
## 1 8 48mean(knn.pred != mpg01.test)## [1] 0.1016949Error rate get better when K=3. with a # 11.01%
train.k=cbind(cylinders,weight,displacement,horsepower)[train,]
test.k=cbind(cylinders,weight,displacement,horsepower)[-train,]
train.mpg01 =mpg01 [train]
mpg01.test=mpg01 [-train]
dim(train.k)## [1] 274 4dim(test.k)## [1] 118 4set.seed(1)
knn.pred=knn(train.k,test.k,train.mpg01 ,k=5)
table(knn.pred ,mpg01.test)## mpg01.test
## knn.pred 0 1
## 0 56 4
## 1 10 48mean(knn.pred != mpg01.test)## [1] 0.1186441When K=5,the error rate get worse. 12.71%
train.k=cbind(cylinders,weight,displacement,horsepower)[train,]
test.k=cbind(cylinders,weight,displacement,horsepower)[-train,]
train.mpg01 =mpg01 [train]
mpg01.test=mpg01 [-train]
dim(train.k)## [1] 274 4dim(test.k)## [1] 118 4set.seed(1)
knn.pred=knn(train.k,test.k,train.mpg01 ,k=10)
table(knn.pred ,mpg01.test)## mpg01.test
## knn.pred 0 1
## 0 56 2
## 1 10 50mean(knn.pred != mpg01.test)## [1] 0.1016949With error rate 14.41% by K=10, we see that the best prediction test error rate is when K=3
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')
names(Boston)## [1] "crim" "zn" "indus" "chas" "nox" "rm" "age"
## [8] "dis" "rad" "tax" "ptratio" "black" "lstat" "medv"crime_m = median(Boston$crim)
crime_m## [1] 0.25651crime_rate=ifelse(Boston$crim > crime_m,1,0)
as.factor(crime_rate)## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0
## [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
## [112] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## [186] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0 0 0 1 1
## [223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 1
## [260] 1 1 1 1 1 1 1 1 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [297] 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 1 0 1 1 0 0 1 1 1 0 1 0 0 0 0 0 0
## [334] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [371] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [408] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [445] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [482] 1 1 1 1 1 1 1 0 0 0 0 0 0 1 0 1 1 0 0 0 0 0 0 0 0
## Levels: 0 1Boston.rate=data.frame(Boston,crime_rate)
head(Boston.rate)## crim zn indus chas nox rm age dis rad tax ptratio black lstat
## 1 0.00632 18 2.31 0 0.538 6.575 65.2 4.0900 1 296 15.3 396.90 4.98
## 2 0.02731 0 7.07 0 0.469 6.421 78.9 4.9671 2 242 17.8 396.90 9.14
## 3 0.02729 0 7.07 0 0.469 7.185 61.1 4.9671 2 242 17.8 392.83 4.03
## 4 0.03237 0 2.18 0 0.458 6.998 45.8 6.0622 3 222 18.7 394.63 2.94
## 5 0.06905 0 2.18 0 0.458 7.147 54.2 6.0622 3 222 18.7 396.90 5.33
## 6 0.02985 0 2.18 0 0.458 6.430 58.7 6.0622 3 222 18.7 394.12 5.21
## medv crime_rate
## 1 24.0 0
## 2 21.6 0
## 3 34.7 0
## 4 33.4 0
## 5 36.2 0
## 6 28.7 0set.seed(1)
train <- sample(nrow(Boston), size = nrow(Boston) * 0.7)
train.Boston<-Boston[train, ]
test.Boston<-Boston[-train, ]
crime_rate.test=crime_rate[-train]
dim(Boston)## [1] 506 14dim(train.Boston)## [1] 354 14dim(test.Boston)## [1] 152 14names(Boston)## [1] "crim" "zn" "indus" "chas" "nox" "rm" "age"
## [8] "dis" "rad" "tax" "ptratio" "black" "lstat" "medv"cor(Boston.rate)## crim zn indus chas nox
## crim 1.00000000 -0.20046922 0.40658341 -0.055891582 0.42097171
## zn -0.20046922 1.00000000 -0.53382819 -0.042696719 -0.51660371
## indus 0.40658341 -0.53382819 1.00000000 0.062938027 0.76365145
## chas -0.05589158 -0.04269672 0.06293803 1.000000000 0.09120281
## nox 0.42097171 -0.51660371 0.76365145 0.091202807 1.00000000
## rm -0.21924670 0.31199059 -0.39167585 0.091251225 -0.30218819
## age 0.35273425 -0.56953734 0.64477851 0.086517774 0.73147010
## dis -0.37967009 0.66440822 -0.70802699 -0.099175780 -0.76923011
## rad 0.62550515 -0.31194783 0.59512927 -0.007368241 0.61144056
## tax 0.58276431 -0.31456332 0.72076018 -0.035586518 0.66802320
## ptratio 0.28994558 -0.39167855 0.38324756 -0.121515174 0.18893268
## black -0.38506394 0.17552032 -0.35697654 0.048788485 -0.38005064
## lstat 0.45562148 -0.41299457 0.60379972 -0.053929298 0.59087892
## medv -0.38830461 0.36044534 -0.48372516 0.175260177 -0.42732077
## crime_rate 0.40939545 -0.43615103 0.60326017 0.070096774 0.72323480
## rm age dis rad tax
## crim -0.21924670 0.35273425 -0.37967009 0.625505145 0.58276431
## zn 0.31199059 -0.56953734 0.66440822 -0.311947826 -0.31456332
## indus -0.39167585 0.64477851 -0.70802699 0.595129275 0.72076018
## chas 0.09125123 0.08651777 -0.09917578 -0.007368241 -0.03558652
## nox -0.30218819 0.73147010 -0.76923011 0.611440563 0.66802320
## rm 1.00000000 -0.24026493 0.20524621 -0.209846668 -0.29204783
## age -0.24026493 1.00000000 -0.74788054 0.456022452 0.50645559
## dis 0.20524621 -0.74788054 1.00000000 -0.494587930 -0.53443158
## rad -0.20984667 0.45602245 -0.49458793 1.000000000 0.91022819
## tax -0.29204783 0.50645559 -0.53443158 0.910228189 1.00000000
## ptratio -0.35550149 0.26151501 -0.23247054 0.464741179 0.46085304
## black 0.12806864 -0.27353398 0.29151167 -0.444412816 -0.44180801
## lstat -0.61380827 0.60233853 -0.49699583 0.488676335 0.54399341
## medv 0.69535995 -0.37695457 0.24992873 -0.381626231 -0.46853593
## crime_rate -0.15637178 0.61393992 -0.61634164 0.619786249 0.60874128
## ptratio black lstat medv crime_rate
## crim 0.2899456 -0.38506394 0.4556215 -0.3883046 0.40939545
## zn -0.3916785 0.17552032 -0.4129946 0.3604453 -0.43615103
## indus 0.3832476 -0.35697654 0.6037997 -0.4837252 0.60326017
## chas -0.1215152 0.04878848 -0.0539293 0.1752602 0.07009677
## nox 0.1889327 -0.38005064 0.5908789 -0.4273208 0.72323480
## rm -0.3555015 0.12806864 -0.6138083 0.6953599 -0.15637178
## age 0.2615150 -0.27353398 0.6023385 -0.3769546 0.61393992
## dis -0.2324705 0.29151167 -0.4969958 0.2499287 -0.61634164
## rad 0.4647412 -0.44441282 0.4886763 -0.3816262 0.61978625
## tax 0.4608530 -0.44180801 0.5439934 -0.4685359 0.60874128
## ptratio 1.0000000 -0.17738330 0.3740443 -0.5077867 0.25356836
## black -0.1773833 1.00000000 -0.3660869 0.3334608 -0.35121093
## lstat 0.3740443 -0.36608690 1.0000000 -0.7376627 0.45326273
## medv -0.5077867 0.33346082 -0.7376627 1.0000000 -0.26301673
## crime_rate 0.2535684 -0.35121093 0.4532627 -0.2630167 1.00000000From the correlations matrix we can see none of the variables has a numeric significant relationship with crim rate. only the “chas” variable has the lowest P value of 0.07 among other variables but still about the 0.05 cut off point.
logistic regression(with one predictor)
glm.Boston<- glm(crime_rate ~chas,
data = Boston,family=binomial,subset = train )
glm.probs <- predict(glm.Boston,test.Boston,type="response")
glm.pred <- rep("0" ,152)
glm.pred[glm.probs >0.5]="1"
table(glm.pred, crime_rate.test)## crime_rate.test
## glm.pred 0 1
## 0 72 74
## 1 1 5mean(glm.pred!=crime_rate.test )## [1] 0.4934211mean(glm.pred==crime_rate.test )## [1] 0.5065789I only used the “chas” variable in this model since it has the lowest P value in the correlated matrix. it came out with a test error rate of 49.34% which means the model with one variable did not do a good job on average, but it accurately predicted 98.63% (which is 72/73) when the crime rate is low than median.
logistic regression(with multiple predictors)
glm.Boston<- glm(crime_rate ~chas+zn+rm+age+tax+ptratio+lstat+medv,
data = Boston,family=binomial,subset = train )
glm.probs <- predict(glm.Boston,test.Boston,type="response")
glm.pred <- rep("0" ,152)
glm.pred[glm.probs >0.5]="1"
table(glm.pred, crime_rate.test)## crime_rate.test
## glm.pred 0 1
## 0 65 17
## 1 8 62mean(glm.pred!=crime_rate.test )## [1] 0.1644737mean(glm.pred==crime_rate.test )## [1] 0.8355263logistic regression has test error rate with 16.48%. it did a better job with multiple predictors.model accuracy rate is 83.56% overall
LDA model (with one predictor)
lda.Boston=lda(crime_rate~medv, data = Boston,subset = train)
lda.pred=predict(lda.Boston,test.Boston)
names(lda.pred)## [1] "class" "posterior" "x"lda.class = lda.pred$class
table(lda.class, crime_rate.test)## crime_rate.test
## lda.class 0 1
## 0 38 20
## 1 35 59mean(lda.class!=crime_rate.test)## [1] 0.3618421LDA model with one predictor “medv” has a higher test error rate 36.18%. On the lower than median rate part it did a similar random guessing job(38/35)
LDA model (with more predictors)
lda.Boston=lda(crime_rate~zn+rm+age+tax+ptratio+lstat+medv,
data = Boston,subset = train)
lda.pred=predict(lda.Boston,test.Boston)
names(lda.pred)## [1] "class" "posterior" "x"lda.class = lda.pred$class
table(lda.class, crime_rate.test)## crime_rate.test
## lda.class 0 1
## 0 65 17
## 1 8 62mean(lda.class!=crime_rate.test)## [1] 0.1644737LDA model did better with multiple predictors.The test error rate 16.45% which is the same as losgitic regression model.
NaiveBayes model (with one predictor)
library (e1071)
nb.fit=naiveBayes(crime_rate~ptratio, data = Boston,subset = train)
nb.class=predict(nb.fit ,test.Boston)
table(nb.class ,crime_rate.test)## crime_rate.test
## nb.class 0 1
## 0 59 25
## 1 14 54mean (nb.class ==crime_rate.test)## [1] 0.743421159/(59+14)100=81% (lower than median rate), 54/(25+54)100=68% (higher than median rate). NB model with one predictor has a overall test accuate rate 74.34%. It predicted 81% of time right with “ptratio” variable on “lower than median rate” cases.
NaiveBayes model (with more predictors)
nb.fit=naiveBayes(crime_rate~zn+rm+age+tax+ptratio+lstat+medv,
data = Boston,subset = train)
nb.class=predict(nb.fit ,test.Boston)
table(nb.class ,crime_rate.test)## crime_rate.test
## nb.class 0 1
## 0 60 13
## 1 13 66mean (nb.class !=crime_rate.test)## [1] 0.1710526NB model with mutiple predictors has a slight higher test error rate than LDA model. it is 17.1%
KNN model, K=1
library(class)
attach(Boston)
train.Boston=cbind(zn,rm,age,tax,ptratio,lstat,medv)[train,]
test.Boston=cbind(zn,rm,age,tax,ptratio,lstat,medv)[-train,]
train.crime_rate =crime_rate[train]
crime_rate.test=crime_rate [-train]
dim(train.Boston)## [1] 354 7dim(test.Boston)## [1] 152 7set.seed(1)
knn.pred=knn(train.Boston,test.Boston,train.crime_rate ,k=1)
table(knn.pred ,crime_rate.test)## crime_rate.test
## knn.pred 0 1
## 0 65 5
## 1 8 74mean(knn.pred == crime_rate.test)## [1] 0.9144737When K=1, KNN model test accuracy rate is pretty good with a number 91.44%
KNN model, K=2
train.Boston=cbind(zn,rm,age,tax,ptratio,lstat,medv)[train,]
test.Boston=cbind(zn,rm,age,tax,ptratio,lstat,medv)[-train,]
train.crime_rate =crime_rate[train]
crime_rate.test=crime_rate [-train]
dim(train.Boston)## [1] 354 7dim(test.Boston)## [1] 152 7set.seed(1)
knn.pred=knn(train.Boston,test.Boston,train.crime_rate ,k=2)
table(knn.pred ,crime_rate.test)## crime_rate.test
## knn.pred 0 1
## 0 61 7
## 1 12 72mean(knn.pred == crime_rate.test)## [1] 0.875When K=2, the test accuracy rate went down to 87.5%
KNN model, K=5
train.Boston=cbind(zn,rm,age,tax,ptratio,lstat,medv)[train,]
test.Boston=cbind(zn,rm,age,tax,ptratio,lstat,medv)[-train,]
train.crime_rate =crime_rate[train]
crime_rate.test=crime_rate [-train]
dim(train.Boston)## [1] 354 7dim(test.Boston)## [1] 152 7set.seed(1)
knn.pred=knn(train.Boston,test.Boston,train.crime_rate ,k=5)
table(knn.pred ,crime_rate.test)## crime_rate.test
## knn.pred 0 1
## 0 63 6
## 1 10 73mean(knn.pred == crime_rate.test)## [1] 0.8947368When K=5, the test rate went up a little to 89.47%
KNN model, K=10,15,20
train.Boston=cbind(zn,rm,age,tax,ptratio,lstat,medv)[train,]
test.Boston=cbind(zn,rm,age,tax,ptratio,lstat,medv)[-train,]
train.crime_rate =crime_rate[train]
crime_rate.test=crime_rate [-train]
dim(train.Boston)## [1] 354 7dim(test.Boston)## [1] 152 7set.seed(1)
knn.pred=knn(train.Boston,test.Boston,train.crime_rate ,k=10)
table(knn.pred ,crime_rate.test)## crime_rate.test
## knn.pred 0 1
## 0 62 7
## 1 11 72mean(knn.pred == crime_rate.test)## [1] 0.8815789train.Boston=cbind(zn,rm,age,tax,ptratio,lstat,medv)[train,]
test.Boston=cbind(zn,rm,age,tax,ptratio,lstat,medv)[-train,]
train.crime_rate =crime_rate[train]
crime_rate.test=crime_rate [-train]
dim(train.Boston)## [1] 354 7dim(test.Boston)## [1] 152 7set.seed(1)
knn.pred=knn(train.Boston,test.Boston,train.crime_rate ,k=15)
table(knn.pred ,crime_rate.test)## crime_rate.test
## knn.pred 0 1
## 0 60 10
## 1 13 69mean(knn.pred == crime_rate.test)## [1] 0.8486842train.Boston=cbind(zn,rm,age,tax,ptratio,lstat,medv)[train,]
test.Boston=cbind(zn,rm,age,tax,ptratio,lstat,medv)[-train,]
train.crime_rate =crime_rate[train]
crime_rate.test=crime_rate [-train]
dim(train.Boston)## [1] 354 7dim(test.Boston)## [1] 152 7set.seed(1)
knn.pred=knn(train.Boston,test.Boston,train.crime_rate ,k=20)
table(knn.pred ,crime_rate.test)## crime_rate.test
## knn.pred 0 1
## 0 58 10
## 1 15 69mean(knn.pred == crime_rate.test)## [1] 0.8355263As K number goes up, test accuracy rate going down.from the result, KNN model has the best test accuracy rate when K = 1 and 5