In this section, I compare the common machine learning models by solving the question 11 at chapter 4 from the An Introduction to Statistical Learning Book.
Data:
Auto data set: A data frame with 392 observations on the following 9 variables.
mpg:miles per gallon
cylinders: Number of cylinders between 4 and 8
displacement:Engine displacement (cu. inches)
horsepower:Engine horsepower
weight:Vehicle weight (lbs.)
acceleration:Time to accelerate from 0 to 60 mph (sec.)
year:Model year (modulo 100)
origin:Origin of car (1. American, 2. European, 3. Japanese)
name:Vehicle name
The orginal data contained 408 observations but 16 observations with missing values were removed.
Objective:
In this problem, I will develop a ML models to predict whether a given car gets high or low gas mileage based on the Auto data set.
Loading Librarires
#Loading necessary libraries
library(tidyverse) # For data manipulation
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.4 v dplyr 1.0.7
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 2.0.1 v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(dplyr)
library(caret) # confusion matrix
## Loading required package: lattice
##
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
##
## lift
library(GGally)
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
library(ISLR)
library(ggplot2)
library(corrplot)
## corrplot 0.90 loaded
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
library(class)
library(bootStepAIC)
#PREPARING WORK SPAcE
# Clear the workspace:
rm(list = ls())
Data Preparation/Exploration
# Load data
names(Auto)
## [1] "mpg" "cylinders" "displacement" "horsepower" "weight"
## [6] "acceleration" "year" "origin" "name"
dim(Auto)
## [1] 392 9
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
head(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 500
(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.
Auto$mpg01 <- rep(0,392)
Auto$mpg01[Auto$mpg>median(Auto$mpg)]=1
head(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 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
#Convert to Categorical Variable
df <- Auto %>%
mutate_at(vars(mpg01), funs(factor))
str(df)
## 'data.frame': 392 obs. of 10 variables:
## $ mpg : num 18 15 18 16 17 15 14 14 14 15 ...
## $ cylinders : num 8 8 8 8 8 8 8 8 8 8 ...
## $ displacement: num 307 350 318 304 302 429 454 440 455 390 ...
## $ horsepower : num 130 165 150 150 140 198 220 215 225 190 ...
## $ weight : num 3504 3693 3436 3433 3449 ...
## $ acceleration: num 12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
## $ year : num 70 70 70 70 70 70 70 70 70 70 ...
## $ origin : num 1 1 1 1 1 1 1 1 1 1 ...
## $ name : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...
## $ mpg01 : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
(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? Scatter plots and boxplots may be useful tools to answer this question.Describe your findings.
#Correlation Matrix
corrplot(cor(Auto[,-9]), method="square")
We can notice that Cylinder, displacement, horsepower and weight are highly negatively correlated.
#Paris Plot
ggpairs(df, # Data frame
columns = 1:8, # Columns
aes(color = mpg01, # Color by group (cat. variable)
alpha = 0.5), # Transparency
upper = list(continuous = wrap("cor", size = 2.4))) # Font size
Green dots represents the MPG data lower than median, and Red dots are the ones higher than median. In most graph, red and green dots are separated already, it is proof of the correlation. Only year and acceleration plots are mixed which is a sign no-correlation.
#Scatter plot
ggplot(data = df) +
geom_point(mapping = aes(x = cylinders, y = horsepower,color = mpg01))
#Box Plot
ggplot(df, aes(mpg01,horsepower, fill=mpg01)) +
geom_boxplot() +
facet_grid(~cylinders) + scale_fill_brewer(palette = "Set1")
You can see the relationship between horsepower and cylinders for the vehicles with MPG consumption. It is hard to make comment about MPG consumption for the vehicles that has 6 cylinders with known horsepower. But it is easy to make comments for the vehicles with other type of cylinders, because the colors are already separated in the plot.
(c) Split the data into a training set and a test set.
ind <- sample(2, nrow(Auto), replace=T, prob =c(0.8, 0.2))
train <- Auto[ind==1,]
test <- Auto[ind==2,]
Mpg01_test = test$mpg01
(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?
set.seed(111)
#Fit Model
lda.fit = lda(mpg01~ cylinders + displacement + horsepower + weight + origin, data=train)
#Predict
lda.pred <-predict(lda.fit, test)
#Assign predicted classes
lda.class <-lda.pred$class
#Create a Confusion Matrix
table(lda.class, Mpg01_test)
## Mpg01_test
## lda.class 0 1
## 0 40 3
## 1 1 42
#Accuracy Rate
(acc<-mean(lda.class==Mpg01_test))
## [1] 0.9534884
#Test Error:
1-acc
## [1] 0.04651163
LDA method predicted test data set with an error rate of 10.8%
(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?
set.seed(111)
#Fit Model
qda.fit <- qda(mpg01~ cylinders + displacement + horsepower + weight + origin, data=train)
#Predict
qda <- predict(qda.fit, test)
#Assign predicted classes
qda.class <- qda$class
#Create a Confusion Matrix
table(qda.class, Mpg01_test)
## Mpg01_test
## qda.class 0 1
## 0 40 4
## 1 1 41
#Accuracy Rate
(acc<-mean(qda.class==Mpg01_test))
## [1] 0.9418605
#Test Error:
1-acc
## [1] 0.05813953
QDA method predicted test data set with an error rate of 8.1%
(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?
set.seed(111)
#Fit Model
glm.fits <- glm(mpg01~ cylinders + displacement + horsepower + weight + origin, data=train, family=binomial)
#Predict
glm.probs = predict(glm.fits, test, type = "response")
#Assign predicted classes
glm.pred <- rep(0,dim(test)[1])
glm.pred[glm.probs>.5]=1
#Create a Confusion Matrix
table(glm.pred,Mpg01_test)
## Mpg01_test
## glm.pred 0 1
## 0 41 5
## 1 0 40
#Accuracy Rate
(acc<-mean(glm.pred==Mpg01_test))
## [1] 0.9418605
#Test Error:
1-acc
## [1] 0.05813953
Logistic Regression method predicted test data set with an error rate of 12.16%
(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?
#Data Partition for KNN Model
nr<- nrow(Auto)
trainingindex <- sample(1:nr, round(.8*nr) )
train <- rep(FALSE,nr)
train [trainingindex]=TRUE
test <- Auto[!train,]
Mpg01_data = Auto$mpg01[!train]
train.KNN <- cbind(Auto$cylinders,Auto$displacement,Auto$horsepower,Auto$weight,Auto$origin)[train,]
test.KNN <- cbind(Auto$cylinders,Auto$displacement,Auto$horsepower,Auto$weight,Auto$origin)[!train,]
train.mpg01 =Auto$mpg01[train]
set.seed(111)
#K=2
knn.pred <- knn(train.KNN,test.KNN, train.mpg01, k=1)
#Error Rate
mean(knn.pred!=Mpg01_data)
## [1] 0.1282051
#K=10
knn.pred <- knn(train.KNN,test.KNN, train.mpg01, k=10)
#Error Rate
mean(knn.pred!=Mpg01_data)
## [1] 0.1410256
#K=20
knn.pred <- knn(train.KNN,test.KNN, train.mpg01, k=20)
#Error Rate
mean(knn.pred!=Mpg01_data)
## [1] 0.1410256
#K=40
knn.pred <- knn(train.KNN,test.KNN, train.mpg01, k=40)
#Error Rate
mean(knn.pred!=Mpg01_data)
## [1] 0.1538462
#K=50
knn.pred <- knn(train.KNN,test.KNN, train.mpg01, k=50)
#Error Rate
mean(knn.pred!=Mpg01_data)
## [1] 0.1538462
KNN method with K=2 predicted test data set with a test error rate of 12.82%. With these results, Quadratic Discriminant Analysis returned the model with a lowest test error rate with 8%.