library(hydroGOF)
library(e1071)
library(caret)
library(dplyr)
library(ggplot2)
library(glmnet)
library(kernlab)
library(rpart)
library(randomForest)
library(gbm)
install.packages("vegan")
library(vegan)

#Creating heat map 
data(mtcars)
mtcars_std<-as.data.frame(scale(mtcars))
mtcars_matrix<-data.matrix(mtcars_std)
heatmap(mtcars_matrix,col=rainbow(256))

#Principal component analysis
data(mtcars)
mtcars_std<-as.data.frame(scale(mtcars))
pca1<-prcomp(mtcars_std)
pca1$sdev/sum(pca1$sdev)
loads<-pca1$rotation
scores<-pca1$x
#Select number of PCs
plot(pca1$sdev)
(pca1$sdev[1]+pca1$sdev[2]+pca1$sdev[3]+pca1$sdev[4])/sum(pca1$sdev)

#VIP Calculation
#update loadings w/ reduced no. of PCs
loads_vip<-loads[,1:4]
property_vip<-loads_vip[1,]
features_vip<-loads_vip[2:11,]
weight_vip<-property_vip*features_vip
#no. of weights should be equal to number of PCs included
vip<-weight_vip[,1]+weight_vip[,2]+weight_vip[,3]+weight_vip[,4]
barplot(vip)

# IsoMap
data(mtcars)
mtcars_std<-as.data.frame(scale(mtcars))
mtcars_matrix<-data.matrix(mtcars_std)
dist_matrix <- dist(data_matrix)
# k = 10 specifies the number of shortest dissimilarities (nearest neighbors) to retain
mtcars_isomap <- isomap(dist_matrix, ndim = 3, k = 25)
# 'net = TRUE' visualizes the underlying neighborhood graph connections
plot(mtcars_isomap, net = TRUE, col = as.factor(mtcars$mpg), pch = 16)

#Linear Regression
cars_cor<-round(cor(mtcars),digits=2)
cars_cor
#Identify the most correlated feature and use that for LR
#For cars data, wt is most correlated (both from cor and VIP)
data_lr<-mtcars[,c("mpg","wt")]
#Divide into training and testing data
set.seed(508)
mtcars_datasort_lr<-sample(1:nrow(data_lr),nrow(data_lr)*.8)
train_lr<-data_lr[mtcars_datasort_lr,]
test_lr<-data_lr[-mtcars_datasort_lr,]
#Build the model and predict the values
mdl_lr<-lm(mpg~wt,data=train_lr)
pred_train_lr<-predict(mdl_lr,train_lr)
pred_test_lr<-predict(mdl_lr,test_lr)
#Get RMSE values
rmse_lr_train<-rmse(pred_train_lr,train_lr$mpg)
rmse_lr_train
rmse_lr_test<-rmse(pred_test_lr,test_lr$mpg)
rmse_lr_test
#R2 value for training data
sst<-sum((train_lr$mpg-mean(train_lr$mpg))^2)
sse<-sum((pred_train_lr-train_lr$mpg)^2)
rsq<-1-sse/sst
rsq
plot(train_lr$mpg,pred_train_lr,xlab="Actual",ylab="Predicted")

#Multiple Linear Regression
data(mtcars)
set.seed(508)
mtcars_datasort_mlr<-sample(1:nrow(mtcars),nrow(mtcars)*.8)
train_mlr<-mtcars[mtcars_datasort_mlr,]
test_mlr<-mtcars[-mtcars_datasort_mlr,]
mdl_mlr<-lm(mpg~.,data=train_mlr)
pred_train_mlr<-predict(mdl_mlr,train_mlr)
pred_test_mlr<-predict(mdl_mlr,test_mlr)
#Get RMSE values
rmse_mlr_train<-rmse(pred_train_mlr,train_mlr$mpg)
rmse_mlr_train
rmse_mlr_test<-rmse(pred_test_mlr,test_mlr$mpg)
rmse_mlr_test
#R2 value for training data
sst<-sum((train_mlr$mpg-mean(train_mlr$mpg))^2)
sse<-sum((pred_train_mlr-train_mlr$mpg)^2)
rsq<-1-sse/sst
rsq
plot(train_mlr$mpg,pred_train_mlr,xlab="Actual",ylab="Predicted")

