{"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":"# Use Lasso find a good subset of features\n# Working with feather data, but still have to be very careful with memory\n# Following: 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))\nsuppressPackageStartupMessages(library(glmnet))\nsuppressPackageStartupMessages(library(broom))\nsuppressPackageStartupMessages(library(corrplot))\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-30T20:28:06.32078Z","iopub.execute_input":"2022-06-30T20:28:06.32328Z","iopub.status.idle":"2022-06-30T20:28:12.469671Z"},"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()\nprint(paste('columns in target data..',colnames(train_Y)))","metadata":{"execution":{"iopub.status.busy":"2022-06-30T20:28:12.483115Z","iopub.execute_input":"2022-06-30T20:28:12.628893Z","iopub.status.idle":"2022-06-30T20:28:15.240717Z"},"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_X <- \n    arrow::read_parquet(file.path(pqt_dir, \"train.parquet\")) %>% \n    mutate(S_2 = lubridate::ymd(S_2)) %>%\n    #\n    #grab latest obs for each group\n    group_by(customer_ID) %>% \n    slice_max(S_2) %>% \n    ungroup()\n\n#do calculations taking 1 row per customer\ntrain_X <- \n    train_X %>%\n    #\n    mutate(year = lubridate::year(S_2)) %>%\n    mutate(month = lubridate::month(S_2)) %>% \n    #\n    #match to target vector\n    right_join(., train_Y, by=\"customer_ID\") %>%\n    #\n    #downsample due to size issues\n    slice_sample(prop=0.5)\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\")\n\n#this is just being super careful to keep the ordering right\ntrain_Y <- train_X %>% select(customer_ID,target)\ntrain_X <- train_X %>% select(-target)","metadata":{"execution":{"iopub.status.busy":"2022-06-30T20:28:15.245733Z","iopub.execute_input":"2022-06-30T20:28:15.248189Z","iopub.status.idle":"2022-06-30T20:32:13.744104Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# diagnostics on loading\n\n# drop columns with too many missings\ntrain_X = train_X[, which(colMeans(is.na(train_X)) < 0.25)]\n\n#check\nprint('percentage of missings...')\ncolSums(is.na(train_X))/nrow(train_X)","metadata":{"execution":{"iopub.status.busy":"2022-06-30T20:32:13.748554Z","iopub.execute_input":"2022-06-30T20:32:13.750734Z","iopub.status.idle":"2022-06-30T20:32:14.473065Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#replace all nas with [fill in blank]\n\nnames <- colnames(train_X) \n\nfor (j in 3:ncol(train_X)){ \n    \n    name <- names[j]    \n    #\n    col <- train_X %>% select(., !!name) %>% pull()    \n    impute_val <- -9\n    new_col <- ifelse(is.na(col),impute_val,col)\n    #    \n    train_X <- train_X %>% mutate( !!name := new_col )\n    \n}\n\n#check\nprint('there should be no missings anymore')\nsum(is.na(train_X))","metadata":{"execution":{"iopub.status.busy":"2022-06-30T20:32:14.476739Z","iopub.execute_input":"2022-06-30T20:32:14.479016Z","iopub.status.idle":"2022-06-30T20:32:17.599304Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#matrix transformations (if you're working with this)\ntrain_X <- as.matrix(train_X[,3:ncol(train_X)])\ntrain_Y <- as.factor(train_Y$target)\nclass(train_X)\nclass(train_Y)","metadata":{"execution":{"iopub.status.busy":"2022-06-30T20:32:17.603173Z","iopub.execute_input":"2022-06-30T20:32:17.605482Z","iopub.status.idle":"2022-06-30T20:32:17.911877Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#let's keep this model pretty simple and see how it goes\n\n#lasso model\nt1 <- Sys.time()\n\nlambda_grid=10^seq(5,-2,length=25)\nmod.cv.glmnet <- cv.glmnet(train_X, train_Y, family = \"binomial\", alpha=1, lambda=lambda_grid, nfolds=10) \n\nt2 <- Sys.time()\n\nprint('Time to fit baseline CV lasso model')\ndifftime(t2,t1, units=\"secs\")\n","metadata":{"execution":{"iopub.status.busy":"2022-06-30T20:32:17.914427Z","iopub.execute_input":"2022-06-30T20:32:17.916286Z","iopub.status.idle":"2022-06-30T20:33:44.604523Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#have a look at CV (looks absolutely horrible... basically just shrink penalty to very small)\nplot(mod.cv.glmnet)\n\n#get best lambda\nbest_lambda <- mod.cv.glmnet$lambda.min\nprint(paste('optimal lambda value', best_lambda))","metadata":{"execution":{"iopub.status.busy":"2022-06-30T20:33:44.607148Z","iopub.execute_input":"2022-06-30T20:33:44.609117Z","iopub.status.idle":"2022-06-30T20:33:44.923001Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#recover the coefficients that LASSO preserves\nmod.glmnet = glmnet(train_X, train_Y, family = \"binomial\", alpha=1, lambda=best_lambda, intercept = F)\ncoefs = predict(mod.glmnet,type=\"coefficients\")\nnonzero_indices = which(coefs!=0)","metadata":{"execution":{"iopub.status.busy":"2022-06-30T20:33:44.925976Z","iopub.execute_input":"2022-06-30T20:33:44.927872Z","iopub.status.idle":"2022-06-30T20:34:26.129377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#coefs (forget about the intercept)\nnonzero_indices = nonzero_indices-1","metadata":{"execution":{"iopub.status.busy":"2022-06-30T20:34:26.132431Z","iopub.execute_input":"2022-06-30T20:34:26.134071Z","iopub.status.idle":"2022-06-30T20:34:26.148007Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#pull those variables out\nselected_vars <- train_X[,nonzero_indices] %>% as_tibble() %>% colnames()\nprint(selected_vars)\n\n#pull those variables out\nnew_df = tibble( target=as.numeric(train_Y)-1, as_tibble(train_X[,nonzero_indices]) )\nhead(new_df)\n\n#clean up\nrm(train_X, train_Y)\ngc()","metadata":{"execution":{"iopub.status.busy":"2022-07-01T16:34:51.530421Z","iopub.execute_input":"2022-07-01T16:34:51.534509Z","iopub.status.idle":"2022-07-01T16:34:51.760153Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#estimate logit\nmodel1 = glm(target ~ ., data=new_df, family=\"binomial\") #summary(model1)\n\n#tidy display, maybe sort or whatever\ntidy(model1) %>%\n    mutate(estimate = round(estimate, digits=3),\n           std.error = round(std.error, digits=3),\n           statistic = round(statistic, digits=3),\n           p.value = round(p.value, digits=3))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#export for later use\nselected_vars_tibble = tibble(vars = selected_vars)\nprint(selected_vars_tibble, n=Inf)\nwrite_csv(selected_vars_tibble, \"lasso_vars.csv\")","metadata":{},"execution_count":null,"outputs":[]}]}