week9_assi5_sol1_dhar.r

# # -*- coding: utf-8 -*-
# """Week9 Assi5 Sol1.ipynb
# 
# Automatically generated by Colaboratory.
# 
# Original file is located at
#     https://colab.research.google.com/drive/1_Et_v9mmzBrmbE8jR1XTGHR_m1d5DieY
# """
# 
###########################################################################
## Week-9, Homework-5, Sol-2
## Sreya Dhar 
## Created: Nov 04, 2020
## Edited: Nov 16, 2020
###########################################################################

# ## installing all the libaries in R kernel
# 
# install.packages("corrplot")
# install.packages("forecast")
# install.packages("zoo")
# install.packages("rsample")
# install.packages("leaps")
# install.packages("car")
# install.packages("caret")
# install.packages("ROCR")
# install.packages("PerformanceAnalytics")
# install.packages("funModeling")
# install.packages("hrbrthemes")
# install.packages("ggthemes")
# install.packages("GGally")
# install.packages("glmnet")
# install.packages("ISLR")
# install.packages("kableExtra")
# install.packages("broom")
# install.packages("knitr")
# install.packages("psych")
# install.packages("aod")
# install.packages("epiDisplay")
# install.packages("e1071")
# install.packages("class")
# install.packages('tree')
# install.packages('rpart')
# install.packages('rattle')
# install.packages("partykit")

## importing the libraries in R kernel

library(ggplot2)
library(dplyr)
library(tidyverse)
library(tidyr)
library(corrplot)
library(repr)
library(reshape2)
library(forecast)
library(zoo)
library(rsample)
library(gplots)
library(ROCR)
library(class)
library(readr)
library(leaps)
library(car)
library(PerformanceAnalytics)
library(funModeling)
library(caret)
library(MASS)
library(Hmisc)
library(hrbrthemes)
library(GGally)
library(glmnet)
library(pROC)
library(ISLR)
library(psych)
library(aod)
library(epiDisplay)
library(e1071)
library(ggthemes)
library(kableExtra)
library(broom)
library(knitr)
library(tree)
library(rpart)
library(rattle)
library(partykit)

rm(list=ls())
## set directory ##
setwd("C:/File E/EAS 506 Statistical Mining I/Week 9/Assignment-5")

## 
# setwd("C:/File E/EAS 506 Statistical Mining I/Week 9/Assignment-5")
# load("C:/File E/EAS 506 Statistical Mining I/Week 9/Assignment-5/vehicle.RData")
# write.csv(vehicle,'vehicle.csv')

vehicle <- read.csv("vehicle.csv", header = TRUE)
# vehicle

glimpse(vehicle)

head(vehicle)

names(vehicle)

glimpse(vehicle)

status(vehicle)

profiling_num(vehicle)

status(vehicle)

describe(vehicle)

summary(vehicle)

data <- vehicle %>% mutate_if(is.character, as.factor)
data_1 <- data %>% mutate_if(is.factor, as.numeric)
data_1 <- data_1[,-19] ## removing the "classdigit" variable
vehic <- data[,-19] ## removing the "classdigit" variable
glimpse(vehic)

glimpse(data_1)

options(repr.plot.width=6, repr.plot.height=6, repr.plot.res = 200)
L <- cor(data_1)
corrplot(L, method = "circle",  type = "lower")

options(repr.plot.width=6, repr.plot.height=6, repr.plot.res = 200)
plot_num(vehic)

## plotting the correlation values on chart matrix which also combined with histogram and scatter plots of different features on scaled dataset.
options(repr.plot.width=10, repr.plot.height=10, repr.plot.res = 200)
chart.Correlation(data_1, histogram=TRUE, pch=15)

options(repr.plot.width=8, repr.plot.height=8, repr.plot.res = 230)

pairs.panels(data[,1:18], main = "Pairs plot on Vehicle dataset, unclassed on Class", pch = 21, bg = c("blue", "green", "pink", "yellow")[unclass(data$class)], hist.col="red")

# Quick display of two cabapilities of GGally, to assess the distribution and correlation of variables 
 
# ggpair plot divided into three groups
options(repr.plot.width=12, repr.plot.height=12, repr.plot.res = 200)
ggpairs(data, columns = 1:18, ggplot2::aes(colour=class))+theme_bw()