# Principal Component Regression (PCR)
prop_PCR<-mtcars$mpg
feature_PCR<-mtcars[,!names(mtcars)%in%c("mpg")]
feature_PCR_std<-as.data.frame(scale(feature_PCR))
pca<-prcomp(feature_PCR_std)
plot(pca$sdev)
#Select number of PCs
#In this case, we'll select 4
scores<-pca$x[,1:4]
loads<-pca$rotation[,1:4]
pcr_data<-cbind(prop_PCR,scores)
#Have now defined new data with property and PCs
#Now repeat the MLR process
set.seed(508)
mtcars_datasort_pcr<-sample(1:nrow(pcr_data),nrow(pcr_data)*.8)
train_pcr<-pcr_data[mtcars_datasort_pcr,]
test_pcr<-pcr_data[-mtcars_datasort_pcr,]
train_pcr_d<-as.data.frame(train_pcr)
test_pcr_d<-as.data.frame(test_pcr)
mdl_pcr<-lm(prop_PCR~.,data=train_pcr_d)
pred_train_pcr<-predict(mdl_pcr,train_pcr_d)
pred_test_pcr<-predict(mdl_pcr,test_pcr_d)
#Get RMSE values
rmse_pcr_train<-rmse(pred_train_pcr,train_pcr_d$prop_PCR)
rmse_pcr_train
rmse_pcr_test<-rmse(pred_test_pcr,test_pcr_d$prop_PCR)
rmse_pcr_test
#R2 value for training data
sst<-sum((train_pcr_d$prop_PCR-mean(train_pcr_d$prop_PCR))^2)
sse<-sum((pred_train_pcr-train_pcr_d$prop_PCR)^2)
rsq<-1-sse/sst
rsq
plot(train_pcr_d$prop_PCR,pred_train_pcr,xlab="Actual",ylab="Predicted")

#Ridge Regression
data(mtcars)
set.seed(508)
datasort<-sample(1:nrow(mtcars),nrow(mtcars)*0.8)
train_ridge<-mtcars[datasort,]
test_ridge<-mtcars[-datasort,]
descriptors_train_ridge<-train_ridge[,! names(train_ridge) %in% c("mpg")]
descriptors_test_ridge<-test_ridge[,! names(test_ridge) %in% c("mpg")]
descriptors_train_ridge<-as.matrix(descriptors_train_ridge)
descriptors_test_ridge<-as.matrix(descriptors_test_ridge)
mdl_ridge<-glmnet(descriptors_train_ridge,train_ridge$mpg,alpha=0)
mdl_ridge_cv<-cv.glmnet(descriptors_train_ridge,train_ridge$mpg,alpha=0)
best_lambda<-mdl_ridge_cv$lambda.min
mdl_ridge_best<-glmnet(descriptors_train_ridge,train_ridge$mpg,alpha=0,lambda=best_lambda)
coef(mdl_ridge_best)
pred_train_ridge<-predict(mdl_ridge,s=best_lambda,newx=descriptors_train_ridge)
pred_test_ridge<-predict(mdl_ridge,s=best_lambda,newx=descriptors_test_ridge)
pred_train_ridge<-as.data.frame(pred_train_ridge)
pred_test_ridge<-as.data.frame(pred_test_ridge)
rmse_ridge_train<-rmse(pred_train_ridge,train_ridge$mpg)
rmse_ridge_test<-rmse(pred_test_ridge,test_ridge$mpg)
rmse_ridge_train
rmse_ridge_test
sst<-sum((train_ridge$mpg-mean(train_ridge$mpg))^2)
sse<-sum((pred_train_ridge-train_ridge$mpg)^2)
rsq<-1-sse/sst
rsq

#Lasso Regression
#The process is exactly the same as Ridge, we just change alpha from 0 to 1
#Lines 132 - 135

