{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Model Fitting\n\nIn this notebook we experiement with different model fitting and tuning to determine the best model for this dataset.","metadata":{}},{"cell_type":"markdown","source":"# Library","metadata":{}},{"cell_type":"code","source":"library(tidyverse)\nlibrary(caret)\nlibrary(gridExtra)\nlibrary(ggplot2)","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2021-11-03T10:59:48.078012Z","iopub.execute_input":"2021-11-03T10:59:48.080286Z","iopub.status.idle":"2021-11-03T10:59:48.105023Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Import data","metadata":{}},{"cell_type":"code","source":"path <- \"../input/eda-fb-data/\"\nlist.files(path)","metadata":{"execution":{"iopub.status.busy":"2021-11-03T10:59:48.109310Z","iopub.execute_input":"2021-11-03T10:59:48.111767Z","iopub.status.idle":"2021-11-03T10:59:48.136083Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data <- read.csv(paste0(path, \"train_processed.csv\"), row.names=1)\ndata$outcome <- as.factor(ifelse(data$outcome == 0, \"Human\", \"Robot\"))\ndata %>% str()\ndata %>% head()","metadata":{"execution":{"iopub.status.busy":"2021-11-03T10:59:48.138372Z","iopub.execute_input":"2021-11-03T10:59:48.139759Z","iopub.status.idle":"2021-11-03T10:59:48.190077Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Train valid split","metadata":{}},{"cell_type":"code","source":"set.seed(123)\n\nPERCENT = 0.8\n\n# save 20% for validation, and 80% for training\n# we will use the training data for selecting best parameter (using CV)\n# we will use validation data to choose the best model out of the models with best param\n# due to the imbalance class, we use cluster split\n\n# selecting 80% of robots as train\nrobot_names <- row.names(data %>% filter(outcome==\"Robot\"))\nrobot_train <- sample(robot_names, PERCENT * length(robot_names))\nrobot_valid <- robot_names[!(robot_names %in% robot_train)]\n\n# selecting 80% of humans as train\nhuman_names <- row.names(data %>% filter(outcome==\"Human\"))\nhuman_train <- sample(human_names, PERCENT * length(human_names))\nhuman_valid <- human_names[!(human_names %in% human_train)]\n\ntrain <- data[row.names(data) %in% c(robot_train, human_train), ]\nvalid <- data[row.names(data) %in% c(robot_valid, human_valid), ]\n\n\n\ncat(\n    \"Total data: \", data %>% nrow(),\n    \"\\nTrain data:\", train %>% nrow(),\n    \"\\nValid data:\", valid %>% nrow()\n )","metadata":{"execution":{"iopub.status.busy":"2021-11-03T10:59:48.192429Z","iopub.execute_input":"2021-11-03T10:59:48.193866Z","iopub.status.idle":"2021-11-03T10:59:48.231797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# sanity check: make sure no repeated observations\nrownames(train) %in% rownames(valid) %>% sum()","metadata":{"execution":{"iopub.status.busy":"2021-11-03T10:59:48.234178Z","iopub.execute_input":"2021-11-03T10:59:48.235576Z","iopub.status.idle":"2021-11-03T10:59:48.251364Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data preprocessing","metadata":{}},{"cell_type":"code","source":"# normalizing based on train\n# note: although we technically can normalize testing data with training data\n# doing so is not practical in deployment\n# we will not be able to re-fit the model everytime when we have a real time update of data\npreProcValues <- preProcess(train[-1], method=c(\"center\", \"scale\"))","metadata":{"execution":{"iopub.status.busy":"2021-11-03T10:59:48.253816Z","iopub.execute_input":"2021-11-03T10:59:48.255220Z","iopub.status.idle":"2021-11-03T10:59:48.274187Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# indeed, training data has mean 0 while validation data has mean close to 0\npredict(preProcValues, train[-1]) %>% summary()\npredict(preProcValues, valid[-1]) %>% summary()","metadata":{"execution":{"iopub.status.busy":"2021-11-03T10:59:48.276546Z","iopub.execute_input":"2021-11-03T10:59:48.277896Z","iopub.status.idle":"2021-11-03T10:59:48.316147Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model fitting with hyper-parameter tunning\n\nWe have learnt quite a few classification models in this modules. Every model can be tuned to improve their performance.","metadata":{}},{"cell_type":"code","source":"# Cross Validation fitting\nfitControl <- trainControl(\n    # repeated Cross Validation\n    method = \"repeatedcv\",\n    number = 10, # k-fold CV\n    repeats = 10, # repeat n times\n    # estimate class probabilities\n    classProbs = TRUE,\n    # two class summary function\n    summaryFunction = twoClassSummary,\n    # sample\n    #sampling = \"down\"\n)","metadata":{"execution":{"iopub.status.busy":"2021-11-03T10:59:48.319701Z","iopub.execute_input":"2021-11-03T10:59:48.321408Z","iopub.status.idle":"2021-11-03T10:59:48.335385Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# twoClassSummary function example\n# we want to retrive the Area under ROC as the metric for comparision\ntwoClassSummary(\n    data.frame(\n        obs=train$outcome,\n        pred=factor(rep(\"Human\", train %>% nrow()), levels=c(\"Human\", \"Robot\")),\n        Human=rep(1.0, train %>% nrow()),\n        Robot=rep(0.0, train %>% nrow())),\n    lev=c(\"Human\", \"Robot\"))","metadata":{"execution":{"iopub.status.busy":"2021-11-03T10:59:48.337805Z","iopub.execute_input":"2021-11-03T10:59:48.339520Z","iopub.status.idle":"2021-11-03T10:59:48.358688Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# custom function to fit in the model with a given name\ntrain_model <- function(model_name, tune=NULL, training=TRUE) {\n    if(training){\n        usefitControl <- fitControl\n        useData <- train\n    } else {\n        usefitControl <- trainControl(method=\"none\", classProbs=TRUE)\n        useData <- data\n    }\n    train(outcome ~ ., data=useData,\n          # Preprocessing techniques\n          preProcess = c(\"center\", \"scale\"),\n          # Training control\n          trControl = usefitControl,\n          # metric: Area under ROC curve\n          metric = \"ROC\",\n          # specify model name here\n          method = model_name,\n          # model tunning\n          tuneGrid = tune\n         )\n}","metadata":{"execution":{"iopub.status.busy":"2021-11-03T10:59:48.361130Z","iopub.execute_input":"2021-11-03T10:59:48.362519Z","iopub.status.idle":"2021-11-03T10:59:48.375993Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# models learnt in this module\n# R doesn't have one major package that establish the standard for model implementation\n# we use Caret, which is a wrapper for different libraries\n\n####################################\n# Models that doesn't need tuning #\n###################################\n# Linear Discriminant Analysis\nmodel_lda <- train_model(\"lda\")\n\n# Quadratic Discriminant Analysis\nmodel_qda <- train_model(\"qda\")\n\n# Logistic Regerssion (under Generalised Linear Models)\n# We do not spend time fitting interaction term and polynomial terms\n# since we also includes non-linear model such as RF\nmodel_glm <- train_model(\"glm\")\n\n# Naive Bayes\n# no tuning but we add laplace as 1\ngrid_nb <- data.frame(laplace=1, usekernel=FALSE, adjust=0)\nmodel_nb <- train_model(\"naive_bayes\", tune=grid_nb)\n\n##########################\n# Model require tunning #\n#########################\n# K-Nearest Neighbourhood\n# Tune: K (number of nearest data points)\ngrid_knn <- data.frame(k=c(3, 5, 10, 15, 20, 25, 30, 40, 50))\nmodel_knn <- train_model(\"knn\", tune=grid_knn)\n\n# Support Vector (RBF Kernel or Gaussian kernel)\n# tune: C (Regularization parameter)\n# we use default sigma by the author\ngrid_svmRadial <- expand.grid(C = c(0.5, 1.0, 1.5, 3.0, 5.0), sigma = c(1.0, 5.0, 10.0))\nmodel_svmRadial <- train_model(\"svmRadial\", tune=grid_svmRadial)\n\n# Random Forest (Bagging)\n# tune: mtry (Number of parameters to try)\ngrid_rf <- data.frame(mtry=c(\n    # sqrt of the columns\n    round(sqrt(train %>% ncol() - 1)),\n    # half of the columns\n    round(0.5 * (train %>% ncol() - 1)),\n    # all the columns (bagging)\n    train %>% ncol()\n))\nmodel_rf <- train_model(\"rf\")","metadata":{"execution":{"iopub.status.busy":"2021-11-03T10:59:48.378853Z","iopub.execute_input":"2021-11-03T10:59:48.380356Z","iopub.status.idle":"2021-11-03T11:15:33.921410Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# comparing the repeated CV results across models\nresults <- resamples(\n    list(LDA = model_lda,\n         KNN = model_knn,\n         QDA = model_qda,\n         GLM = model_glm,\n         NB = model_nb,\n         svmRadial = model_svmRadial,\n         RF = model_rf\n        )\n)","metadata":{"execution":{"iopub.status.busy":"2021-11-03T11:15:33.923884Z","iopub.execute_input":"2021-11-03T11:15:33.925281Z","iopub.status.idle":"2021-11-03T11:15:34.660785Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plotting tunning results\ngrid.arrange(\n    plot(model_rf, main=\"Random Forest\"),\n    plot(model_knn, main=\"KNN\"),\n    plot(model_svmRadial, main=\"SVM Radial\")\n)","metadata":{"execution":{"iopub.status.busy":"2021-11-03T11:15:34.663185Z","iopub.execute_input":"2021-11-03T11:15:34.664563Z","iopub.status.idle":"2021-11-03T11:15:34.953320Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# We focus on ROC\nsummary(results)","metadata":{"execution":{"iopub.status.busy":"2021-11-03T11:15:34.955890Z","iopub.execute_input":"2021-11-03T11:15:34.957396Z","iopub.status.idle":"2021-11-03T11:15:34.979465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Focus on ROC\nbwplot(results)","metadata":{"execution":{"iopub.status.busy":"2021-11-03T11:15:34.981954Z","iopub.execute_input":"2021-11-03T11:15:34.983419Z","iopub.status.idle":"2021-11-03T11:15:35.181000Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# checking parameters\nmodel_rf$finalModel","metadata":{"execution":{"iopub.status.busy":"2021-11-03T11:15:35.184272Z","iopub.execute_input":"2021-11-03T11:15:35.186073Z","iopub.status.idle":"2021-11-03T11:15:35.207971Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model_knn$finalModel","metadata":{"execution":{"iopub.status.busy":"2021-11-03T11:15:35.211991Z","iopub.execute_input":"2021-11-03T11:15:35.213825Z","iopub.status.idle":"2021-11-03T11:15:35.232273Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model_svmRadial$finalModel","metadata":{"execution":{"iopub.status.busy":"2021-11-03T11:15:35.234896Z","iopub.execute_input":"2021-11-03T11:15:35.236376Z","iopub.status.idle":"2021-11-03T11:15:35.252251Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Estimate results with validation split\n\nWe will choose the best model based on parameters found using training set.\\\nWe will use validation set to choose the best model.\\\nWe will only compare the top3 models tunned.","metadata":{}},{"cell_type":"code","source":"# final model\n\n####################################\n# Models that doesn't need tuning #\n###################################\n# Linear Discriminant Analysis\nfinal_lda <- train_model(\"lda\", training=FALSE)\n\n# Quadratic Discriminant Analysis\nfinal_qda <- train_model(\"qda\", training=FALSE)\n\n# Logistic Regerssion (under Generalised Linear Models)\n# We do not spend time fitting interaction term and polynomial terms\n# since we also includes non-linear model such as RF\nfinal_glm <- train_model(\"glm\", training=FALSE)\n\n# Naive Bayes\n# no tuning but we add laplace as 1\ngrid_nb <- data.frame(laplace=1, usekernel=FALSE, adjust=0)\nfinal_nb <- train_model(\"naive_bayes\", tune=grid_nb, training=FALSE)\n\n##########################\n# Model require tunning #\n#########################\n# K-Nearest Neighbourhood\n# Tune: K (number of nearest data points)\ngrid_knn <- data.frame(k=c(40))\nfinal_knn <- train_model(\"knn\", tune=grid_knn, training=FALSE)\n\n# Support Vector (RBF Kernel or Gaussian kernel)\n# tune: C (Regularization parameter)\n# we use default sigma by the author\ngrid_svmRadial <- data.frame(C = c(0.5), sigma = c(1.0))\nfinal_svmRadial <- train_model(\"svmRadial\", tune=grid_svmRadial, training=FALSE)\n\n# Random Forest (Bagging)\n# tune: mtry (Number of parameters to try)\ngrid_rf <- data.frame(mtry=c(2))\nfinal_rf <- train_model(\"rf\", training=FALSE)","metadata":{"execution":{"iopub.status.busy":"2021-11-03T11:15:35.254843Z","iopub.execute_input":"2021-11-03T11:15:35.256255Z","iopub.status.idle":"2021-11-03T11:15:40.384523Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pred_ROC <- function(model, X, prob=FALSE){\n    pred_prob = predict(model, X, type=\"prob\")\n    pred <- data.frame(\n        obs=X$outcome,\n        pred=predict(model, X),\n        Human=pred_prob[\"Human\"],\n        Robot=pred_prob[\"Robot\"]\n    )\n    result <- twoClassSummary(\n        pred,\n        lev=c(\"Human\", \"Robot\")\n    )[[\"ROC\"]]\n    if(prob){\n        pred\n    } else {\n        result\n    }\n}","metadata":{"execution":{"iopub.status.busy":"2021-11-03T11:15:40.386964Z","iopub.execute_input":"2021-11-03T11:15:40.388420Z","iopub.status.idle":"2021-11-03T11:15:40.400700Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data.frame(\n    Methods=c(\"RF\", \"KNN\", \"LDA\", \"QDA\", \"GLM\", \"NB\", \"SVM\"),\n    ROC=c(\n        pred_ROC(final_rf, valid),\n        pred_ROC(final_knn, valid),\n        pred_ROC(final_lda, valid),\n        pred_ROC(final_qda, valid),\n        pred_ROC(final_glm, valid),\n        pred_ROC(final_nb, valid),\n        pred_ROC(final_svmRadial, valid)\n    )\n) %>% arrange(desc(ROC))","metadata":{"execution":{"iopub.status.busy":"2021-11-03T11:15:40.403033Z","iopub.execute_input":"2021-11-03T11:15:40.404474Z","iopub.status.idle":"2021-11-03T11:15:40.740150Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Our best model is RandomForest with 2 random features per split.","metadata":{}},{"cell_type":"markdown","source":"# Error analysis\n\nLet's take a look at our rf prediction","metadata":{}},{"cell_type":"code","source":"pred_rf <- pred_ROC(final_rf, valid, prob=TRUE)\n# col is for true, row is for predict\ntable(pred_rf$pred, pred_rf$obs)","metadata":{"execution":{"iopub.status.busy":"2021-11-03T11:15:40.748421Z","iopub.execute_input":"2021-11-03T11:15:40.754241Z","iopub.status.idle":"2021-11-03T11:15:40.828744Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rf_feature <- final_rf$finalModel$importance %>% t()","metadata":{"execution":{"iopub.status.busy":"2021-11-03T11:15:40.831125Z","iopub.execute_input":"2021-11-03T11:15:40.832520Z","iopub.status.idle":"2021-11-03T11:15:40.844699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rf_feature %>% barplot(las=2, horiz=FALSE)","metadata":{"execution":{"iopub.status.busy":"2021-11-03T11:15:40.847242Z","iopub.execute_input":"2021-11-03T11:15:40.848746Z","iopub.status.idle":"2021-11-03T11:15:40.917140Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"summary(final_glm$finalModel)","metadata":{"execution":{"iopub.status.busy":"2021-11-03T11:15:40.919678Z","iopub.execute_input":"2021-11-03T11:15:40.921101Z","iopub.status.idle":"2021-11-03T11:15:40.939944Z"},"trusted":true},"execution_count":null,"outputs":[]}]}