Assignment3
Matthew Westley
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
library(ISLR)
library(MASS)
library(class)Question 10
(10a) Produce some numerical and graphical summaries of the Weekly data. Do there appear to be any patterns?
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
##
##
##
## pairs(Weekly)
The only trend noticeable is the Volume of shares traded over time. The trading Volume grows at an increasing pace as the years pass.
plot(Weekly$Volume)
(10b) 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?
Lag2 appears to be statistically significant based off of the p-value of 0.0296, which indicates that Lag2 has a statistically significant relationship with Direction.
fit.direction <- glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume, data = Weekly, family = binomial)
summary(fit.direction)##
## 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: 4(10c) 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.
The model correctly predicted the Up/Down trends 56.10652% of the time, which is only slightly better than if the model had predicted Up every time (55.56%), so the model is not great.
ConfProb <- predict(fit.direction, type = 'response')
ConfPred <- rep('Down', length(ConfProb))
ConfPred[ConfProb > 0.5] = 'Up'
table(ConfPred, Weekly$Direction)##
## ConfPred Down Up
## Down 54 48
## Up 430 557(54+557)/(54+557+430+48)## [1] 0.5610652(10d) 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). This model correctly predicted the outcome 62.5% of the time.
training <- (Weekly$Year < 2009)fit.lag2 <- glm(Direction~Lag2, data = Weekly, subset = training, family = 'binomial')
summary(fit.lag2)##
## Call:
## glm(formula = Direction ~ Lag2, family = "binomial", data = Weekly,
## subset = training)
##
## 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: 4testing <- Weekly[!training,]
Lag2Prob <- predict(fit.lag2, testing, type = 'response')
Lag2Pred <- rep('Down', length(Lag2Prob))
Lag2Pred[Lag2Prob > 0.5] = 'Up'
Dirtest <- Weekly$Direction[!training]
table(Lag2Pred, Dirtest)## Dirtest
## Lag2Pred Down Up
## Down 9 5
## Up 34 56(9+56)/(9+5+34+56)## [1] 0.625(10e) Repeat (d) using LDA.
The LDA model returned the same results as the Logistic Regression model, and therefore had the same 62.5% correct prediction.
ldaFit.lag2 <- lda(Direction~Lag2, data = Weekly, subset = training, family = 'binomial')
ldaPred.lag2 <- predict(ldaFit.lag2, testing, type = 'response')
table(ldaPred.lag2$class, Dirtest)## Dirtest
## Down Up
## Down 9 5
## Up 34 56(10f) Repeat (d) using QDA.
The QDA model perfomed worse than the LDA and Logistic Regression models, scoring a lower prediction percentage (58.65385%)
qdaFit.lag2 <- qda(Direction~Lag2, data = Weekly, subset = training, family = 'binomial')
qdaPred.lag2 <- predict(qdaFit.lag2, testing, type = 'response')
table(qdaPred.lag2$class, Dirtest)## Dirtest
## Down Up
## Down 0 0
## Up 43 6161 / 104## [1] 0.5865385(10g) Repeat (d) using KNN with K = 1.
The KNN model with K = 1 is correct 50% of the time.
trainingKNN.data <- Weekly[Weekly$Year<2009,]
testKNN.data <- Weekly[Weekly$Year>2008,]
train.X <- cbind(trainingKNN.data$Lag2)
test.X <- cbind(testKNN.data$Lag2)
train.Y <- cbind(trainingKNN.data$Direction)
set.seed(1)
knn.pred <- knn(train.X, test.X, train.Y, k=1)
table(knn.pred, testKNN.data$Direction)##
## knn.pred Down Up
## 1 21 30
## 2 22 31(21 + 31) / (21+30+22+31)## [1] 0.5(10h) Which of these methods appears to provide the best results on this data?
The Logistic Regression and LDA models appeared to provide the best results (62.5%) on this data.
(10i) 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. The best of the various models I ran with different variables and K values was KNN Lag 4 with K = 4 and LDA with Lag1
ldaFit.lag1 <- lda(Direction~Lag1, data = Weekly, subset = training, family = 'binomial')
ldaPred.lag1 <- predict(ldaFit.lag1, testing, type = 'response')
table(ldaPred.lag1$class, Dirtest)## Dirtest
## Down Up
## Down 4 6
## Up 39 55(4 + 55) / (4 + 6 + 39 + 55)## [1] 0.5673077GLM with Lag4, Lag3, Lag2, and Volume
glm.lag3Vol <- lda(Direction~Lag4 + Lag3 + Lag2 + Volume, data = Weekly, subset = training, family = 'binomial')
glmPred.lag3Vol <- predict(glm.lag3Vol, testing, type = 'response')
table(glmPred.lag3Vol$class, Dirtest)## Dirtest
## Down Up
## Down 26 35
## Up 17 26(26 + 26) / (26 + 35 + 17 + 26)## [1] 0.5KNN K = 4
train.X = cbind(trainingKNN.data$Lag4)
test.X = cbind(testKNN.data$Lag4)
train.Y = cbind(trainingKNN.data$Direction)
set.seed(1)
knn.pred = knn(train.X, test.X, train.Y, k=4)
table(knn.pred, testKNN.data$Direction)##
## knn.pred Down Up
## 1 20 22
## 2 23 39(20 + 39) / (20+22+23+39)## [1] 0.567307711. 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.
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) :365(11a) 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.
auto <- Auto
mpg01 <- rep(0, dim(auto)[1])
mpg01[auto$mpg > median(auto$mpg)] = 1
auto = data.frame(Auto, mpg01)head(auto, n = 20)## 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
## 7 14 8 454 220 4354 9.0 70 1
## 8 14 8 440 215 4312 8.5 70 1
## 9 14 8 455 225 4425 10.0 70 1
## 10 15 8 390 190 3850 8.5 70 1
## 11 15 8 383 170 3563 10.0 70 1
## 12 14 8 340 160 3609 8.0 70 1
## 13 15 8 400 150 3761 9.5 70 1
## 14 14 8 455 225 3086 10.0 70 1
## 15 24 4 113 95 2372 15.0 70 3
## 16 22 6 198 95 2833 15.5 70 1
## 17 18 6 199 97 2774 15.5 70 1
## 18 21 6 200 85 2587 16.0 70 1
## 19 27 4 97 88 2130 14.5 70 3
## 20 26 4 97 46 1835 20.5 70 2
## 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 0
## 7 chevrolet impala 0
## 8 plymouth fury iii 0
## 9 pontiac catalina 0
## 10 amc ambassador dpl 0
## 11 dodge challenger se 0
## 12 plymouth 'cuda 340 0
## 13 chevrolet monte carlo 0
## 14 buick estate wagon (sw) 0
## 15 toyota corona mark ii 1
## 16 plymouth duster 0
## 17 amc hornet 0
## 18 ford maverick 0
## 19 datsun pl510 1
## 20 volkswagen 1131 deluxe sedan 1(11b) 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.
Cylinders, displacement, weight, acceleration, and year seem likely to be useful in predicting mpg01.
pairs(auto)
par(mfrow=c(2,3))
boxplot(auto$mpg01, auto$horsepower, main = "HP and mpg01")
boxplot(auto$mpg01, auto$cylinders, main = "Cylinders and mpg01")
boxplot(auto$mpg01, auto$weight, main = "Weight and mpg01")
boxplot(auto$mpg01, auto$displacement, main = "Displacement and mpg01")
boxplot(auto$mpg01, auto$acceleration, main = "Acceleration and mpg01")
(11c) Split the data into a training set and a test set.
attach(Auto)Auto <- data.frame(mpg01, apply(cbind(cylinders, weight, displacement, horsepower), 2, scale), year)train <- (year %% 2 == 0)
test <- !train
Auto.train <- Auto[train,]
Auto.test <- Auto[test,]
mpg01.test <- mpg01[test](11d) 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? The test error for this LDA model is 12.63736%
Autolda.fit <- lda(mpg01 ~ cylinders + weight + displacement + horsepower, data = Auto, subset = train)
Autolda.pred <- predict(Autolda.fit, Auto.test)
table(Autolda.pred$class, mpg01.test)## mpg01.test
## 0 1
## 0 86 9
## 1 14 73mean(Autolda.pred$class != mpg01.test)## [1] 0.1263736(11e) 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?
The test error of the QDA model is 13.18681%
Autoqda.fit <- qda(mpg01 ~ cylinders + weight + displacement + horsepower, data = Auto, subset = train)
Autoqda.pred <- predict(Autoqda.fit, Auto.test)
table(Autoqda.pred$class, mpg01.test)## mpg01.test
## 0 1
## 0 89 13
## 1 11 69mean(Autoqda.pred$class != mpg01.test)## [1] 0.1318681(11f) Perform logistic regression on the training data in order to pre- dict mpg01 using the variables that seemed most associated with mpg01 in (b). What is the test error of the model obtained?
The test error of the logistic regression model is 12.63736%
Autoglm.fit <- glm(mpg01 ~ cylinders + weight + displacement + horsepower, data = Auto, subset = train, family = binomial, )
Autoglm.prob <- predict(Autoglm.fit, Auto.test)
Autoglm.pred <- rep(0, length(Autoglm.prob))
Autoglm.pred[Autoglm.prob > 0.5] = 1
table(Autoglm.pred, mpg01.test)## mpg01.test
## Autoglm.pred 0 1
## 0 90 13
## 1 10 69mean(Autoglm.pred != mpg01.test)## [1] 0.1263736Autoglm.prob <- predict(Autoglm.fit, Auto[-train, ], type = "response")
Autoglm.pred <- rep(0, dim(Auto[-train, ])[1])
Autoglm.pred[Autoglm.prob > 0.5] = 1(11g) 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?
11.54% for K=1, 12.64% for K=10, 11.54% for K=5, 9.89% for K=4
KNN K=4 seems to have the lowest percentage of test errors.
KNNtrain.X <- cbind(cylinders, weight, displacement, horsepower)[train,]
KNNtest.X <- cbind(cylinders, weight, displacement, horsepower)[test,]
train.mpg01 <- mpg01[train]
set.seed(1)knn.pred <- knn(KNNtrain.X, KNNtest.X, train.mpg01, k = 1)
mean(knn.pred != mpg01.test)## [1] 0.1538462knn.pred <- knn(KNNtrain.X, KNNtest.X, train.mpg01, k = 10)
mean(knn.pred != mpg01.test)## [1] 0.1648352knn.pred <- knn(KNNtrain.X, KNNtest.X, train.mpg01, k = 5)
mean(knn.pred != mpg01.test)## [1] 0.1483516knn.pred <- knn(KNNtrain.X, KNNtest.X, train.mpg01, k = 4)
mean(knn.pred != mpg01.test)## [1] 0.1483516Question 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 sub- sets of the predictors. Describe your findings. The Logistic Regression with variables (dis + nox + rad) had a fairly low error rate. KNN with K=1 and the variables (dis, nox, age, medv) had a pretty high error rate, but the error rate improved greatly when K=5. Of the models I ran, LDA with the variables (dis + nox + age + tax + rad + medv) had the lowest error rate.
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.00attach(Boston)crimes <- rep(0, length(crim))
crimes[crim > median(crim)] <- 1
Boston <- data.frame(Boston, crimes)Ctrain <- 1:(dim(Boston)[1]/2)
Ctest <- (dim(Boston)[1]/2 + 1):dim(Boston)[1]
Bost.train <- Boston[Ctrain, ]
Bost.test <- Boston[Ctest, ]
crimes.test <- crimes[Ctest]par(mfrow=c(2,3))
boxplot(Boston$crimes, Boston$dis, main = "dis and crimes")
boxplot(Boston$crimes, Boston$nox, main = "nox and crimes")
boxplot(Boston$crimes, Boston$age, main = "age and crimes")
boxplot(Boston$crimes, Boston$tax, main = "tax and crimes")
boxplot(Boston$crimes, Boston$rad, main = "rad and crimes")
boxplot(Boston$crimes, Boston$medv, main = "medv and crimes")
Logistic Regression
set.seed(1)
BostonGLM.fit <- glm(crimes~dis + nox + rad, data = Bost.train, family = binomial)
Bost.prob <- predict(BostonGLM.fit, Bost.test, type = "response")
Bost.pred <- rep(0, length(Bost.prob))
Bost.pred[Bost.prob > 0.5] = 1
table(Bost.pred, crimes.test)## crimes.test
## Bost.pred 0 1
## 0 80 16
## 1 10 147mean(Bost.pred != crimes.test)## [1] 0.1027668LDA
Bostlda.fit <- lda(crimes ~ dis + nox + age + tax + rad + medv, data = Bost.train, family=binomial)
Bostlda.pred <- predict(Bostlda.fit, Bost.test)
table(Bostlda.pred$class, crimes.test)## crimes.test
## 0 1
## 0 81 16
## 1 9 147mean(Bostlda.pred$class != crimes.test)## [1] 0.09881423KNN K=1
KNtrain <- cbind(dis, nox, age, medv)[Ctrain, ]
KNtest <- cbind(dis, nox, age, medv)[Ctest, ]
BostKNN.pred <- knn(KNtrain, KNtest, crimes.test, k=1)
table(BostKNN.pred, crimes.test)## crimes.test
## BostKNN.pred 0 1
## 0 39 32
## 1 51 131mean(BostKNN.pred != crimes.test)## [1] 0.3280632KNN K=5
KNtrain <- cbind(dis, nox, age, medv)[Ctrain, ]
KNtest <- cbind(dis, nox, age, medv)[Ctest, ]
BostKNN.pred <- knn(KNtrain, KNtest, crimes.test, k=5)
table(BostKNN.pred, crimes.test)## crimes.test
## BostKNN.pred 0 1
## 0 51 28
## 1 39 135mean(Bost.pred != crimes.test)## [1] 0.1027668