rm(list=ls())
####################################################
####           Load required packages           ####
####################################################
installIfAbsentAndLoad <- function(neededVector) {
  for(thispackage in neededVector) {
    if( ! require(thispackage, character.only = T) )
    { install.packages(thispackage)}
    require(thispackage, character.only = T)
  }
}
needed <- c('e1071','MLeval','caret','MLmetrics','ROCR','PRROC','InformationValue',
            'rpart','tidyverse','data.table','GGally','corrplot','verification',
            'ROCR','maptree','glmnet','gridExtra','randomForest','mgcv','nnet','pROC',
            'gbm','DT',"NeuralNetTools",'rpart.plot','Matrix','tree')  
installIfAbsentAndLoad(needed)
require(methods)

########################################################
#### Importing and preprocessing the data           ####
########################################################
set.seed(21)
df <- read.table("train.csv", header=T, sep=',')
df<-df[, -1]
# Create test and train data sets 
n <- nrow(df)
train.indices <- sample(n, .8 * n)
df.train <- df[train.indices, ] # training data 80% split
df.train$target <- as.factor(df.train$target) #train y

df.test <- df[-train.indices, ] # test set 20%
df.test$target <- as.factor(df.test$target) # test y

df.train.scale <- df.train # copy of original training set
df.train.scale[,2:201] <- scale(df.train.scale[,2:201])  # scale the x variables 
df.test.scale <- df.test # copy of original test set
df.test.scale[,2:201] <- scale(df.test.scale[,2:201])


performance_colnames <- c("model", "accuracy", "f1", "precision", "recall",
                          "rocAUC", "prAUC")

########################################################
####       Implementing different models            ####
########################################################

# logistic regresion
## Standard using non scaled x features
set.seed(21)
model.name <-  "Logistic"
logistic.fit <- glm(target~., data= df.train, family=binomial)
logistic.pred.probs <- predict(logistic.fit, df.test, type="response")  

#optCutOff <- optimalCutoff(df.test$target, logistic.pred.probs)[1] #generate a cutoff threshold that gives minimum misclassification error
logistic.model.pred <-  ifelse(logistic.pred.probs > 0.5, 1, 0)
# Evaluate performance of the model on the test set
accuracy <- Accuracy(logistic.model.pred, df.test$target)
f1 <- F1_Score(logistic.model.pred, df.test$target)
prec <- Precision(logistic.model.pred, df.test$target)
rec <- Recall(logistic.model.pred, df.test$target)
rocauc <- AUC(logistic.model.pred, df.test$target)
prauc <- PRAUC(logistic.model.pred, df.test$target)
performance <- data.frame( model.name, accuracy, f1, prec, rec, 
                           rocauc, prauc, stringsAsFactors = F)
colnames(performance) <- performance_colnames

## Ridge using the scaled data
set.seed(21)
model.name <- "Ridge"
x.data.train <- model.matrix(target~., df.train.scale)[,-1] # convert Xs into a usable form for glmnet
x.data.test <- model.matrix(target~., df.test.scale)[,-1]  # convert Xs into a usable form for glmnet

  ## run cv to determine the best lambda  
lambdagrid <- 10^seq(10,-3, length=100)
cv.out.class <-  cv.glmnet(x.data.train, df.train.scale$target, alpha =0 , 
                           lambda = lambdagrid,  family = 'binomial')
(bestlam <-  cv.out.class$lambda.min) # best lambda is 0.001
bestlam=0.001
#run again with best lambda =0.01
mod.ridge <-  glmnet(x.data.train, df.train.scale$target, alpha = 0, 
                     lambda=bestlam, family = 'binomial')
ridge.model.pred<-  rep(0, length(df.test.scale$target)) 
prob  <-  predict(mod.ridge,s = bestlam, newx = x.data.test, type = 'response')

#optCutOff <- optimalCutoff(df.test$target, prob)[1]
ridge.model.pred[prob>0.5] <- 1