# ## min-max scaling on vehicle dataset predictors prior to classification
normalized<-function(y) {
  
  x<-y[!is.na(y)]
  
  x<-(x - min(x)) / (max(x) - min(x))
  
  y[!is.na(y)]<-x
  
  return(y)
}

data_1[,1:18]<-apply(data_1[,c(1:18)],2,normalized)
vehic[,1:18]<-apply(vehic[,c(1:18)],2,normalized)
head(data_1)

## splitting the dataset into train and test sets
set.seed(1236) ## seeding the sampling
data_split <- initial_split(vehic, prop = 0.75) ## spliting the data by library 'rsample'
data_train <- training(data_split)
data_test  <- testing(data_split)

dim(data_train)

############### Classification from Tree library ################
#################################################################

class_tree = tree(class~., data = data_train, split = "deviance")
summary(class_tree)

## description of classification tree
class_tree

options(repr.plot.width=20, repr.plot.height=10, repr.plot.res = 200)
plot(class_tree)
text(class_tree ,pretty =0)

# Distributional prediction on train set
tree_predict_tr <- predict(class_tree, data_train, type='class') # gives the probability for each class
conf_tree_tr <- confusionMatrix(tree_predict_tr , data_train$class)
conf_tree_tr

# Distributional prediction on test set
tree_predict <- predict(class_tree, data_test, type='class') # gives the probability for each class
conf_tree <- confusionMatrix(tree_predict , data_test$class)
conf_tree

# prune the tree
prune_tree<- prune.tree(class_tree, best=5)
options(repr.plot.width=7, repr.plot.height=7, repr.plot.res = 200)
plot(prune_tree)
text(prune_tree)

### CV from tree model


cv_model <- cv.tree(class_tree,  FUN=prune.tree)

options(repr.plot.width=4, repr.plot.height=4, repr.plot.res = 200)
plot(cv_model, col='black', type='l')
points(cv_model$size, cv_model$dev, col='red')
# line(prune_tree$size, prune_tree$dev, type='l')
# points(prune_tree$size, prune_tree$dev, col='red')

cv_model

cv_model$dev  # gives the deviance for each K (small is better)

best_size <- cv_model$size[which(cv_model$dev==min(cv_model$dev))] # which size is better?
best_size

# let's refit the tree model (the number of leafs will be no more than best.size)
cv_model_pruned <- prune.misclass(class_tree, best=best_size)
summary(cv_model_pruned)

options(repr.plot.width=10, repr.plot.height=10, repr.plot.res = 200)
plot(cv_model_pruned)
text(cv_model_pruned)

# Distributional prediction on test set
cv_pred_pr_tr <- predict(cv_model_pruned, data_train, type='class') # gives the probability for each class
cv_tree_pr_tr <- confusionMatrix(cv_pred_pr_tr , data_train$class)
cv_tree_pr_tr

# Distributional prediction on test set
cv_pred_prune <- predict(cv_model_pruned, data_test, type='class') # gives the probability for each class
cv_tree_pr <- confusionMatrix(cv_pred_prune , data_test$class)
cv_tree_pr

########################################################################
##################### Classification Tree with rpart ###################
########################################################################

# grow a tree
model.control <- rpart.control(minsplit = 5, cp = 0,  xval = 10)
rpart_fit <- rpart(class~., method="class", data=vehic, control = model.control)
printcp(rpart_fit) # display the results

options(repr.plot.width=4, repr.plot.height=4, repr.plot.res = 200)
plotcp(rpart_fit) # visualize cross-validation results

options(repr.plot.width=4, repr.plot.height=4, repr.plot.res = 200)
plot(rpart_fit$cptable[,4], main = "Cp for model selection", ylab = "cv error", type='l', ylim=c(0.3,1.2))
points(rpart_fit$cptable[,4], col="red", cex=1,pch=20)

# summary(rpart_fit) # detailed summary of splits

#### variable importance plot from rpart #####

varimp_rpart<- data.frame( imp= rpart_fit$variable.importance)
varimp_rpart_plot <- varimp_rpart %>% 
  tibble::rownames_to_column() %>% 
  dplyr::rename("variable" = rowname) %>% 
  dplyr::arrange("Importance"=imp) %>%
  dplyr::mutate(variable = forcats::fct_inorder(variable))

