{"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":"code","source":"# Basic Logit Submission\n# Working with feather data, but want to see how far basic tools can get me\n# Still have to be very careful with memory, can't really handle test data in memory (even as feather)\n# Started with: https://www.kaggle.com/code/stautxie/amex-default-prediction-model-in-r-part-1\n\nsuppressPackageStartupMessages(library(data.table)) \nsuppressPackageStartupMessages(library(tidyverse))\nsuppressPackageStartupMessages(library(dtplyr)) #data.table with tidy syntax\nsuppressPackageStartupMessages(library(arrow))\n\ndir(\"..\")\nprint('available files...')\nlist.files(path = \"../input/amex-default-prediction\") %>% print()\n\npqt_dir <- '../input/amex-data-integer-dtypes-parquet-format'\ncsv_dir <- '../input/amex-default-prediction'\ndt_threads <- getDTthreads()\ncat(paste(\"Number of threads for data.table: \", dt_threads, \"\\n\", sep=\"\"))","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2022-06-24T17:50:23.008409Z","iopub.execute_input":"2022-06-24T17:50:23.011137Z","iopub.status.idle":"2022-06-24T17:50:25.010032Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#pull in training labels\nt1 <- Sys.time()\ntrain_Y <- fread(\"../input/amex-default-prediction/train_labels.csv\") %>% as_tibble()\nt2 <- Sys.time()\n\nprint('Time to load training labels..')\ndifftime(t2,t1, units=\"secs\")\nprint('Number of rows')\ntrain_Y %>% nrow()\nprint('Number of IDs')\ntrain_Y %>% distinct(customer_ID) %>% nrow()","metadata":{"execution":{"iopub.status.busy":"2022-06-24T17:50:25.013133Z","iopub.execute_input":"2022-06-24T17:50:25.05059Z","iopub.status.idle":"2022-06-24T17:50:28.470199Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#load subset of the efficiently-stored training data, collapse to 1 row customer on import\nt1 <- Sys.time()\n\ntrain_df <- \n    arrow::read_parquet(file.path(pqt_dir, \"train.parquet\"), col_select = 1:10) %>% \n    mutate(S_2 = lubridate::ymd(S_2)) %>%\n    group_by(customer_ID) %>% \n    slice_max(S_2) %>% \n    ungroup()\n\nt2 <- Sys.time()\n\nprint('Time to load training parquet file with arrow and grab latest obs per group')\ndifftime(t2,t1, units=\"secs\")","metadata":{"execution":{"iopub.status.busy":"2022-06-24T17:50:28.474619Z","iopub.execute_input":"2022-06-24T17:50:28.477213Z","iopub.status.idle":"2022-06-24T17:53:42.035927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#attach training labels to training data\ntrain_df <- left_join(train_Y,train_df,by=c(\"customer_ID\"))\n\nnrow(train_df)\ncolSums(is.na(train_df))\n\nprint('Proportion of defaults')\ntrain_df %>% count(target) %>% mutate(perc = n/sum(n)) %>% print()","metadata":{"execution":{"iopub.status.busy":"2022-06-24T17:53:42.038933Z","iopub.execute_input":"2022-06-24T17:53:42.040757Z","iopub.status.idle":"2022-06-24T17:53:42.555243Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#estimate logit on the first few variables\nmodel1 = glm(target ~ P_2 + B_1 + B_2 + R_1 + D_39 + S_3, data=train_df, family=\"binomial\")\nsummary(model1)\n\n#predict\npredictions = predict(model1, newdata=train_df, type=\"response\") %>% replace_na(0) #fallback should definitely be a 0\nhead(predictions)","metadata":{"execution":{"iopub.status.busy":"2022-06-24T17:53:42.557442Z","iopub.execute_input":"2022-06-24T17:53:42.558732Z","iopub.status.idle":"2022-06-24T17:53:45.337871Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#visualize how competition metric varies with prediction threshold\n    #might be interesting to do this with accuracy and specificity as well\n\n# thank you: https://www.kaggle.com/code/igjit1/r-package-for-amex-competition-metric\nremotes::install_github(\"igjit/amexmetric\", quiet = TRUE)\n\n#calc loss for different thresholds\nthreshold = seq(0.1,0.9,0.1)\nmetric=numeric()\n\nfor (i in 1:length(threshold)){\n    metric[i] = amexmetric::amex_metric( pull(train_Y,target), (predictions > threshold[i]))\n}\n\n#viz loss for different thresholds\ntibble(threshold,metric) %>%\n    ggplot(aes(x=threshold,y=metric))+\n    geom_point()+\n    geom_line()+\n    geom_text(label=round(metric,3),vjust=-1)\n\n#best threshold\nbest_thresh = tibble(threshold,metric) %>% slice_max(metric) %>% select(threshold) %>% pull()","metadata":{"execution":{"iopub.status.busy":"2022-06-24T18:20:35.060754Z","iopub.execute_input":"2022-06-24T18:20:35.062689Z","iopub.status.idle":"2022-06-24T18:20:39.804891Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#tidy up workspace (hopefully will be able to keep the model and the test data in memory at the same time...)\nrm(train_df,train_Y,predictions)\ngc()","metadata":{"execution":{"iopub.status.busy":"2022-06-24T15:57:27.067927Z","iopub.execute_input":"2022-06-24T15:57:27.07003Z","iopub.status.idle":"2022-06-24T15:57:27.521564Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#load the efficiently-stored test data, collapse to 1 row per customer on import\nt1 <- Sys.time()\n\ntest_X <- \n    arrow::read_parquet(file.path(pqt_dir, \"test.parquet\"), col_select = 1:10) %>% \n    mutate(S_2 = lubridate::ymd(S_2)) %>%\n    group_by(customer_ID) %>% \n    slice_max(S_2) %>% \n    ungroup()\n\nt2 <- Sys.time()\n\nprint('Time to load test parquet file with arrow and grab latest obs per group')\ndifftime(t2,t1, units=\"secs\")\n\nclass(test_X)\nnrow(test_X)\n#test_X","metadata":{"execution":{"iopub.status.busy":"2022-06-24T15:57:27.524156Z","iopub.execute_input":"2022-06-24T15:57:27.52561Z","iopub.status.idle":"2022-06-24T16:04:22.780673Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#make predictions on test data\n\n#predict\npredictions = predict(model1, newdata=test_X, type=\"response\") %>% replace_na(0) #fallback should definitely be a 0\nhead(predictions)\n\n#check distribution\nhist(predictions)","metadata":{"execution":{"iopub.status.busy":"2022-06-24T16:10:08.348447Z","iopub.execute_input":"2022-06-24T16:10:08.351772Z","iopub.status.idle":"2022-06-24T16:10:10.601194Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#bind test_IDs and predictions\ntest_df <- test_X %>% mutate(prediction=as.numeric((predictions>best_thresh))) %>% select(customer_ID,prediction)\nnrow(test_df)\nhead(test_df)","metadata":{"execution":{"iopub.status.busy":"2022-06-24T16:15:12.22818Z","iopub.execute_input":"2022-06-24T16:15:12.230243Z","iopub.status.idle":"2022-06-24T16:15:12.282829Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#write to submission.csv \nfwrite(x=test_df, file=\"submission.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-06-24T16:15:44.773968Z","iopub.execute_input":"2022-06-24T16:15:44.775885Z","iopub.status.idle":"2022-06-24T16:15:44.917477Z"},"trusted":true},"execution_count":null,"outputs":[]}]}