#calculate performance measures 
accuracy <- Accuracy(ridge.model.pred, df.test$target)
f1 <- F1_Score(ridge.model.pred, df.test$target)
prec <- Precision(ridge.model.pred, df.test$target)
rec <- Recall(ridge.model.pred, df.test$target)
rocauc <- AUC(ridge.model.pred, df.test$target)
prauc <- PRAUC(ridge.model.pred, df.test$target)
performance <- rbind(performance, c( model.name, accuracy, f1, prec, rec, rocauc, prauc))
colnames(performance) <- performance_colnames

## Lasso using scaled data
set.seed(21)
model.name <- "Lasso"

# run cv to determine the best lambda  
cv.out.class <-  cv.glmnet(x.data.train, df.train.scale$target, alpha =1 , 
                           lambda = lambdagrid,  family = 'binomial') #run
(bestlam <-  cv.out.class$lambda.min) # best lambda is 0.001
bestlam <- 0.001
#run again with best lambda
mod.lasso <-  glmnet(x.data.train, df.train.scale$target, alpha = 1, 
                     lambda=bestlam, family = 'binomial')
lasso.model.pred<-  rep(0, length(df.test.scale$target))
prob  <-  predict(mod.lasso,s = bestlam, newx = x.data.test, type = 'response')
#optCutOff <- optimalCutoff(df.test$target, prob)[1]
lasso.model.pred[prob>.5] <- 1

#calculate performance measures 
accuracy <- Accuracy(lasso.model.pred, df.test$target)
f1 <- F1_Score(lasso.model.pred, df.test$target)
prec <- Precision(lasso.model.pred, df.test$target)
rec <- Recall(lasso.model.pred, df.test$target)
rocauc <- AUC(lasso.model.pred, df.test$target)
prauc <- PRAUC(lasso.model.pred, df.test$target)
performance <- rbind(performance, c( model.name, accuracy, f1, prec, rec, rocauc, prauc))
colnames(performance) <- performance_colnames

# Tree Methods
## Tree classification
set.seed(21)
model.name <- "Tree Classification"
tree.mod<- tree(target~., df.train)
tree.mod.pred<-predict(tree.mod, df.test, type = 'class')

#calculate performance measures 
accuracy <- Accuracy(tree.mod.pred, df.test$target)
f1 <- F1_Score(tree.mod.pred, df.test$target)
prec <- Precision(tree.mod.pred, df.test$target)
rec <- Recall(tree.mod.pred, df.test$target)
rocauc <- AUC(tree.mod.pred, df.test$target)
prauc <- PRAUC(tree.mod.pred, df.test$target)
performance <- rbind(performance, c( model.name, accuracy, f1, prec, rec, rocauc, prauc))
colnames(performance) <- performance_colnames


## Random Forest
set.seed(21)
model.name <- "Random Forrest"
library(rfUtilities)
# Train on subset of data
c1 <- randomForest(target ~ ., df.train[1:40000,], ntree=50, norm.votes=FALSE, importance=TRUE)
c2 <- randomForest(target ~ ., df.train[40001:80000,], ntree=50, norm.votes=FALSE, importance=TRUE)
c3 <- randomForest(target ~ ., df.train[80001:120000,], ntree=50, norm.votes=FALSE, importance=TRUE)
c4 <- randomForest(target ~ ., df.train[120001:160000,], ntree=50, norm.votes=FALSE, importance=TRUE)
# Combine all the random forests
class.combine <- rf.combine(c1, c2, c3, c4) 
# Predict on the combined random forest
pred.rf=predict(class.combine, newdata=df.test, type="response")
# Get the scoring metrics
accuracy <- Accuracy(pred.rf, df.test$target)
f1 <- F1_Score(pred.rf, df.test$target)
prec <- Precision(pred.rf, df.test$target)
rec <- Recall(pred.rf, df.test$target)
rocauc <- AUC(pred.rf, df.test$target)
prauc <- PRAUC(pred.rf, df.test$target)
performance <- rbind(performance, c( model.name, accuracy, f1, prec, rec, rocauc, prauc))
colnames(performance) <- performance_colnames

#write.csv(performance, "model_performance.csv", row.names = FALSE)
#performance <- read.csv("model_performance.csv")




