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.
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.
Explore the data graphically in order to investigate the associ- ation between mpg01 and the other features. Which of the other features seem most likely to be useful in predicting mpg01? Scat- terplots and boxplots may be useful tools to answer this ques- tion. Describe your findings.
Split the data into a training set and a test set.
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?
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?
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?
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(ISLR)
library(MASS)
library(class) # knn
Column names
names(Auto)
## [1] "mpg" "cylinders" "displacement" "horsepower" "weight"
## [6] "acceleration" "year" "origin" "name"
Creating a mpg01 predictor that returns 1 if greater than median of mpg
Creating a vector that returns 1 if >= median(Auto$mpg) else 0
mpg01 = rep(0, nrow(Auto))
mpg01[Auto$mpg >= median(Auto$mpg)] = 1
new.auto = data.frame(Auto, mpg01)
Graphically viewing the correlation of mpg01 with other predictors pairs(new.auto)
Comments: There seems to be some association with horsepower, weight, acceleration, and possibly displacement. There is a distinction between values when mpg01 = 1 or mpg01 = 0.
mpg01 by displacement, acceleration, horsepower, weight
boxplot(displacement~mpg01,data=new.auto, xlab="MPG > Median", ylab="Displacement")
boxplot(acceleration~mpg01,data=new.auto, xlab="MPG > Median", ylab="Acceleration")
boxplot(horsepower~mpg01,data=new.auto, xlab="MPG > Median", ylab="Horsepower")
boxplot(weight~mpg01,data=new.auto, xlab="MPG > Median", ylab="Weight")
Comments: Based on the boxplots, acceleration might be hard to predict if mpg is higher than median. The values of acceleration are too close to tell apart.
Creating sample size for train data (80% of the data)
set.seed(1)
sample.size = floor(0.8*nrow(new.auto))
# Get random sample of indices for train data
train = sample(seq_len(nrow(new.auto)), size=sample.size)
# Creating the train and test data
train.data = new.auto[train,]
test.data = new.auto[-train, ]
# Creating responses for train and test
y.train = train.data$mpg01
y.test = test.data$mpg01
lda.fit = lda(mpg01 ~ displacement + horsepower + weight, data=new.auto, subset = train)
lda.pred = predict(lda.fit, test.data)
lda.class = lda.pred$class
table(lda.class, y.test)
## y.test
## lda.class 0 1
## 0 33 1
## 1 9 36
mean(lda.class != y.test) # test error rate: 12.7%
## [1] 0.1265823
qda.fit = qda(mpg01 ~ displacement + horsepower + weight, data=new.auto, subset = train)
qda.pred = predict(qda.fit, test.data)
qda.class = qda.pred$class
table(qda.class, y.test)
## y.test
## qda.class 0 1
## 0 36 2
## 1 6 35
mean(qda.class != y.test) # test-error rate: 10.1%
## [1] 0.1012658
lr.fit = glm(mpg01 ~ displacement + horsepower + weight,
data=new.auto,
subset = train,
family=binomial)
lr.fit.probs = predict(lr.fit, test.data, type="response")
lr.fit.pred = rep(0, nrow(test.data))
lr.fit.pred[lr.fit.probs > .5] = 1
table(lr.fit.pred, y.test)
## y.test
## lr.fit.pred 0 1
## 0 38 1
## 1 4 36
mean(lr.fit.pred != y.test) # test-error rate: 6.33%
## [1] 0.06329114
Extracting the predictors from training data
train.X = cbind(train.data$displacement + train.data$horsepower + train.data$weight)
train.y = train.data$mpg01
test.X = cbind(test.data$displacement + test.data$horsepower + test.data$weight)
scaled.train.X = scale(train.X)
scaled.test.X = scale(test.X)
knn.pred.1 = knn(scaled.train.X, scaled.test.X, train.y, k=1)
mean(y.test != knn.pred.1) # 15.1% (changed from 13.9% after scaling)
## [1] 0.1518987
knn.pred.3 = knn(scaled.train.X, scaled.test.X, train.y, k=3)
mean(y.test != knn.pred.3) # 12.7% (changed from 13.9% after scaling)
## [1] 0.1265823
knn.pred.5 = knn(scaled.train.X, scaled.test.X, train.y, k=5)
mean(y.test != knn.pred.5) # 10.1%
## [1] 0.1012658
knn.pred.10 = knn(scaled.train.X, scaled.test.X, train.y, k=10)
mean(y.test != knn.pred.10) # 10.1%
## [1] 0.1392405