week12_assi6_sol2_Dhar.R

# # -*- coding: utf-8 -*-
# """Week12 Assi6 Sol2.ipynb
# 
# Automatically generated by Colaboratory.
# 
# Original file is located at
#     https://colab.research.google.com/drive/1q-y9n5pJ16jmsw-7mCrxsKYUma6s3djS
# """
###########################################################################
## Week-12, Homework-6, Sol-2 
## Sreya Dhar 
## Created: Nov 24, 2020
## Edited: Dec 3, 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")
# install.packages("randomForest")
# install.packages("party")
# install.packages("kernlab")
# install.packages("gbm")
# install.packages("plotmo")
# """# New Section"""

## 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(rpart)
library(rattle)
library(partykit)
library(randomForest)
library(tree)
library(party)
library(kernlab)
library(gbm)
library(plotmo)
library(tree)


## set directory ##
setwd("C:/File G/EAS 506 Statistical Mining I/Week 12/Assignment-6")
load("C:/File G/EAS 506 Statistical Mining I/Week 12/Assignment-6/pima.RData")
write.csv(pima,'pima.csv')

## uploading the data
pima = read.csv('pima.csv',header = TRUE)

head(pima)

pima <- pima[,-8]
head(pima)

glimpse(pima)

names(pima)

status(pima)

profiling_num(pima)

describe(pima)

summary(pima)

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

pairs.panels(pima[,1:7], main = "Pairs plot on pima dataset, unclassed on class", pch = 21, bg = c("blue", "green")[unclass(pima$class)], hist.col="red")

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

pima <- pima %>% mutate_if(is.character, as.factor)
pima_n <- pima %>% mutate_if(is.factor, as.numeric)

# heatmap matrix 
options(repr.plot.width=5, repr.plot.height=5, repr.plot.res = 230)
data_h <- as.data.frame(scale(pima_n,center=TRUE,scale=TRUE))
heatmap.2(as.matrix(pima_n), scale = "none", col = bluered(100), trace = "none", density.info = "none")

## 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(pima_n, histogram=TRUE, pch=15)

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


########################################################################
##################### Classification Tree with rpart ###################
########################################################################
# Normalizing the predictor variable 

# ## 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)
}

pima[,1:7]<-apply(pima[,c(1:7)],2,normalized)
head(pima)

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

############### Classification from Tree library ################
#################################################################
tree_model <- tree(class ~ ., data=data_train, split = "deviance")
# tree_model

summary(tree_model)

options(repr.plot.width=17, repr.plot.height=6, repr.plot.res = 200)
plot(tree_model)
text(tree_model)

# Distributional prediction on train set
tree_pred_tr <- predict(tree_model, data_train, type="class") # gives the probability for each class

confusionMatrix(table(tree_pred_tr, data_train$class))

# Distributional prediction on test set
tree_pred <- predict(tree_model, data_test, type="class") # gives the probability for each class

confusionMatrix(table(tree_pred, data_test$class))