#Support vector regression
data(mtcars)
set.seed(508)
svr_datasort<-sample(1:nrow(mtcars),nrow(mtcars)*0.8)
train_svr<-mtcars[svr_datasort,]
test_svr<-mtcars[-svr_datasort,]
train_svr_d<-data.frame(train_svr)
descriptors_train_svr<-train_svr[,! names(train_svr) %in% c("mpg")]
descriptors_test_svr<-test_svr[,! names(test_svr) %in% c("mpg")]
descriptors_train_svr<-as.matrix(descriptors_train_svr)
descriptors_test_svr<-as.data.frame(descriptors_test_svr)
prop_train_svr<-train_svr$mpg
prop_test_svr<-test_svr$mpg
mdl_svr<-tune(svm,prop_train_svr~descriptors_train_svr,ranges=list(epsilon=seq(0,1,0.1),cost=1:10))
BstModel<-mdl_svr$best.model
summary(BstModel)
#Update the regression model with the selections from BstModel (kernel, cost, gamma, epsilon)
svmfit <- svm(train_svr$mpg ~., data = train_svr, method="eps-regression",kernel = 'radial', cost = 3, gamma=0.1,epsilon=.3,scale=FALSE)
pred_train_svr<-predict(svmfit, data=descriptors_train_svr)
pred_test_svr<-predict(svmfit,newdata=descriptors_test_svr)
rmse_SVR_train<-rmse(pred_train_svr,prop_train_svr)
rmse_SVR_test<-rmse(pred_test_svr,prop_test_svr)
rmse_SVR_train
rmse_SVR_test
sst<-sum((train_svr$mpg-mean(train_svr$mpg))^2)
sse<-sum((pred_train_svr-train_svr$mpg)^2)
rsq<-1-sse/sst
rsq

#Gaussian Process Regression
data(mtcars)
set.seed(508)
datasort<-sample(1:nrow(mtcars),nrow(mtcars)*0.8)
train_gpr<-mtcars[datasort,]
test_gpr<-mtcars[-datasort,]
descriptors_train_gpr<-train_gpr[,! names(train_gpr) %in% c("mpg")]
descriptors_test_gpr<-test_gpr[,! names(test_gpr) %in% c("mpg")]
mdl_gpr<-gausspr(descriptors_train_gpr,train_gpr$mpg)
pred_train_gpr<-predict(mdl_gpr,descriptors_train_gpr)
pred_test_gpr<-predict(mdl_gpr,descriptors_test_gpr)
rmse_gpr_train<-rmse(pred_train_gpr,as.matrix(train_gpr$mpg))
rmse_gpr_test<-rmse(pred_test_gpr,as.matrix(test_gpr$mpg))
rmse_gpr_train
rmse_gpr_test
sst<-sum((train_gpr$mpg-mean(train_gpr$mpg))^2)
sse<-sum((pred_train_gpr-train_gpr$mpg)^2)
rsq<-1-sse/sst
rsq


#Random Forest Regression
data(mtcars)
set.seed(508)
datasort<-sample(1:nrow(mtcars),nrow(mtcars)*0.8)
train_rf<-mtcars[datasort,]
test_rf<-mtcars[-datasort,]
model_rf<-randomForest(train_rf$mpg~.,data=train_rf,mtry=3,importance=TRUE,na.action=na.omit)
pred_train_rf<-predict(model_rf,train_rf)
pred_test_rf<-predict(model_rf,newdata=test_rf)
rmse_rf_train<-rmse(pred_train_rf,train_rf$mpg)
rmse_rf_test<-rmse(pred_test_rf,test_rf$mpg)
rmse_rf_train
rmse_rf_test
sst<-sum((train_rf$mpg-mean(train_rf$mpg))^2)
sse<-sum((pred_train_rf-train_rf$mpg)^2)
rsq<-1-sse/sst
rsq
plot(model_rf)

########################
# CATEGORICAL REGRESSIONS

#Support Vector Machine
#Categorize the Data
hist(mtcars$mpg)
#from dataset, we can pick a cutoff.  For this, we will define 22.5 as it is in the middle.  Obviously not particularly rigorous in our selection, but we'll go with it for the example.
High=ifelse(mtcars$mpg<=22.5,-1,1)
class=data.frame(mtcars,High)
class=class[,-1]
set.seed(508)
sort<-sample(1:nrow(class),nrow(class)*0.8)
train_svm<-class[sort,]
test_svm<-class[-sort,]
svmfit=svm(train_svm$High~.,data=train_svm,kernel="linear",epsilon=0.1,cost=1,scale=FALSE)
svm.probs_train <- predict(svmfit,train_svm,type = "response")
svm.pred_train <- ifelse(svm.probs_train > 0.5, "1", "-1")
(misclass <- table(svm.pred_train, truth = train_svm$High))
svm.probs_test <- predict(svmfit,newdata=test_svm,type = "response")
svm.pred_test <- ifelse(svm.probs_test > 0.5, "1", "-1")
(misclass <- table(svm.pred_test, truth = test_svm$High))

