{"cells":[{"metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","trusted":true},"cell_type":"code","source":"## Importing packages\n\n# This R environment comes with all of CRAN and many other helpful packages preinstalled.\n# You can see which packages are installed by checking out the kaggle/rstats docker image: \n# https://github.com/kaggle/docker-rstats\n\nlibrary(tidyverse) # metapackage with lots of helpful functions\n\n## Running code\n\n# In a notebook, you can run a single code cell by clicking in the cell and then hitting \n# the blue arrow to the left, or by clicking in the cell and pressing Shift+Enter. In a script, \n# you can run code by highlighting the code you want to run and then clicking the blue arrow\n# at the bottom of this window.\n\n## Reading in files\n\n# You can access files from datasets you've added to this kernel in the \"../input/\" directory.\n# You can see the files added to this kernel by running the code below. \n\nlist.files(path = \"../input/\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Check what data files we have"},{"metadata":{"trusted":true},"cell_type":"code","source":"list.files(path=\"../input/preanalysis/\")\ntest_data <- read.csv(\"../input/preanalysis/analysis/analysis.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Check missing values (NA)\n# sapply(test_data, function(x) mean(is.na(x)))\ndim(test_data)\n# Remove records with NA\ntest_data <- test_data[complete.cases(test_data), ]\ndim(test_data)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Fundamental Analysis"},{"metadata":{"trusted":true},"cell_type":"code","source":"ynames <- c(\"DM_M1\",\"DM_M7\",\"DM_M28\",\"DM_M42\",\"BodyPart\",\"Surface\",\"PlayerGame\",\"BodyPart_Ankle\",\"BodyPart_Foot\",\"BodyPart_Heel\",\"BodyPart_Knee\",\"BodyPart_Toes\",\"Surface_Natural\",\"Surface_Synthetic\")\nkeys <- c(\"PlayerKey\",\"GameID\",\"PlayKey\")\ncategorical_names <- c(\"RosterPosition\",\"PlayerDay\",\"StadiumType\",\"FieldType\",\"Temperature\",\"Weather\",\"PlayType\",\"PlayerGamePlay\",\"Position\",\"PositionGroup\")\ndummies <- colnames(test_data)[!names(test_data) %in% c(ynames, keys, categorical_names)] \n\ntest_data$tem <- floor((test_data$Temperature - 0)/10)\ntest_data$pdg <- floor((test_data$PlayerDay - 0)/50)\n\ninjury_data <- test_data[test_data$DM_M1 == 1,]\n\ncompare_data <- function(col_name){\n\ta <- table(test_data[,col_name])/nrow(test_data)\n\tb <- table(injury_data[,col_name])/nrow(injury_data)\n\tadd_lists = setdiff(names(a),names(b))\n\tif (length(add_lists) > 0) {\n\t\tfor (l in add_lists){\n\t\t\tb[l] <- 0\n\t\t}\n\t} \n\tb <- b[names(a)]\n\tagg <- rbind(a,b)\n\tbarplot(agg, main=paste0(\"Distribution of \", col_name),\n\t  xlab=col_name, ylab=\"Percentage\", col=c(\"darkblue\",\"lightblue\"),\n\t  args.legend = list(x = \"topright\", bty = \"n\", inset=c(0.15, 0)),\n\t  legend = c(\"all\",\"injury\"), beside=TRUE)\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's check some different patterns between injuries and all data"},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_data(\"RosterPosition\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_data(\"StadiumType\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_data(\"FieldType\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_data(\"Position\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_data('PositionGroup')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_data('PlayType')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_data('tem') ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_data('pdg')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_data('Weather')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"compare_data('Cell_19')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"library(lattice)\nhistogram(~ a_99th | DM_M1, data=test_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"histogram(~ s_99th | DM_M1, data=test_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"histogram(~ dir_chg_99th | DM_M1, data=test_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"histogram(~ o_chg_99th | DM_M1, data=test_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"histogram(~ a_avg | DM_M1, data=test_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"histogram(~ s_avg | DM_M1, data=test_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"histogram(~ dir_chg_avg | DM_M1, data=test_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"histogram(~ o_chg_avg | DM_M1, data=test_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Let's check some injury data\nsummary(injury_data[, colnames(injury_data) %in% c(\"DM_M1\",\"DM_M7\",\"DM_M28\",\"DM_M42\")])\nys <- apply(injury_data[, colnames(injury_data) %in% c(\"DM_M1\",\"DM_M7\",\"DM_M28\",\"DM_M42\")],2,sum)\nbarplot(ys)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"histogram(~ DM_M7 | Surface, data=injury_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"histogram(~ DM_M28 | Surface, data=injury_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"histogram(~ DM_M42 | Surface, data=injury_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"histogram(~ DM_M7 | BodyPart, data=injury_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"histogram(~ DM_M28 | BodyPart, data=injury_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"histogram(~ DM_M42 | BodyPart, data=injury_data)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"library(rpart)\nlibrary(FNN)\nlibrary(neuralnet)\nlibrary(gbm)\nlibrary(glmnet)\nlibrary(randomForest)\nlibrary(caret)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"raw","source":"Predictive Analysis"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Remove some data without injury to speed up\nset.seed(123)\nnoinjury <- test_data[ sample( which( test_data$DM_M1 == 0 ) , 10000 ) , ]\ninjury <- test_data[test_data$DM_M1 == 1,]\nnew_data <- rbind(noinjury, injury)\n\n# Remove columns with constant values\nnew_data <- new_data[,!apply(new_data, MARGIN = 2, function(x) max(x, na.rm = TRUE) == min(x, na.rm = TRUE))]\n\ndummies <- dummies[dummies %in% colnames(new_data)]\n\n# Remove one of highly correlated pair\nthreshold <- 0.99\nremove_list <- c()\ncounter <- 2\nfor (a_dummy in dummies[1:(length(dummies)-1)]){\n#\tprint(a_dummy)\n\tfor (b_dummy in dummies[counter:length(dummies)]){\n#\t\tprint(b_dummy)\n\t\tif (abs(cor(new_data[,a_dummy],new_data[,b_dummy])) > threshold){\n\t\t\tremove_list <- c(remove_list, a_dummy)\n\t\t\tprint(paste0(a_dummy, \" will be removed. Correlation between \", a_dummy, \" and \", b_dummy, \" is \", cor(new_data[,a_dummy],new_data[,b_dummy])))\n\t\t\tbreak\n\t\t}\n\t}\n\tcounter <- counter + 1\n}\n\nnew_data <- new_data[, !colnames(new_data) %in% remove_list]\n\ndummies <- dummies[dummies %in% colnames(new_data)]\n\n# Spliting between training data and validation data\nset.seed(6)\nidx <- sample(seq(1, 2), size = nrow(new_data), replace = TRUE, prob = c(.8, .2))\n\ntraining <- new_data[idx==1,]\nvalidation <- new_data[idx==2,]\nprint(dim(training))\nprint(dim(validation))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"yname <- 'DM_M1'\n#formula\nf <-as.formula(paste(yname,\"~\",paste(dummies,collapse=\"+\")))\nprint(f)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"modeloutput <- data.frame(y=character(),\n                 Precision_train=double(), #precision based on training data\n                 Recall_train=double(),\t   #recall based on training data\n                 FMeasure_train=double(),  #F measure based on training data\n                 Precision_valid=double(), #precision based on validation data\n                 Recall_valid=double(),    #recall based on validation data\n                 FMeasure_valid=double(),  #F measure based on validation data\n\t\t\t\t tp=double(),              #true positive\n\t\t\t\t tn=double(),\t\t\t   #true negative\n\t\t\t\t fp=double(),\t\t\t   #false positive\n\t\t\t\t fn=double(),\t\t\t   #false negative\n                 Models=character(),\t   #Model type\n                 stringsAsFactors=FALSE)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Logistic Model\nstart_time <- Sys.time()\n\nglmr <- glm(f, data=training, family=binomial)\nglmr <- glm(f, data=training, family=binomial)\nglmpredict <- predict(glmr, training)\npredict <- ifelse(glmpredict>0.5, 1, 0)\npredictyes <- sum(predict==1)\npredictno<-sum(predict==0)\nactual <- training[,yname]\nactualwhenpredyes <- actual[predict==1]\nactualwhenpredno <- actual[predict==0]\ntp <- sum(actualwhenpredyes==1)\ntn <- sum(actualwhenpredno==0)\nfp <- length(actualwhenpredyes) - tp\nfn <- length(actualwhenpredno) - tn\nprecision<-tp/predictyes\nrecall<-tp/(tp+fn)\nF = 2*precision*recall/(precision+recall)\npaste(\"Training: The precision is\",precision)\npaste(\"Training: The recall is\",recall)\npaste(\"Training: The F is\",F)\n\nglmpredict <- predict(glmr, validation)\npredict <- ifelse(glmpredict>0.5, 1, 0)\npredictyes <- sum(predict==1)\npredictno<-sum(predict==0)\nactual <- validation[,yname]\nactualwhenpredyes <- actual[predict==1]\nactualwhenpredno <- actual[predict==0]\ntp <- sum(actualwhenpredyes==1)\ntn <- sum(actualwhenpredno==0)\nfp <- length(actualwhenpredyes) - tp\nfn <- length(actualwhenpredno) - tn\nprecision_v<-tp/predictyes\nrecall_v<-tp/(tp+fn)\nF_v = 2*precision_v*recall_v/(precision_v+recall_v)\npaste(\"Validation: The precision is\",precision_v)\npaste(\"Validation: The recall is\",recall_v)\npaste(\"Validation: The F is\",F_v)\n\nprint(paste0(\"tp = \",tp))\nprint(paste0(\"tn = \",tn))\nprint(paste0(\"fp = \",fp))\nprint(paste0(\"fn = \",fn))\n\nend_time <- Sys.time()\nprint(paste0(\"total run time is \", end_time - start_time, \"mins\"))\n\nmodeloutput[nrow(modeloutput)+1,] <- c(\"DM_M1\", precision, recall, F, precision_v, recall_v, F_v, tp, tn, fp, fn, 'Logistic')\nwrite.csv(glmr$coefficients,paste0(yname,'glmr_coef.csv'))\n\n\n# #remove variables that caused issues in glm\n# remove_list <- dummies[!dummies %in% rownames(summary(glmr)$coefficients)]\n# dummies <- dummies[!dummies %in% remove_list]\n# yname <- 'DM_M1'\n# #formula\n# f <-as.formula(paste(yname,\"~\",paste(dummies,collapse=\"+\")))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#cart\nstart_time <- Sys.time()\n\ncart = rpart(f, data = training, cp = 10^(-3),minsplit = 10,, method = \"class\")\ncartpredict <- predict(cart, training)\npredict <- ifelse(cartpredict[,2]>0.5, 1, 0)\npredictyes <- sum(predict==1)\npredictno<-sum(predict==0)\nactual <- training[,yname]\nactualwhenpredyes <- actual[predict==1]\nactualwhenpredno <- actual[predict==0]\ntp <- sum(actualwhenpredyes==1)\ntn <- sum(actualwhenpredno==0)\nfp <- length(actualwhenpredyes) - tp\nfn <- length(actualwhenpredno) - tn\nprecision<-tp/predictyes\nrecall<-tp/(tp+fn)\nF = 2*precision*recall/(precision+recall)\npaste(\"Training: The precision is\",precision)\npaste(\"Training: The recall is\",recall)\npaste(\"Training: The F is\",F)\n\ncartpredict <- predict(cart, validation)\npredict <- ifelse(cartpredict[,2]>0.5, 1, 0)\npredictyes <- sum(predict==1)\npredictno<-sum(predict==0)\nactual <- validation[,yname]\nactualwhenpredyes <- actual[predict==1]\nactualwhenpredno <- actual[predict==0]\ntp <- sum(actualwhenpredyes==1)\ntn <- sum(actualwhenpredno==0)\nfp <- length(actualwhenpredyes) - tp\nfn <- length(actualwhenpredno) - tn\nprecision_v<-tp/predictyes\nrecall_v<-tp/(tp+fn)\nF_v = 2*precision_v*recall_v/(precision_v+recall_v)\npaste(\"Validation: The precision is\",precision_v)\npaste(\"Validation: The recall is\",recall_v)\npaste(\"Validation: The F is\",F_v)\n\nprint(paste0(\"tp = \",tp))\nprint(paste0(\"tn = \",tn))\nprint(paste0(\"fp = \",fp))\nprint(paste0(\"fn = \",fn))\n\nend_time <- Sys.time()\nprint(paste0(\"total run time is \", end_time - start_time, \"mins\"))\n\nmodeloutput[nrow(modeloutput)+1,] <- c(\"DM_M1\", precision, recall, F, precision_v, recall_v, F_v, tp, tn, fp, fn, 'CART')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#knn\nstart_time <- Sys.time()\n\nknnreg <- knn.reg(train = training[,names(training) %in% dummies] , test=training[,names(training) %in% dummies], y=training[,yname], k=10, algorithm = \"kd_tree\")\npredict <- ifelse(knnreg$pred>0.5, 1, 0)\npredictyes <- sum(predict==1)\npredictno<-sum(predict==0)\nactual <- training[,yname]\nactualwhenpredyes <- actual[predict==1]\nactualwhenpredno <- actual[predict==0]\ntp <- sum(actualwhenpredyes==1)\ntn <- sum(actualwhenpredno==0)\nfp <- length(actualwhenpredyes) - tp\nfn <- length(actualwhenpredno) - tn\nprecision<-tp/predictyes\nrecall<-tp/(tp+fn)\nF = 2*precision*recall/(precision+recall)\npaste(\"Training: The precision is\",precision)\npaste(\"Training: The recall is\",recall)\npaste(\"Training: The F is\",F)\n\nknnreg <- knn.reg(train = training[,names(training) %in% dummies] , test=validation[,names(validation) %in% dummies], y=training[,yname], k=10, algorithm = \"kd_tree\")\npredict <- ifelse(knnreg$pred>0.5, 1, 0)\npredictyes <- sum(predict==1)\npredictno<-sum(predict==0)\nactual <- validation[,yname]\nactualwhenpredyes <- actual[predict==1]\nactualwhenpredno <- actual[predict==0]\ntp <- sum(actualwhenpredyes==1)\ntn <- sum(actualwhenpredno==0)\nfp <- length(actualwhenpredyes) - tp\nfn <- length(actualwhenpredno) - tn\nprecision_v<-tp/predictyes\nrecall_v<-tp/(tp+fn)\nF_v = 2*precision_v*recall_v/(precision_v+recall_v)\npaste(\"Validation: The precision is\",precision_v)\npaste(\"Validation: The recall is\",recall_v)\npaste(\"Validation: The F is\",F_v)\n\nprint(paste0(\"tp = \",tp))\nprint(paste0(\"tn = \",tn))\nprint(paste0(\"fp = \",fp))\nprint(paste0(\"fn = \",fn))\n\nend_time <- Sys.time()\nprint(paste0(\"total run time is \", end_time - start_time, \"mins\"))\n\nmodeloutput[nrow(modeloutput)+1,] <- c(\"DM_M1\", precision, recall, F, precision_v, recall_v, F_v, tp, tn, fp, fn, 'KNN')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#linear regression\nstart_time <- Sys.time()\n\nlr<-lm(f, data=training)\npredict <- ifelse(lr$fitted.values>0.5, 1, 0)\npredictyes <- sum(predict==1)\npredictno<-sum(predict==0)\nactual <- training[,yname]\nactualwhenpredyes <- actual[predict==1]\nactualwhenpredno <- actual[predict==0]\ntp <- sum(actualwhenpredyes==1)\ntn <- sum(actualwhenpredno==0)\nfp <- length(actualwhenpredyes) - tp\nfn <- length(actualwhenpredno) - tn\nprecision<-tp/predictyes\nrecall<-tp/(tp+fn)\nF = 2*precision*recall/(precision+recall)\npaste(\"Training: The precision is\",precision)\npaste(\"Training: The recall is\",recall)\npaste(\"Training: The F is\",F)\n\npredict <- ifelse(predict(lr,validation)>0.5, 1, 0)\npredictyes <- sum(predict==1)\npredictno<-sum(predict==0)\nactual <- validation[,yname]\nactualwhenpredyes <- actual[predict==1]\nactualwhenpredno <- actual[predict==0]\ntp <- sum(actualwhenpredyes==1)\ntn <- sum(actualwhenpredno==0)\nfp <- length(actualwhenpredyes) - tp\nfn <- length(actualwhenpredno) - tn\nprecision_v<-tp/predictyes\nrecall_v<-tp/(tp+fn)\nF_v = 2*precision_v*recall_v/(precision_v+recall_v)\npaste(\"Validation: The precision is\",precision_v)\npaste(\"Validation: The recall is\",recall_v)\npaste(\"Validation: The F is\",F_v)\n\nprint(paste0(\"tp = \",tp))\nprint(paste0(\"tn = \",tn))\nprint(paste0(\"fp = \",fp))\nprint(paste0(\"fn = \",fn))\n\nend_time <- Sys.time()\nprint(paste0(\"total run time is \", end_time - start_time, \"mins\"))\n\nmodeloutput[nrow(modeloutput)+1,] <- c(\"DM_M1\", precision, recall, F, precision_v, recall_v, F_v, tp, tn, fp, fn, 'Linear')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#lasso\nlambda = 0.5\nstart_time <- Sys.time()\n\nlasso<-glmnet(as.matrix(training[,dummies]), as.matrix(training[,yname]), alpha = 1, lambda = lambda)\ncv.out <- cv.glmnet(as.matrix(training[,dummies]), as.matrix(training[,yname]), alpha = 1)\nlambda <- cv.out$lambda.min\nlassopredict <- predict(lasso, s = lambda, newx = as.matrix(training[dummies]))\npredict <- ifelse(lassopredict>0.5, 1, 0)\npredictyes <- sum(predict==1)\npredictno<-sum(predict==0)\nactual <- training[,yname]\nactualwhenpredyes <- actual[predict==1]\nactualwhenpredno <- actual[predict==0]\ntp <- sum(actualwhenpredyes==1)\ntn <- sum(actualwhenpredno==0)\nfp <- length(actualwhenpredyes) - tp\nfn <- length(actualwhenpredno) - tn\nprecision<-tp/predictyes\nrecall<-tp/(tp+fn)\nF = 2*precision*recall/(precision+recall)\npaste(\"Training: The precision is\",precision)\npaste(\"Training: The recall is\",recall)\npaste(\"Training: The F is\",F)\n\nlassopredict <- predict(lasso, s = lambda, newx = as.matrix(validation[dummies]))\npredict <- ifelse(lassopredict>0.5, 1, 0)\npredictyes <- sum(predict==1)\npredictno<-sum(predict==0)\nactual <- validation[,yname]\nactualwhenpredyes <- actual[predict==1]\nactualwhenpredno <- actual[predict==0]\ntp <- sum(actualwhenpredyes==1)\ntn <- sum(actualwhenpredno==0)\nfp <- length(actualwhenpredyes) - tp\nfn <- length(actualwhenpredno) - tn\nprecision_v<-tp/predictyes\nrecall_v<-tp/(tp+fn)\nF_v = 2*precision_v*recall_v/(precision_v+recall_v)\npaste(\"Validation: The precision is\",precision_v)\npaste(\"Validation: The recall is\",recall_v)\npaste(\"Validation: The F is\",F_v)\n\nprint(paste0(\"tp = \",tp))\nprint(paste0(\"tn = \",tn))\nprint(paste0(\"fp = \",fp))\nprint(paste0(\"fn = \",fn))\n\nend_time <- Sys.time()\nprint(paste0(\"total run time is \", end_time - start_time, \"mins\"))\n\nmodeloutput[nrow(modeloutput)+1,] <- c(\"DM_M1\", precision, recall, F, precision_v, recall_v, F_v, tp, tn, fp, fn, 'Lasso')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#ridge\nstart_time <- Sys.time()\n\nridge<-glmnet(as.matrix(training[,dummies]), as.matrix(training[,yname]), alpha = 0, lambda = lambda)\ncv.out <- cv.glmnet(as.matrix(training[,dummies]), as.matrix(training[,yname]), alpha = 0)\nlambda <- cv.out$lambda.min\nridgepredict <- predict(ridge, s = lambda, newx = as.matrix(training[dummies]))\npredict <- ifelse(ridgepredict>0.5, 1, 0)\npredictyes <- sum(predict==1)\npredictno<-sum(predict==0)\nactual <- training[,yname]\nactualwhenpredyes <- actual[predict==1]\nactualwhenpredno <- actual[predict==0]\ntp <- sum(actualwhenpredyes==1)\ntn <- sum(actualwhenpredno==0)\nfp <- length(actualwhenpredyes) - tp\nfn <- length(actualwhenpredno) - tn\nprecision<-tp/predictyes\nrecall<-tp/(tp+fn)\nF = 2*precision*recall/(precision+recall)\npaste(\"Training: The precision is\",precision)\npaste(\"Training: The recall is\",recall)\npaste(\"Training: The F is\",F)\n\nridgepredict <- predict(ridge, s = lambda, newx = as.matrix(validation[dummies]))\npredict <- ifelse(ridgepredict>0.5, 1, 0)\npredictyes <- sum(predict==1)\npredictno<-sum(predict==0)\nactual <- validation[,yname]\nactualwhenpredyes <- actual[predict==1]\nactualwhenpredno <- actual[predict==0]\ntp <- sum(actualwhenpredyes==1)\ntn <- sum(actualwhenpredno==0)\nfp <- length(actualwhenpredyes) - tp\nfn <- length(actualwhenpredno) - tn\nprecision_v<-tp/predictyes\nrecall_v<-tp/(tp+fn)\nF_v = 2*precision_v*recall_v/(precision_v+recall_v)\npaste(\"Validation: The precision is\",precision_v)\npaste(\"Validation: The recall is\",recall_v)\npaste(\"Validation: The F is\",F_v)\n\nprint(paste0(\"tp = \",tp))\nprint(paste0(\"tn = \",tn))\nprint(paste0(\"fp = \",fp))\nprint(paste0(\"fn = \",fn))\n\nend_time <- Sys.time()\nprint(paste0(\"total run time is \", end_time - start_time, \"mins\"))\n\nmodeloutput[nrow(modeloutput)+1,] <- c(\"DM_M1\", precision, recall, F, precision_v, recall_v, F_v, tp, tn, fp, fn, 'Ridge')\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#gbm\nstart_time <- Sys.time()\n\nn_gbm_trees = 1000\n\nset.seed(123) #reset random seed as gbm may use random subset\ngbmreg<-gbm(f, data=training, distribution = \"bernoulli\", interaction.depth=10, n.minobsinnode = 1, bag.fraction=0.7, n.trees = n_gbm_trees, shrinkage = 0.1)\ngbmpredict <- predict(gbmreg, training, n.trees = n_gbm_trees)\npredict <- ifelse(gbmpredict>0.5, 1, 0)\npredictyes <- sum(predict==1)\npredictno<-sum(predict==0)\nactual <- training[,yname]\nactualwhenpredyes <- actual[predict==1]\nactualwhenpredno <- actual[predict==0]\ntp <- sum(actualwhenpredyes==1)\ntn <- sum(actualwhenpredno==0)\nfp <- length(actualwhenpredyes) - tp\nfn <- length(actualwhenpredno) - tn\nprecision<-tp/predictyes\nrecall<-tp/(tp+fn)\nF = 2*precision*recall/(precision+recall)\npaste(\"Training: The precision is\",precision)\npaste(\"Training: The recall is\",recall)\npaste(\"Training: The F is\",F)\n\ngbmpredict <- predict(gbmreg, validation, n.trees = n_gbm_trees)\npredict <- ifelse(gbmpredict>0.5, 1, 0)\npredictyes <- sum(predict==1)\npredictno<-sum(predict==0)\nactual <- validation[,yname]\nactualwhenpredyes <- actual[predict==1]\nactualwhenpredno <- actual[predict==0]\ntp <- sum(actualwhenpredyes==1)\ntn <- sum(actualwhenpredno==0)\nfp <- length(actualwhenpredyes) - tp\nfn <- length(actualwhenpredno) - tn\nprecision_v<-tp/predictyes\nrecall_v<-tp/(tp+fn)\nF_v = 2*precision_v*recall_v/(precision_v+recall_v)\npaste(\"Validation: The precision is\",precision_v)\npaste(\"Validation: The recall is\",recall_v)\npaste(\"Validation: The F is\",F_v)\n\nprint(paste0(\"tp = \",tp))\nprint(paste0(\"tn = \",tn))\nprint(paste0(\"fp = \",fp))\nprint(paste0(\"fn = \",fn))\n\nsummary(gbmreg,n.trees = n_gbm_trees, plot=FALSE)[1:30,]\n\nend_time <- Sys.time()\nprint(paste0(\"total run time is \", end_time - start_time, \"mins\"))\n\nmodeloutput[nrow(modeloutput)+1,] <- c(\"DM_M1\", precision, recall, F, precision_v, recall_v, F_v, tp, tn, fp, fn, 'GBM')\nwrite.csv(summary(gbmreg,n.trees = n_gbm_trees, plot=FALSE),paste0(yname,\"_gbm_fi.csv\"))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#random forest\nstart_time <- Sys.time()\n\ntraining[,yname] <- as.factor(training[,yname])\n\nset.seed(123) #reset random seed as gbm may use random subset\nrfreg<-randomForest(f, data=training, mtry = 100, ntree = 500, importance=TRUE)\nrfpredict <- predict(rfreg, training)\npredict <- rfpredict\npredictyes <- sum(predict==1)\npredictno<-sum(predict==0)\nactual <- training[,yname]\nactualwhenpredyes <- actual[predict==1]\nactualwhenpredno <- actual[predict==0]\ntp <- sum(actualwhenpredyes==1)\ntn <- sum(actualwhenpredno==0)\nfp <- length(actualwhenpredyes) - tp\nfn <- length(actualwhenpredno) - tn\nprecision<-tp/predictyes\nrecall<-tp/(tp+fn)\nF = 2*precision*recall/(precision+recall)\npaste(\"Training: The precision is\",precision)\npaste(\"Training: The recall is\",recall)\npaste(\"Training: The F is\",F)\n\nrfpredict <- predict(rfreg, validation)\npredict <- rfpredict\npredictyes <- sum(predict==1)\npredictno<-sum(predict==0)\nactual <- validation[,yname]\nactualwhenpredyes <- actual[predict==1]\nactualwhenpredno <- actual[predict==0]\ntp <- sum(actualwhenpredyes==1)\ntn <- sum(actualwhenpredno==0)\nfp <- length(actualwhenpredyes) - tp\nfn <- length(actualwhenpredno) - tn\nprecision_v<-tp/predictyes\nrecall_v<-tp/(tp+fn)\nF_v = 2*precision_v*recall_v/(precision_v+recall_v)\npaste(\"Validation: The precision is\",precision_v)\npaste(\"Validation: The recall is\",recall_v)\npaste(\"Validation: The F is\",F_v)\n\nprint(paste0(\"tp = \",tp))\nprint(paste0(\"tn = \",tn))\nprint(paste0(\"fp = \",fp))\nprint(paste0(\"fn = \",fn))\n\n#varImp(rfreg)\nvarImpPlot(rfreg,type=2)\n\nend_time <- Sys.time()\nprint(paste0(\"total run time is \", end_time - start_time, \"mins\"))\n\nmodeloutput[nrow(modeloutput)+1,] <- c(\"DM_M1\", precision, recall, F, precision_v, recall_v, F_v, tp, tn, fp, fn, 'randomForest')\nwrite.csv(varImp(rfreg),paste0(yname,\"_rf_fi.csv\"))\n\nwrite.csv(modeloutput,\"injury_predict_model_output.csv\", row.names=FALSE)\n\n#get top variables from gbm and rf\n# n_tops <- 200\n# gbmvars <- as.vector(summary(gbmreg,n.trees = n_gbm_trees, plot=FALSE)[1:n_tops,]$var)\n# rfvarimp <- varImp(rfreg)\n# rfvarimp$var <- rownames(rfvarimp)\n\n# rfvars <- rfvarimp[order(-rfvarimp[,\"1\"]),]\n# rfvars <- as.vector(rfvars$var[1:n_tops])\n\n# dummies <- unique(c(gbmvars, rfvars))\n\n#I have tried using less dummies based on their importance to see if they improve the results\n#However, so far no luck","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Now let's see if we can predict body part, surface, and length of injury based on injury cases\n#Let's do surface first\n\n# Remove columns with constant values\nnew_injury_data <- injury_data[,!apply(injury_data, MARGIN = 2, function(x) max(x, na.rm = TRUE) == min(x, na.rm = TRUE))]\n\ndummies <- dummies[dummies %in% colnames(new_injury_data)]\n\n# Remove one of highly correlated pair\nthreshold <- 0.99\nremove_list <- c()\ncounter <- 2\nfor (a_dummy in dummies[1:(length(dummies)-1)]){\n#\tprint(a_dummy)\n\tfor (b_dummy in dummies[counter:length(dummies)]){\n#\t\tprint(b_dummy)\n\t\tif (abs(cor(new_injury_data[,a_dummy],new_injury_data[,b_dummy])) > threshold){\n\t\t\tremove_list <- c(remove_list, a_dummy)\n\t\t\tprint(paste0(a_dummy, \" will be removed. Correlation between \", a_dummy, \" and \", b_dummy, \" is \", cor(new_injury_data[,a_dummy],new_injury_data[,b_dummy])))\n\t\t\tbreak\n\t\t}\n\t}\n\tcounter <- counter + 1\n}\n\nnew_injury_data <- new_injury_data[, !colnames(new_injury_data) %in% remove_list]\n\ndim(new_injury_data)\n\nynames <- c(\"DM_M1\",\"DM_M7\",\"DM_M28\",\"DM_M42\",\"BodyPart\",\"Surface\",\"PlayerGame\",\"BodyPart_Ankle\",\"BodyPart_Foot\",\"BodyPart_Heel\",\"BodyPart_Knee\",\"BodyPart_Toes\",\"Surface_Natural\",\"Surface_Synthetic\",\"BodyPart\")\nkeys <- c(\"PlayerKey\",\"GameID\",\"PlayKey\")\ncategorical_names <- c(\"RosterPosition\",\"PlayerDay\",\"StadiumType\",\"FieldType\",\"Temperature\",\"Weather\",\"PlayType\",\"PlayerGamePlay\",\"Position\",\"PositionGroup\")\ndummies <- colnames(test_data)[!names(test_data) %in% c(ynames, keys, categorical_names)] \n\ndummies <- dummies[dummies %in% colnames(new_injury_data)]\n\n# Spliting between training data and validation data\nset.seed(6)\nidx <- sample(seq(1, 2), size = nrow(new_injury_data), replace = TRUE, prob = c(.8, .2))\n\ntraining <- new_injury_data[idx==1,]\nvalidation <- new_injury_data[idx==2,]\nprint(dim(training))\nprint(dim(validation))\n\nmodeloutput <- data.frame(y=character(),\n                 Precision_train=double(), #precision based on training data\n                 Recall_train=double(),\t   #recall based on training data\n                 FMeasure_train=double(),  #F measure based on training data\n                 Precision_valid=double(), #precision based on validation data\n                 Recall_valid=double(),    #recall based on validation data\n                 FMeasure_valid=double(),  #F measure based on validation data\n\t\t\t\t tp=double(),              #true positive\n\t\t\t\t tn=double(),\t\t\t   #true negative\n\t\t\t\t fp=double(),\t\t\t   #false positive\n\t\t\t\t fn=double(),\t\t\t   #false negative\n                 Models=character(),\t   #Model type\n                 stringsAsFactors=FALSE)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ynames <- c(\"DM_M7\",\"DM_M28\",\"DM_M42\",\"BodyPart_Ankle\",\"BodyPart_Foot\",\"BodyPart_Heel\",\"BodyPart_Knee\",\"BodyPart_Toes\")\n\nfor (yname in ynames){\n\t#formula\n\tf <-as.formula(paste(yname,\"~\",paste(dummies,collapse=\"+\")))\n\n\t# Logistic Model\n\tstart_time <- Sys.time()\n\n\tglmr <- glm(f, data=training, family=binomial)\n\tglmr <- glm(f, data=training, family=binomial)\n\tglmpredict <- predict(glmr, training)\n\tpredict <- ifelse(glmpredict>0.5, 1, 0)\n\tpredictyes <- sum(predict==1)\n\tpredictno<-sum(predict==0)\n\tactual <- training[,yname]\n\tactualwhenpredyes <- actual[predict==1]\n\tactualwhenpredno <- actual[predict==0]\n\ttp <- sum(actualwhenpredyes==1)\n\ttn <- sum(actualwhenpredno==0)\n\tfp <- length(actualwhenpredyes) - tp\n\tfn <- length(actualwhenpredno) - tn\n\tprecision<-tp/predictyes\n\trecall<-tp/(tp+fn)\n\tF = 2*precision*recall/(precision+recall)\n\tpaste(\"Training: The precision is\",precision)\n\tpaste(\"Training: The recall is\",recall)\n\tpaste(\"Training: The F is\",F)\n\n\tglmpredict <- predict(glmr, validation)\n\tpredict <- ifelse(glmpredict>0.5, 1, 0)\n\tpredictyes <- sum(predict==1)\n\tpredictno<-sum(predict==0)\n\tactual <- validation[,yname]\n\tactualwhenpredyes <- actual[predict==1]\n\tactualwhenpredno <- actual[predict==0]\n\ttp <- sum(actualwhenpredyes==1)\n\ttn <- sum(actualwhenpredno==0)\n\tfp <- length(actualwhenpredyes) - tp\n\tfn <- length(actualwhenpredno) - tn\n\tprecision_v<-tp/predictyes\n\trecall_v<-tp/(tp+fn)\n\tF_v = 2*precision_v*recall_v/(precision_v+recall_v)\n\tpaste(\"Validation: The precision is\",precision_v)\n\tpaste(\"Validation: The recall is\",recall_v)\n\tpaste(\"Validation: The F is\",F_v)\n\n\tprint(paste0(\"tp = \",tp))\n\tprint(paste0(\"tn = \",tn))\n\tprint(paste0(\"fp = \",fp))\n\tprint(paste0(\"fn = \",fn))\n\n\twrite.csv(glmr$coefficients,paste0(yname,'glmr_coef.csv'))\n\tend_time <- Sys.time()\n\tprint(paste0(\"total run time is \", end_time - start_time, \"mins\"))\n\n\tmodeloutput[nrow(modeloutput)+1,] <- c(yname, precision, recall, F, precision_v, recall_v, F_v, tp, tn, fp, fn, 'Logistic')\n\n\n\t#cart\n\tstart_time <- Sys.time()\n\n\tcart = rpart(f, data = training, cp = 10^(-3),minsplit = 10,, method = \"class\")\n\tcartpredict <- predict(cart, training)\n\tpredict <- ifelse(cartpredict[,2]>0.5, 1, 0)\n\tpredictyes <- sum(predict==1)\n\tpredictno<-sum(predict==0)\n\tactual <- training[,yname]\n\tactualwhenpredyes <- actual[predict==1]\n\tactualwhenpredno <- actual[predict==0]\n\ttp <- sum(actualwhenpredyes==1)\n\ttn <- sum(actualwhenpredno==0)\n\tfp <- length(actualwhenpredyes) - tp\n\tfn <- length(actualwhenpredno) - tn\n\tprecision<-tp/predictyes\n\trecall<-tp/(tp+fn)\n\tF = 2*precision*recall/(precision+recall)\n\tpaste(\"Training: The precision is\",precision)\n\tpaste(\"Training: The recall is\",recall)\n\tpaste(\"Training: The F is\",F)\n\n\tcartpredict <- predict(cart, validation)\n\tpredict <- ifelse(cartpredict[,2]>0.5, 1, 0)\n\tpredictyes <- sum(predict==1)\n\tpredictno<-sum(predict==0)\n\tactual <- validation[,yname]\n\tactualwhenpredyes <- actual[predict==1]\n\tactualwhenpredno <- actual[predict==0]\n\ttp <- sum(actualwhenpredyes==1)\n\ttn <- sum(actualwhenpredno==0)\n\tfp <- length(actualwhenpredyes) - tp\n\tfn <- length(actualwhenpredno) - tn\n\tprecision_v<-tp/predictyes\n\trecall_v<-tp/(tp+fn)\n\tF_v = 2*precision_v*recall_v/(precision_v+recall_v)\n\tpaste(\"Validation: The precision is\",precision_v)\n\tpaste(\"Validation: The recall is\",recall_v)\n\tpaste(\"Validation: The F is\",F_v)\n\n\tprint(paste0(\"tp = \",tp))\n\tprint(paste0(\"tn = \",tn))\n\tprint(paste0(\"fp = \",fp))\n\tprint(paste0(\"fn = \",fn))\n\n\tend_time <- Sys.time()\n\tprint(paste0(\"total run time is \", end_time - start_time, \"mins\"))\n\n\tmodeloutput[nrow(modeloutput)+1,] <- c(yname, precision, recall, F, precision_v, recall_v, F_v, tp, tn, fp, fn, 'cart')\n\n\t#knn\n\tstart_time <- Sys.time()\n\n\t# knnreg <- knn.reg(train = training[,names(training) %in% dummies] , test=training[,names(training) %in% dummies], y=training[,yname], k=10, algorithm = \"kd_tree\")\n\t# predict <- ifelse(knnreg$pred>0.5, 1, 0)\n\t# predictyes <- sum(predict==1)\n\t# predictno<-sum(predict==0)\n\t# actual <- training[,yname]\n\t# actualwhenpredyes <- actual[predict==1]\n\t# actualwhenpredno <- actual[predict==0]\n\t# tp <- sum(actualwhenpredyes==1)\n\t# tn <- sum(actualwhenpredno==0)\n\t# fp <- length(actualwhenpredyes) - tp\n\t# fn <- length(actualwhenpredno) - tn\n\t# precision<-tp/predictyes\n\t# recall<-tp/(tp+fn)\n\t# F = 2*precision*recall/(precision+recall)\n\t# paste(\"Training: The precision is\",precision)\n\t# paste(\"Training: The recall is\",recall)\n\t# paste(\"Training: The F is\",F)\n\n\t# knnreg <- knn.reg(train = training[,names(training) %in% dummies] , test=validation[,names(validation) %in% dummies], y=training[,yname], k=10, algorithm = \"kd_tree\")\n\t# predict <- ifelse(knnreg$pred>0.5, 1, 0)\n\t# predictyes <- sum(predict==1)\n\t# predictno<-sum(predict==0)\n\t# actual <- validation[,yname]\n\t# actualwhenpredyes <- actual[predict==1]\n\t# actualwhenpredno <- actual[predict==0]\n\t# tp <- sum(actualwhenpredyes==1)\n\t# tn <- sum(actualwhenpredno==0)\n\t# fp <- length(actualwhenpredyes) - tp\n\t# fn <- length(actualwhenpredno) - tn\n\t# precision_v<-tp/predictyes\n\t# recall_v<-tp/(tp+fn)\n\t# F_v = 2*precision_v*recall_v/(precision_v+recall_v)\n\t# paste(\"Validation: The precision is\",precision_v)\n\t# paste(\"Validation: The recall is\",recall_v)\n\t# paste(\"Validation: The F is\",F_v)\n\n\t# print(paste0(\"tp = \",tp))\n\t# print(paste0(\"tn = \",tn))\n\t# print(paste0(\"fp = \",fp))\n\t# print(paste0(\"fn = \",fn))\n\n\t# end_time <- Sys.time()\n\t# print(paste0(\"total run time is \", end_time - start_time, \"mins\"))\n\n\t# modeloutput[nrow(modeloutput)+1,] <- c(yname, precision, recall, F, precision_v, recall_v, F_v, tp, tn, fp, fn, 'knn')\n\n\t#linear regression\n\tstart_time <- Sys.time()\n\n\tlr<-lm(f, data=training)\n\tpredict <- ifelse(lr$fitted.values>0.5, 1, 0)\n\tpredictyes <- sum(predict==1)\n\tpredictno<-sum(predict==0)\n\tactual <- training[,yname]\n\tactualwhenpredyes <- actual[predict==1]\n\tactualwhenpredno <- actual[predict==0]\n\ttp <- sum(actualwhenpredyes==1)\n\ttn <- sum(actualwhenpredno==0)\n\tfp <- length(actualwhenpredyes) - tp\n\tfn <- length(actualwhenpredno) - tn\n\tprecision<-tp/predictyes\n\trecall<-tp/(tp+fn)\n\tF = 2*precision*recall/(precision+recall)\n\tpaste(\"Training: The precision is\",precision)\n\tpaste(\"Training: The recall is\",recall)\n\tpaste(\"Training: The F is\",F)\n\n\tpredict <- ifelse(predict(lr,validation)>0.5, 1, 0)\n\tpredictyes <- sum(predict==1)\n\tpredictno<-sum(predict==0)\n\tactual <- validation[,yname]\n\tactualwhenpredyes <- actual[predict==1]\n\tactualwhenpredno <- actual[predict==0]\n\ttp <- sum(actualwhenpredyes==1)\n\ttn <- sum(actualwhenpredno==0)\n\tfp <- length(actualwhenpredyes) - tp\n\tfn <- length(actualwhenpredno) - tn\n\tprecision_v<-tp/predictyes\n\trecall_v<-tp/(tp+fn)\n\tF_v = 2*precision_v*recall_v/(precision_v+recall_v)\n\tpaste(\"Validation: The precision is\",precision_v)\n\tpaste(\"Validation: The recall is\",recall_v)\n\tpaste(\"Validation: The F is\",F_v)\n\n\tprint(paste0(\"tp = \",tp))\n\tprint(paste0(\"tn = \",tn))\n\tprint(paste0(\"fp = \",fp))\n\tprint(paste0(\"fn = \",fn))\n\n\tend_time <- Sys.time()\n\tprint(paste0(\"total run time is \", end_time - start_time, \"mins\"))\n\n\tmodeloutput[nrow(modeloutput)+1,] <- c(yname, precision, recall, F, precision_v, recall_v, F_v, tp, tn, fp, fn, 'linear')\n\n\t# #lasso\n\t# start_time <- Sys.time()\n\t# lambda <- 0.5\n\t# lasso<-glmnet(as.matrix(training[,dummies]), as.matrix(training[,yname]), alpha = 1, lambda = lambda)\n\t# cv.out <- cv.glmnet(as.matrix(training[,dummies]), as.matrix(training[,yname]), alpha = 1)\n\t# lambda <- cv.out$lambda.min\n\t# lassopredict <- predict(lasso, s = lambda, newx = as.matrix(training[dummies]))\n\t# predict <- ifelse(lassopredict>0.5, 1, 0)\n\t# predictyes <- sum(predict==1)\n\t# predictno<-sum(predict==0)\n\t# actual <- training[,yname]\n\t# actualwhenpredyes <- actual[predict==1]\n\t# actualwhenpredno <- actual[predict==0]\n\t# tp <- sum(actualwhenpredyes==1)\n\t# tn <- sum(actualwhenpredno==0)\n\t# fp <- length(actualwhenpredyes) - tp\n\t# fn <- length(actualwhenpredno) - tn\n\t# precision<-tp/predictyes\n\t# recall<-tp/(tp+fn)\n\t# F = 2*precision*recall/(precision+recall)\n\t# paste(\"Training: The precision is\",precision)\n\t# paste(\"Training: The recall is\",recall)\n\t# paste(\"Training: The F is\",F)\n\n\t# lassopredict <- predict(lasso, s = lambda, newx = as.matrix(validation[dummies]))\n\t# predict <- ifelse(lassopredict>0.5, 1, 0)\n\t# predictyes <- sum(predict==1)\n\t# predictno<-sum(predict==0)\n\t# actual <- validation[,yname]\n\t# actualwhenpredyes <- actual[predict==1]\n\t# actualwhenpredno <- actual[predict==0]\n\t# tp <- sum(actualwhenpredyes==1)\n\t# tn <- sum(actualwhenpredno==0)\n\t# fp <- length(actualwhenpredyes) - tp\n\t# fn <- length(actualwhenpredno) - tn\n\t# precision_v<-tp/predictyes\n\t# recall_v<-tp/(tp+fn)\n\t# F_v = 2*precision_v*recall_v/(precision_v+recall_v)\n\t# paste(\"Validation: The precision is\",precision_v)\n\t# paste(\"Validation: The recall is\",recall_v)\n\t# paste(\"Validation: The F is\",F_v)\n\n\t# print(paste0(\"tp = \",tp))\n\t# print(paste0(\"tn = \",tn))\n\t# print(paste0(\"fp = \",fp))\n\t# print(paste0(\"fn = \",fn))\n\n\t# end_time <- Sys.time()\n\t# print(paste0(\"total run time is \", end_time - start_time, \"mins\"))\n\n\t# modeloutput[nrow(modeloutput)+1,] <- c(yname, precision, recall, F, precision_v, recall_v, F_v, tp, tn, fp, fn, 'lasso')\n\n\t# #ridge\n\t# start_time <- Sys.time()\n\n\t# ridge<-glmnet(as.matrix(training[,dummies]), as.matrix(training[,yname]), alpha = 0, lambda = lambda)\n\t# cv.out <- cv.glmnet(as.matrix(training[,dummies]), as.matrix(training[,yname]), alpha = 0)\n\t# lambda <- cv.out$lambda.min\n\t# ridgepredict <- predict(ridge, s = lambda, newx = as.matrix(training[dummies]))\n\t# predict <- ifelse(ridgepredict>0.5, 1, 0)\n\t# predictyes <- sum(predict==1)\n\t# predictno<-sum(predict==0)\n\t# actual <- training[,yname]\n\t# actualwhenpredyes <- actual[predict==1]\n\t# actualwhenpredno <- actual[predict==0]\n\t# tp <- sum(actualwhenpredyes==1)\n\t# tn <- sum(actualwhenpredno==0)\n\t# fp <- length(actualwhenpredyes) - tp\n\t# fn <- length(actualwhenpredno) - tn\n\t# precision<-tp/predictyes\n\t# recall<-tp/(tp+fn)\n\t# F = 2*precision*recall/(precision+recall)\n\t# paste(\"Training: The precision is\",precision)\n\t# paste(\"Training: The recall is\",recall)\n\t# paste(\"Training: The F is\",F)\n\n\t# ridgepredict <- predict(ridge, s = lambda, newx = as.matrix(validation[dummies]))\n\t# predict <- ifelse(ridgepredict>0.5, 1, 0)\n\t# predictyes <- sum(predict==1)\n\t# predictno<-sum(predict==0)\n\t# actual <- validation[,yname]\n\t# actualwhenpredyes <- actual[predict==1]\n\t# actualwhenpredno <- actual[predict==0]\n\t# tp <- sum(actualwhenpredyes==1)\n\t# tn <- sum(actualwhenpredno==0)\n\t# fp <- length(actualwhenpredyes) - tp\n\t# fn <- length(actualwhenpredno) - tn\n\t# precision_v<-tp/predictyes\n\t# recall_v<-tp/(tp+fn)\n\t# F_v = 2*precision_v*recall_v/(precision_v+recall_v)\n\t# paste(\"Validation: The precision is\",precision_v)\n\t# paste(\"Validation: The recall is\",recall_v)\n\t# paste(\"Validation: The F is\",F_v)\n\n\t# print(paste0(\"tp = \",tp))\n\t# print(paste0(\"tn = \",tn))\n\t# print(paste0(\"fp = \",fp))\n\t# print(paste0(\"fn = \",fn))\n\n\t# end_time <- Sys.time()\n\t# print(paste0(\"total run time is \", end_time - start_time, \"mins\"))\n\n\t# modeloutput[nrow(modeloutput)+1,] <- c(yname, precision, recall, F, precision_v, recall_v, F_v, tp, tn, fp, fn, 'ridge')\n\n\t#gbm\n\tstart_time <- Sys.time()\n\n\tn_gbm_trees = 500\n\n\tset.seed(123) #reset random seed as gbm may use random subset\n\tgbmreg<-gbm(f, data=training, distribution = \"bernoulli\", interaction.depth=10, n.minobsinnode = 1, bag.fraction=0.7, n.trees = n_gbm_trees, shrinkage = 0.1)\n\tgbmpredict <- predict(gbmreg, training, n.trees = n_gbm_trees)\n\tpredict <- ifelse(gbmpredict>0.5, 1, 0)\n\tpredictyes <- sum(predict==1)\n\tpredictno<-sum(predict==0)\n\tactual <- training[,yname]\n\tactualwhenpredyes <- actual[predict==1]\n\tactualwhenpredno <- actual[predict==0]\n\ttp <- sum(actualwhenpredyes==1)\n\ttn <- sum(actualwhenpredno==0)\n\tfp <- length(actualwhenpredyes) - tp\n\tfn <- length(actualwhenpredno) - tn\n\tprecision<-tp/predictyes\n\trecall<-tp/(tp+fn)\n\tF = 2*precision*recall/(precision+recall)\n\tpaste(\"Training: The precision is\",precision)\n\tpaste(\"Training: The recall is\",recall)\n\tpaste(\"Training: The F is\",F)\n\n\tgbmpredict <- predict(gbmreg, validation, n.trees = n_gbm_trees)\n\tpredict <- ifelse(gbmpredict>0.5, 1, 0)\n\tpredictyes <- sum(predict==1)\n\tpredictno<-sum(predict==0)\n\tactual <- validation[,yname]\n\tactualwhenpredyes <- actual[predict==1]\n\tactualwhenpredno <- actual[predict==0]\n\ttp <- sum(actualwhenpredyes==1)\n\ttn <- sum(actualwhenpredno==0)\n\tfp <- length(actualwhenpredyes) - tp\n\tfn <- length(actualwhenpredno) - tn\n\tprecision_v<-tp/predictyes\n\trecall_v<-tp/(tp+fn)\n\tF_v = 2*precision_v*recall_v/(precision_v+recall_v)\n\tpaste(\"Validation: The precision is\",precision_v)\n\tpaste(\"Validation: The recall is\",recall_v)\n\tpaste(\"Validation: The F is\",F_v)\n\n\tprint(paste0(\"tp = \",tp))\n\tprint(paste0(\"tn = \",tn))\n\tprint(paste0(\"fp = \",fp))\n\tprint(paste0(\"fn = \",fn))\n\n\tsummary(gbmreg,n.trees = n_gbm_trees, plot=FALSE)[1:30,]\n\twrite.csv(summary(gbmreg,n.trees = n_gbm_trees, plot=FALSE),paste0(yname,\"_gbm_fi.csv\"))\n\n\tend_time <- Sys.time()\n\tprint(paste0(\"total run time is \", end_time - start_time, \"mins\"))\n\n\tmodeloutput[nrow(modeloutput)+1,] <- c(yname, precision, recall, F, precision_v, recall_v, F_v, tp, tn, fp, fn, 'gbm')\n\n\t#random forest\n\tstart_time <- Sys.time()\n\n\ttraining[,yname] <- as.factor(training[,yname])\n\n\tset.seed(123) #reset random seed as gbm may use random subset\n\trfreg<-randomForest(f, data=training, mtry = 100, ntree = 500, importance=TRUE)\n\trfpredict <- predict(rfreg, training)\n\tpredict <- rfpredict\n\tpredictyes <- sum(predict==1)\n\tpredictno<-sum(predict==0)\n\tactual <- training[,yname]\n\tactualwhenpredyes <- actual[predict==1]\n\tactualwhenpredno <- actual[predict==0]\n\ttp <- sum(actualwhenpredyes==1)\n\ttn <- sum(actualwhenpredno==0)\n\tfp <- length(actualwhenpredyes) - tp\n\tfn <- length(actualwhenpredno) - tn\n\tprecision<-tp/predictyes\n\trecall<-tp/(tp+fn)\n\tF = 2*precision*recall/(precision+recall)\n\tpaste(\"Training: The precision is\",precision)\n\tpaste(\"Training: The recall is\",recall)\n\tpaste(\"Training: The F is\",F)\n\n\trfpredict <- predict(rfreg, validation)\n\tpredict <- rfpredict\n\tpredictyes <- sum(predict==1)\n\tpredictno<-sum(predict==0)\n\tactual <- validation[,yname]\n\tactualwhenpredyes <- actual[predict==1]\n\tactualwhenpredno <- actual[predict==0]\n\ttp <- sum(actualwhenpredyes==1)\n\ttn <- sum(actualwhenpredno==0)\n\tfp <- length(actualwhenpredyes) - tp\n\tfn <- length(actualwhenpredno) - tn\n\tprecision_v<-tp/predictyes\n\trecall_v<-tp/(tp+fn)\n\tF_v = 2*precision_v*recall_v/(precision_v+recall_v)\n\tpaste(\"Validation: The precision is\",precision_v)\n\tpaste(\"Validation: The recall is\",recall_v)\n\tpaste(\"Validation: The F is\",F_v)\n\n\tprint(paste0(\"tp = \",tp))\n\tprint(paste0(\"tn = \",tn))\n\tprint(paste0(\"fp = \",fp))\n\tprint(paste0(\"fn = \",fn))\n\n\twrite.csv(varImp(rfreg),paste0(yname,\"_rf_fi.csv\"))\n\tvarImpPlot(rfreg,type=2)\n\n\tend_time <- Sys.time()\n\tprint(paste0(\"total run time is \", end_time - start_time, \"mins\"))\n\n\tmodeloutput[nrow(modeloutput)+1,] <- c(yname, precision, recall, F, precision_v, recall_v, F_v, tp, tn, fp, fn, 'randomForest')\n}\n\nwrite.csv(modeloutput,\"injury_model_output.csv\", row.names=FALSE)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"R","language":"R","name":"ir"},"language_info":{"mimetype":"text/x-r-source","name":"R","pygments_lexer":"r","version":"3.4.2","file_extension":".r","codemirror_mode":"r"}},"nbformat":4,"nbformat_minor":1}