options(repr.plot.width=6, repr.plot.height=6, repr.plot.res = 200)
ggplot2::ggplot(varimp_rpart_plot) +
  geom_col(aes(x = variable, y = imp),
           col = "green", show.legend = F) +
  coord_flip() +
  scale_fill_grey() +
  theme_bw()

options(repr.plot.width=5, repr.plot.height=5, repr.plot.res = 200)

ggplot2::ggplot(varimp_rpart_plot) +
  geom_segment(aes(x = variable, y = 0, xend = variable, yend = imp), 
               size = 1.5, alpha = 0.7) +
  geom_point(aes(x = variable, y = imp, col = variable), 
             size = 4, show.legend = F) +
  coord_flip() +
  theme_bw()

options(repr.plot.width=25, repr.plot.height=10, repr.plot.res = 200)
# plot tree
plot(rpart_fit, uniform=TRUE,
   main="Classification Tree for vehicle")
text(rpart_fit, use.n=TRUE, all=TRUE, cex=.8)

# create attractive postscript plot of tree
post(rpart_fit, file = "tree.ps", title = "Classification Tree for vehicle")

rpart_predict_tr <- predict(rpart_fit, data_train, type="class") # give the predicted class
# table(rpart_predict, data_test$class)
confusionMatrix(rpart_predict_tr, data_train$class)

rpart_predict <- predict(rpart_fit, data_test, type="class") # give the predicted class
# table(rpart_predict, data_test$class)
confusionMatrix(rpart_predict, data_test$class)

# prune the tree
min_cp = which.min(rpart_fit$cptable[,4])
prune_fit<- prune(rpart_fit, cp=rpart_fit$cptable[min_cp,1])
options(repr.plot.width=22, repr.plot.height=8, repr.plot.res = 200)
# plot the pruned tree
plot(prune_fit, uniform=TRUE,
   main="Pruned Classification Tree for Vehicle")
text(prune_fit, use.n=TRUE, all=TRUE, cex=.8)
post(prune_fit, file = "pruned tree.ps",
   title = "Pruned Classification Tree for Vehicle")

summary(prune_fit)

prune_fit$variable.importance

#### variable importance plot from pruned rpart #####

varimp_prpart<- data.frame( imp= prune_fit$variable.importance)
varimp_prpart_plot <- varimp_prpart %>% 
  tibble::rownames_to_column() %>% 
  dplyr::rename("variable" = rowname) %>% 
  dplyr::arrange("Importance"=imp) %>%
  dplyr::mutate(variable = forcats::fct_inorder(variable))

options(repr.plot.width=5, repr.plot.height=5, repr.plot.res = 200)

ggplot2::ggplot(varimp_prpart_plot) +
  geom_segment(aes(x = variable, y = 0, xend = variable, yend = imp), 
               size = 1.5, alpha = 0.7) +
  geom_point(aes(x = variable, y = imp, col = variable), 
             size = 4, show.legend = F) +
  coord_flip() +
  theme_bw()

pruned_predict_tr <- predict(prune_fit, data_train, type="class") # give the predicted class
# table(pruned_predict, data_test$class)
confusionMatrix(pruned_predict_tr, data_train$class)

pruned_predict <- predict(prune_fit, data_test, type="class") # give the predicted class
# table(pruned_predict, data_test$class)
confusionMatrix(pruned_predict, data_test$class)

### Plot for rparty tree at least of 5 cp values #### 
options(repr.plot.width=12, repr.plot.height=10, repr.plot.res = 200)
min_cp = which.min(rpart_fit$cptable[,4])
rpart_prune5 <- prune(rpart_fit, cp = rpart_fit$cptable[min_cp-5,1])
plot(rpart_prune5, uniform=TRUE, branch=0.6, margin=0.05)
text(rpart_prune5, all=TRUE, use.n=TRUE)
title("Training Set's Classification Tree")

summary(rpart_prune5)

options(repr.plot.width=12, repr.plot.height=8, repr.plot.res = 230)
# plot decision tree
fancyRpartPlot(rpart_prune5, main="class")

rparty_tree <- as.party(rpart_prune5)
rparty_tree

options(repr.plot.width=20, repr.plot.height=6, repr.plot.res = 230)
plot(rparty_tree)

### end ###