options(repr.plot.width=5, repr.plot.height=5, repr.plot.res = 200)
cv_model <- cv.tree(tree_model, FUN = prune.misclass)
plot(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(tree_model, best=best_size)
summary(cv_model_pruned)

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

pruned_tr_cv <- predict(cv_model_pruned, data_train, type="class") # give the predicted class
confusionMatrix(table(Prediction=pruned_tr_cv, Reference=data_train$class))

pruned_pred_cv <- predict(cv_model_pruned, data_test, type="class") # give the predicted class
confusionMatrix(table(Prediction=pruned_pred_cv, Reference=data_test$class))




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

options(repr.plot.width=30, repr.plot.height=10, repr.plot.res = 200)
model.control <- rpart.control(minsplit = 5, cp = 0, xval = 10)
rpart_tree <- rpart(class ~ ., data=data_train, control =model.control)
plot(rpart_tree, uniform=TRUE, branch=0.6, margin=0.05)
text(rpart_tree, all=TRUE, use.n=TRUE)
title("Classification Tree on Train set")

# summary(rpart_tree )

options(repr.plot.width=20, repr.plot.height=10, repr.plot.res = 200)
# plot decision tree
fancyRpartPlot(rpart_tree, main="pima")

## rpart predition for train set
rpart_pred_tr <- predict(rpart_tree, data_train, type="class")
confusionMatrix(table(data_train$class, rpart_pred_tr))

rpart_pred <- predict(rpart_tree, data_test, type="class")
confusionMatrix(table(data_test$class, rpart_pred))

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

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

min_cp = which.min(rpart_tree$cptable[,4])
prune_rpart<- prune(rpart_tree, cp=rpart_tree$cptable[min_cp, 1])
 # pruning the tree
options(repr.plot.width=10, repr.plot.height=10, repr.plot.res = 200)
plot(prune_rpart, uniform=TRUE, branch=0.6)
text(prune_rpart, all=TRUE, use.n=TRUE)

# summary(prune_rpart)

rpart_pr_pred_tr <- predict(prune_rpart, data_train, type="class")
confusionMatrix(table(data_train$class, rpart_pr_pred_tr))

rpart_pr_pred_te <- predict(prune_rpart, data_test, type="class")
confusionMatrix(table(data_test$class, rpart_pr_pred_te))


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

varimp_rpart<- data.frame( imp= rpart_tree$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=4, 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() + labs(title = "Variable Importance from rpart tree")
  theme_bw()

# rparty_tree <- as.party(rpart_tree)
# rparty_tree

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

########################################
### Conditional Inference Tree ###
########################################
options(repr.plot.width=10, repr.plot.height=6, repr.plot.res = 200)
inf_model <- ctree(class ~ . , data = data_train)
plot(inf_model)

inf_pred_tr<- predict(inf_model, data_train[,-8])
inf_tab_tr <- table(Prediction=inf_pred_tr, Reference=data_train$class)
confusionMatrix(inf_tab_tr)

inf_pred<- predict(inf_model, data_test[,-8])
inf_tab <- table(Prediction=inf_pred, Reference=data_test$class)
confusionMatrix(inf_tab)

# get the probabilities from the barplots showen above:
tapply(treeresponse(inf_model), where(inf_model), unique)

# The package is able to format the plot tree. Eg:
innerWeights <- function(node){
  grid.circle(gp = gpar(fill = "White", col = 1))
  mainlab <- paste( node$psplit$variableName, "\n(n = ")
  mainlab <- paste(mainlab, sum(node$weights),")" , sep = "")
  grid.text(mainlab,gp = gpar(col='red'))
}

options(repr.plot.width=15, repr.plot.height=5, repr.plot.res = 200)
plot(inf_model, type='simple', inner_panel = innerWeights)

##############################################################
#################### Random Forest ###########################
##############################################################
rf <- randomForest(class ~., data=data_train, importance=TRUE, do.trace=100, ntree=500, mtry = 3)
print(rf)

rf$err.rate[500]

options(repr.plot.width=10, repr.plot.height=5, repr.plot.res = 200)
varImpPlot(rf, main="Variable Importance Accuracy and Gini coeff. from Random Forest")

importance(rf)

## Prediction on train set
rf_pred_tr <- predict(rf, data_train, type = "response")
rf_tab <- confusionMatrix(table(data_train$class, rf_pred_tr))
rf_tab

## Prediction on test set
rf_predict <- predict(rf, data_test,type = "response")
rf_tab <- confusionMatrix(table(data_test$class, rf_predict))
rf_tab

# options(repr.plot.width=4, repr.plot.height=4, repr.plot.res = 200)
# # next function gives a graphical depiction of the marginal effect of a variable on the class probability (classification) 
# partialPlot(rf, data_train, all, "normal")

rf_tree <- getTree(rf, k=2) # show the second tree
# print(rf_tree)

treesize(rf) # size of trees of the ensemble

mean(treesize(rf))

options(repr.plot.width=4, repr.plot.height=4, repr.plot.res = 200)
hist(treesize (rf), col = "red")

######################################### Iteration on mtry for RF #########################################

rf.c <- list()
oob.err <- list()
yhat.rf <-list()
misclass_rf<-list()
for ( i in 1:7 ) {
  set.seed(4444)
  rf.c<-randomForest(class ~ ., data = data_train, mtry = i, importance = TRUE, ntree = 500)
  oob.err[i] = rf.c$err.rate[500] #Error of all Trees fitted
  yhat.rf<-predict(rf.c, newdata = data_test[,-8])
  misclass_rf[i]<- mean(yhat.rf != data_test$class)
}

options(repr.plot.width=6, repr.plot.height=4, repr.plot.res = 200)
matplot(1:7, misclass_rf, xlab = 'No. of Variables (mtry)', ylab = 'Misclassification error', main = "Subset size Vs. Misclassification error")
lines(1:7, misclass_rf, type = "o")

options(repr.plot.width=7, repr.plot.height=7, repr.plot.res = 200)
matplot(1:7 , cbind(oob.err, misclass_rf), pch=19 , col=c("red","blue"), type=c("o","o"), ylab="Misclassification Error rate", xlab="Number of Predictors Considered at each Split, mtry")
legend("topright", legend=c("Out of Bag Error","Test Error"), pch=19, col=c("red","blue"))

misclass_rf ## min for mtry=4,7 on test set

################ Tuning on RF model ###############

library("e1071") # to 'tune' the rf model

tuned_rf <- tune(randomForest, train.x = class ~ ., data = data_train, validation.x = data_test, mtry=4, ntree=500)

best_rf <- tuned_rf$best.model
best_rf_pred <- predict(best_rf, data_test, type='response')
tab_best_rf <- table(Reference=data_test$class, Prediction=best_rf_pred)
confusionMatrix(tab_best_rf)

tuned_rf

best_rf

options(repr.plot.width=10, repr.plot.height=5, repr.plot.res = 200)
par(mfrow=c(1,2))
plot(rf)
legend("topright", 150, legend=c("OOB Error", "diabetic oob error", "normal oob error"), col=c("black", "green", "red"), lty=1:2, cex=1.0)
plot(best_rf)

# computing overall error:
err_rf <- 1 - sum(diag(as.matrix(tab_best_rf))) / sum(tab_best_rf)
err_rf

install.packages("pdp")
library(pdp)

# Two Variables partial dependence plot
options(repr.plot.width=7, repr.plot.height=5, repr.plot.res = 250)
par.glucose.age <- partial(rf, pred.var = c("glucose", "age"), chull = TRUE)
plot.glucose.age <- autoplot(par.glucose.age, contour = TRUE, 
               legend.title = "Partial\ndependence")

grid.arrange(plot.glucose.age )

head(data_train)

options(repr.plot.width=7, repr.plot.height=5, repr.plot.res = 250)
par.bmi.npreg <- partial(rf, pred.var = c("bmi", "npregnant"), chull = TRUE)
plot.bmi.npreg <- autoplot(par.bmi.npreg, contour = TRUE, 
               legend.title = "Partial\ndependence")

grid.arrange(plot.bmi.npreg )

options(repr.plot.width=7, repr.plot.height=5, repr.plot.res = 250)
par.gluc.npreg <- partial(rf, pred.var = c("glucose", "npregnant"), chull = TRUE)
plot.gluc.npreg <- autoplot(par.gluc.npreg, contour = TRUE, 
               legend.title = "Partial\ndependence")

grid.arrange(plot.gluc.npreg )

options(repr.plot.width=7, repr.plot.height=5, repr.plot.res = 250)
par.ped.age <- partial(rf, pred.var = c("bmi", "age"), chull = TRUE)
plot.ped.age <- autoplot(par.ped.age, contour = TRUE, 
               legend.title = "Partial\ndependence")

grid.arrange(plot.ped.age )

options(repr.plot.width=6, repr.plot.height=6, repr.plot.res = 230)
par(mfrow=c(2,2))
# next function gives a graphical depiction of the marginal effect of a variable on the class probability (classification) 
partialPlot(rf, data_train, glucose, "diabetic")
partialPlot(rf, data_train, age, "diabetic")
partialPlot(rf, data_train, bmi, "diabetic")
partialPlot(rf, data_train, npregnant, "diabetic")

plotmo(rf,  pmethod="partdep", all1=TRUE, all2=TRUE)

options(repr.plot.width=3, repr.plot.height=3, repr.plot.res = 230)
gg1 <- rf %>%  # the %>% operator is read as "and then"
  partial(pred.var = "glucose") %>%
  autoplot(smooth = TRUE, ylab = expression(f(glucose))) +
  theme_light() +
  ggtitle("ggplot2-based PDP for Glucose")
gg1

options(repr.plot.width=3, repr.plot.height=3, repr.plot.res = 230)
gg2 <- rf %>%  # the %>% operator is read as "and then"
  partial(pred.var = "bmi") %>%
  autoplot(smooth = TRUE, ylab = expression(f(bmi))) +
  theme_light() +
  ggtitle("ggplot2-based PDP for BMI")
gg2

data_train_n <- data_train %>% mutate_if(is.factor, as.numeric)-1
data_test_n <- data_test %>% mutate_if(is.factor, as.numeric)-1

###############################################
# Boosting ##
###############################################
mod_boost = gbm(class ~.,
              data = data_train_n,
              distribution = "adaboost",
              cv.folds = 10,
              shrinkage = .1,
              n.minobsinnode = 10,
              n.trees = 200)

print(mod_boost)

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

mod_boost1 = gbm(class ~.,
              data = data_train,
              distribution = "multinomial",
              cv.folds = 10,
              shrinkage = .1,
              n.minobsinnode = 10,
              n.trees = 200)

print(mod_boost1)

### prediction on train set from boosting ###
pred_boost_tr = predict.gbm(object = mod_boost1,
                            newdata = data_train,
                            n.trees = 200,
                            type = "response")

labels_tr = colnames(pred_boost_tr)[apply(pred_boost_tr, 1, which.max)]
result_tr = data.frame(data_train$class, labels_tr)
confusionMatrix(data_train$class, as.factor(labels_tr))

### prediction on test set from boosting ###
pred_boost_te = predict.gbm(object = mod_boost1,
                   newdata = data_test,
                   n.trees = 200,
                   type = "response")


labels = colnames(pred_boost_te)[apply(pred_boost_te, 1, which.max)]
result = data.frame(data_test$class, labels)
# print(result)

confusionMatrix(data_test$class, as.factor(labels))

#### Shrinkage factor iteration for boosting ###
shrink <- c(.1, .4, .6, .8)
max_iter <- 5000
store_error <- c()
for (i in 1:length(shrink)){
	boost.fit <- gbm(class~., data = data_train_n, n.trees = max_iter, shrinkage = shrink[i], interaction.depth = 3, distribution = "adaboost")
	temp <- c()
	for (j in 1:max_iter){
		y_hat <- predict(boost.fit, newdat = data_test_n, n.trees = j, type = "response")
		misclass_boost <- mean(y_hat != data_test_n$class)
		temp <- c(temp, misclass_boost)
	}
	store_error <- cbind(store_error, temp) # max_iter x length(shrink)
}

colnames(store_error) <- paste("shrinkage", shrink, sep = ":")

options(repr.plot.width=7, repr.plot.height=7, repr.plot.res = 200)
plot(store_error[,1], type = "l", main = "Error Profiles", ylab = "error", xlab = "boosting iterations", ylim=c(0.4,1))
lines(store_error[,2], col = "red")
lines(store_error[,3], col = "blue")
lines(store_error[,4], col = "green")
legend("bottomleft", legend=c("shrinkage = 0.1","shrinkage = 0.4", "shrinkage = 0.6", "shrinkage = 0.8"), pch=19, col=c('black',"red","blue", "green"))

options(repr.plot.width=8, repr.plot.height=8, repr.plot.res = 230)
plotmo(mod_boost,  pmethod="partdep", all1=TRUE, all2=TRUE)

### end ###