#Decision Tree
#Categorize the Data
hist(mtcars$mpg)
#from dataset, we can pick a cutoff.  For this, we will define 22.5 as it is in the middle.  Obviously not particularly rigorous in our selection, but we'll go with it for the example.
High=ifelse(mtcars$mpg<=22.5,"No","Yes")
class=data.frame(mtcars,High)
class=class[,-1]
set.seed(508)
sort<-sample(1:nrow(class),nrow(class)*0.8)
train_dt<-class[sort,]
test_dt<-class[-sort,]
fit<-rpart(High~.,train_dt)
tree.pred_train=predict(fit,train_dt,type="class")
with(train_dt,table(tree.pred_train,High))
tree.pred_test=predict(fit,test_dt,type="class")
with(test_dt,table(tree.pred_test,High))
accuracy_tune<-function(fit){
  predict_unseen<-predict(fit,test_dt,type='class')
  table_mat<-table(test_dt$High,predict_unseen)
  accuracy_Test<-sum(diag(table_mat))/sum(table_mat)
  accuracy_Test
}
#Allow us to change the parameters (minsplit,minbucket,maxdepth,cp) and see how the model changes
control<-rpart.control(minsplit=4,minbucket=2,maxdepth=3,cp=0)
tune_fit<-rpart(High~.,data=test_dt,method='class',control=control)
accuracy_tune(tune_fit)
#Identify the tuned parameters
tune_fit$control
#Pruning the tree as a function of cp (complexity parameter)
pfit<-prune(tune_fit,cp=0)
plot(pfit, uniform=TRUE,main="Pruned Tree")
text(pfit, use.n=TRUE, all=TRUE, cex=.8)
tree.pred_train=predict(pfit,train_dt,type="class")
with(train_dt,table(tree.pred_train,High))
tree.pred_test=predict(pfit,test_dt,type="class")
with(test_dt,table(tree.pred_test,High))


#Bagging
cvcontrol<-trainControl(method="repeatedcv",number=10,allowParallel=TRUE)
train.bagg<-train(as.factor(High)~.,data=train_dt,method="treebag",trControl=cvcontrol,importance=TRUE,nbagg=50,minsplit=4,minbucket=2,maxdepth=3,cp=0)
train.bagg
plot(varImp(train.bagg))
tree.pred_train=predict(train.bagg,train_dt)
with(train_dt,table(tree.pred_train,High))
tree.pred_test=predict(train.bagg,test_dt)
with(test_dt,table(tree.pred_test,High))


#Boosting
#Random Forest
train.rf <- train(as.factor(High) ~ ., data=train_dt,method="rf",trControl=cvcontrol,importance=TRUE)
train.rf
tree.pred_train=predict(train.rf,train_dt)
with(train_dt,table(tree.pred_train,High))
tree.pred_test=predict(train.rf,test_dt)
with(test_dt,table(tree.pred_test,High))
#Random Forest Boosting
train.gbm <- train(as.factor(High) ~ ., data=train_dt,method="gbm",verbose=F,trControl=cvcontrol)
train.gbm
tree.pred_train=predict(train.gbm,train_dt)
with(train_dt,table(tree.pred_train,High))
tree.pred_test=predict(train.gbm,test_dt)
with(test_dt,table(tree.pred_test,High))


#Logistic Regression
#Break into training and test data
High=ifelse(mtcars$mpg<=22.5,0,1)
class=data.frame(mtcars,High)
class=class[,-1]
set.seed(508)
sort<-sample(1:nrow(class),nrow(class)*0.8)
train_logit<-class[sort,]
test_logit<-class[-sort,]
#Develop Logistic Regression Model on the Training Data
glm.fit <- glm(train_logit$High ~ ., data = train_logit, family = binomial)
glm.probs_train <- predict(glm.fit,train_logit,type = "response")
glm.pred_train <- ifelse(glm.probs_train > 0.5, "1", "0")
#Assess prediction result as table as yes / no accuracy
(misclass <- table(glm.pred_train, truth = train_logit$High))
#Apply model and repeat on training data
glm.probs_test <- predict(glm.fit,test_logit,type = "response")
glm.pred_test <- ifelse(glm.probs_test > 0.5, "1", "0")
(misclass <- table(glm.pred_test, truth = test_logit$High))