Assignment #3
Alliah Kelley
2022-05-09
library(ISLR)
library(MASS)
library(class)
attach(Weekly)
attach(Boston)
attach(Auto)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?
pairs(Weekly[,1:8])
cor(Weekly[,1:8])## 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 looks to be a positive correlation between Year
and Volume with Volume increases over the
years.
(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?
Directionfit=glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume, data=Weekly, family="binomial")
summary(Directionfit)##
## 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: 4Lag2 is the only predictor that is 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.
Correctprob=predict(Directionfit,type="response")
Correctpred=rep("Down", length(Correctprob))
Correctpred[Correctprob>0.5]="Up"
table(Correctpred, Direction)## Direction
## Correctpred Down Up
## Down 54 48
## Up 430 557When the predictions are correct the percentage is 56.11% which is (54+557)/(54+557+48+430). Also when the market goes up its correct 92.07% (557/(557+48)) and then its incorrect 11.16% (54/(430+54)) most 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).
Training=(Year<2009)
Weekly0910=Weekly[!Training, ]
Trainingfit=glm(Direction~Lag2,data=Weekly,family=binomial,subset=Training)
Correctprob=predict(Trainingfit,Weekly0910,type="response")
Correctpred=rep("Down",length(Correctprob))
Correctpred[Correctprob>0.5]="Up"
Direction0910=Direction[!Training]
table(Correctpred,Direction0910)## Direction0910
## Correctpred Down Up
## Down 9 5
## Up 34 56mean(Correctpred==Direction0910)## [1] 0.625This shows that the correct predictions for the test data is 62.5%.
(e) Repeat (d) using LDA.
TrainingLDA=lda(Direction~Lag2,data=Weekly,subset=Training)
PredLDA=predict(TrainingLDA,Weekly0910)
table(PredLDA$class,Direction0910)## Direction0910
## Down Up
## Down 9 5
## Up 34 56mean(PredLDA$class==Direction0910)## [1] 0.625Using LDA shows that the correct predictions for the test data is also 62.5%.
(f) Repeat (d) using QDA.
TrainingQDA=qda(Direction~Lag2,data=Weekly,subset=Training)
PredQDA=predict(TrainingQDA,Weekly0910)$class
table(PredQDA,Direction0910)## Direction0910
## PredQDA Down Up
## Down 0 0
## Up 43 61mean(PredQDA==Direction0910)## [1] 0.5865385Using QDA shows that the correct predictions for the test data is 58.65%.
(g) Repeat (d) using KNN with K = 1.
Training1=as.matrix(Lag2[Training])
Test1=as.matrix(Lag2[!Training])
Trainingdirect=Direction[Training]
set.seed(1)
PredKNN=knn(Training1,Test1,Trainingdirect,k=1)
table(PredKNN,Direction0910)## Direction0910
## PredKNN Down Up
## Down 21 30
## Up 22 31mean(PredKNN==Direction0910)## [1] 0.5Using KNN shows that the correct preditions for the test data is 50%.
(h) Which of these methods appears to provide the best results on this data? Logistical regression and LDA have the best results with 62.5%.
(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:Lag3
Directionfit2=glm(Direction~Lag2:Lag3,data=Weekly,family=binomial,subset=Training)
Prob2=predict(Directionfit2,Weekly0910,type="response")
Pred2=rep("Down",length(Prob2))
Pred2[Prob2>0.5]="Up"
Direction0910=Direction[!Training]
table(Pred2, Direction0910)## Direction0910
## Pred2 Down Up
## Up 43 61mean(Pred2==Direction0910)## [1] 0.5865385# LDA with Lag2 and Lag3
TrainingLDA2=lda(Direction~Lag2:Lag3,data=Weekly,subset=Training)
PredLDA2=predict(TrainingLDA2,Weekly0910)
mean(PredLDA2$class==Direction0910)## [1] 0.5865385# QDA with sqrt of Lag3
TrainQDA2=qda(Direction~Lag3+sqrt(abs(Lag3)),data=Weekly,subset=Training)
PredQDA2=predict(TrainQDA2,Weekly0910)
table(PredQDA2$class,Direction0910)## Direction0910
## Down Up
## Down 9 8
## Up 34 53mean(PredQDA2$class==Direction0910)## [1] 0.5961538#KNN where k=5
PredKNN2=knn(Training1,Test1,Trainingdirect,k=5)
table(PredKNN2,Direction0910)## Direction0910
## PredKNN2 Down Up
## Down 15 20
## Up 28 41mean(PredKNN2==Direction0910)## [1] 0.5384615#KNN where k=50
PredKNN3=knn(Training1,Test1,Trainingdirect,k=50)
table(PredKNN3,Direction0910)## Direction0910
## PredKNN3 Down Up
## Down 21 20
## Up 22 41mean(PredKNN3==Direction0910)## [1] 0.5961538Using all these different combinations the Logistic Regression and LDA still have the best test error rates at 62.5%.
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.
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)
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(acceleration~mpg01,data=Auto,main = "Acceleration vs mpg01")
boxplot(year~mpg01,data=Auto,main="Year vs mpg01")
Looking at the boxplots show that there is some associations between
mpg01, cylinders, weight,
displacement and horsepower.
(c) Split the data into a training set and a test set.
training=(year %% 2 == 0)
Autotraining=Auto[training, ]
Autotesting=Auto[!training, ]
mpg01test=mpg01[!training](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?
trainingLDA=lda(mpg01~cylinders+weight+displacement+horsepower,data=Auto,subset=training)
trainingLDA## Call:
## lda(mpg01 ~ cylinders + weight + displacement + horsepower, data = Auto,
## subset = training)
##
## 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.0059035377predLDA=predict(trainingLDA, Autotesting)
table(predLDA$class, mpg01test)## mpg01test
## 0 1
## 0 86 9
## 1 14 73mean(predLDA$class != mpg01test)## [1] 0.1263736The test error rate of this model is 12.63%
(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?
trainingQDA=qda(mpg01~cylinders+weight+displacement+horsepower,data=Auto,subset=training)
trainingQDA## Call:
## qda(mpg01 ~ cylinders + weight + displacement + horsepower, data = Auto,
## subset = training)
##
## 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.92105predQDA=predict(trainingQDA,Autotesting)
table(predQDA$class,mpg01test)## mpg01test
## 0 1
## 0 89 13
## 1 11 69mean(predQDA$class != mpg01test)## [1] 0.1318681The test error rate of this model is 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?
trainingfit=glm(mpg01~cylinders+weight+displacement+horsepower,data=Auto,family=binomial,subset=training)
summary(trainingfit)##
## Call:
## glm(formula = mpg01 ~ cylinders + weight + displacement + horsepower,
## family = binomial, data = Auto, subset = training)
##
## 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: 7trainingprob=predict(trainingfit, Autotesting,type="response")
trainingpred=rep(0, length(trainingprob))
trainingpred[trainingprob>0.5]=1
table(trainingpred, mpg01test)## mpg01test
## trainingpred 0 1
## 0 89 11
## 1 11 71mean(trainingpred != mpg01test)## [1] 0.1208791The test error rate of this model is 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?
training1=cbind(cylinders,weight,displacement,horsepower)[training, ]
testing1=cbind(cylinders, weight, displacement, horsepower)[!training, ]
trainingmpg01=mpg01[training]
set.seed(1)
predknn=knn(training1,testing1,trainingmpg01,k=1)
table(predknn, mpg01test)## mpg01test
## predknn 0 1
## 0 83 11
## 1 17 71mean(predknn != mpg01test)## [1] 0.1538462The test error rate of this model is 15.38% with K = 1
predknn=knn(training1,testing1,trainingmpg01,k=10)
table(predknn, mpg01test)## mpg01test
## predknn 0 1
## 0 77 7
## 1 23 75mean(predknn != mpg01test)## [1] 0.1648352The test error rate of this model is 16.48% with K = 10
predknn=knn(training1,testing1,trainingmpg01,k=100)
table(predknn, mpg01test)## mpg01test
## predknn 0 1
## 0 81 7
## 1 19 75mean(predknn != mpg01test)## [1] 0.1428571The test error of this model is 14.29% with K = 100 and shows that a if K is 100 then it will perform the best.
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.
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
## Min. : 1.73 Min. : 5.00
## 1st Qu.: 6.95 1st Qu.:17.02
## Median :11.36 Median :21.20
## Mean :12.65 Mean :22.53
## 3rd Qu.:16.95 3rd Qu.:25.00
## Max. :37.97 Max. :50.00crime13=rep(0,length(crim))
crime13[crim>median(crim)]=1
Boston=data.frame(Boston, crime13)train13=1:(length(crim) /2)
test13=(length(crim) /2 + 1):length(crim)
Bostontraining=Boston[train13, ]
Bostontesting=Boston[test13, ]
crime13testing=crime13[test13]Bostonfit=glm(crime13 ~ . - crime13 - crim,data=Boston,family=binomial,subset=train13)#Logistic Regression
prob13=predict(Bostonfit,Bostontesting,type="response")
pred13=rep(0,length(prob13))
pred13[prob13>0.5]=1table(pred13,crime13testing)## crime13testing
## pred13 0 1
## 0 68 24
## 1 22 139mean(pred13 != crime13testing)## [1] 0.1818182This shows that we have a test error rate of 18.18% for logistic regression.
Bostonfit=glm(crime13 ~ . - crime13-crim-chas-nox,data=Boston,family=binomial,subset=train13)#Logistic Regression
prob13=predict(Bostonfit, Bostontesting,type="response")
pred13=rep(0, length(prob13))
pred13[prob13>0.5]=1mean(pred13 != crime13testing)## [1] 0.1581028This shows that we have a test error rate of 15.81% for logistic regression.
#LDA
BostonLDA=lda(crime13 ~ . - crime13-crim,data=Boston,subset=train13)
pred13LDA=predict(BostonLDA,Bostontesting)
table(pred13LDA$class, crime13testing)## crime13testing
## 0 1
## 0 80 24
## 1 10 139mean(pred13LDA$class != crime13testing)## [1] 0.1343874This shows that we have a test error rate of 13.44% for LDA.
BostonLDA=lda(crime13 ~ . - crime13-crim-chas-nox,data=Boston,subset=train13)
pred13LDA=predict(BostonLDA,Bostontesting)
table(pred13LDA$class, crime13testing)## crime13testing
## 0 1
## 0 82 30
## 1 8 133mean(pred13LDA$class != crime13testing)## [1] 0.1501976This shows that we have a test error rate of 15.02% for LDA.
#KNN = 1
KNNtrain=cbind(zn,indus,chas,nox,rm,age,dis,rad,tax,ptratio,black,lstat,medv)[train13, ]
KNNtest=cbind(zn,indus,chas,nox,rm,age,dis,rad,tax,ptratio,black,lstat,medv)[test13, ]
crime13train=crime13[train13]
set.seed(1)
KNNpred=knn(KNNtrain, KNNtest, crime13train, k=1)
table(KNNpred, crime13testing)## crime13testing
## KNNpred 0 1
## 0 85 111
## 1 5 52mean(KNNpred != crime13testing)## [1] 0.458498This shows that we have a test error rate of 45.85% if KNN = 1.
#KNN = 10
KNNpred=knn(KNNtrain, KNNtest, crime13train, k=10)
table(KNNpred, crime13testing)## crime13testing
## KNNpred 0 1
## 0 83 23
## 1 7 140mean(KNNpred != crime13testing)## [1] 0.1185771This shows that we have a test error rate of 11.86% if KNN = 10.
#KNN = 100
KNNpred=knn(KNNtrain, KNNtest, crime13train, k=100)
table(KNNpred, crime13testing)## crime13testing
## KNNpred 0 1
## 0 86 120
## 1 4 43mean(KNNpred != crime13testing)## [1] 0.4901186This shows that we have a test error rate of 11.86% if KNN = 100.
Looking at these different models we can see that using the KNN method actually performed the best and would be a good fit to predict whether a given suburb has a crime rate above or below the median.