{"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))\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-07-13T23:51:57.837609Z","iopub.execute_input":"2022-07-13T23:51:57.879718Z","iopub.status.idle":"2022-07-13T23:51:57.934993Z"},"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-07-13T23:51:57.938190Z","iopub.execute_input":"2022-07-13T23:51:57.939848Z","iopub.status.idle":"2022-07-13T23:51:59.717573Z"},"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    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-07-13T23:51:59.719822Z","iopub.execute_input":"2022-07-13T23:51:59.721297Z","iopub.status.idle":"2022-07-13T23:55:27.644079Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# this is just being super careful to keep the ordering right -- attach training labels to training data\ntrain_X <- left_join(train_Y,train_X,by=c(\"customer_ID\"))\n# separate X,Y\ntrain_Y <- train_X %>% select(customer_ID,target)\ntrain_X <- train_X %>% select(-target)","metadata":{"execution":{"iopub.status.busy":"2022-07-13T23:55:27.646316Z","iopub.execute_input":"2022-07-13T23:55:27.647649Z","iopub.status.idle":"2022-07-13T23:55:28.789332Z"},"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-07-13T23:55:28.791695Z","iopub.execute_input":"2022-07-13T23:55:28.793088Z","iopub.status.idle":"2022-07-13T23:55:30.802276Z"},"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 <- median(col, na.rm=T)\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-07-13T23:55:30.804641Z","iopub.execute_input":"2022-07-13T23:55:30.806066Z","iopub.status.idle":"2022-07-13T23:55:37.495511Z"},"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-07-13T23:55:37.498086Z","iopub.execute_input":"2022-07-13T23:55:37.499566Z","iopub.status.idle":"2022-07-13T23:55:38.180192Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gc()","metadata":{"execution":{"iopub.status.busy":"2022-07-14T00:16:29.931017Z","iopub.execute_input":"2022-07-14T00:16:29.932770Z","iopub.status.idle":"2022-07-14T00:16:30.601524Z"},"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(3,-2,length=10)\nmod.cv.glmnet <- cv.glmnet(train_X, train_Y, family = \"binomial\", alpha=1, lambda=lambda_grid, nfolds=5) \n\nt2 <- Sys.time()\n\nprint('Time to fit baseline CV lasso model')\ndifftime(t2,t1, units=\"secs\")","metadata":{"execution":{"iopub.status.busy":"2022-07-13T23:55:38.927236Z","iopub.execute_input":"2022-07-13T23:55:38.928821Z","iopub.status.idle":"2022-07-13T23:56:28.850539Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#have a look at CV (...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-07-13T23:56:28.853855Z","iopub.execute_input":"2022-07-13T23:56:28.855639Z","iopub.status.idle":"2022-07-13T23:56:29.181812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#\n# placeholder.\n#","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"seq(0.1*best_lambda,0.1, length=10)","metadata":{"execution":{"iopub.status.busy":"2022-07-14T00:22:33.637004Z","iopub.execute_input":"2022-07-14T00:22:33.640765Z","iopub.status.idle":"2022-07-14T00:22:33.674855Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#iterate over some lambdas and find how many coeficients are preserved\nlambda_eps = seq(0.1*best_lambda,0.1, length=10)\n\n#vector of no. of coefs, list of coefs, list of models\nncoef_vec = numeric()\ncoef_list = list()\nmodel_list = list()\n\n\n\n#idk if this will help for reproducibility\nset.seed(1)\n\n\n\nt1 <- Sys.time()\n\nfor (i in 1:length(lambda_eps)){\n    \n    #fit model\n    mod = glmnet( train_X, train_Y, family = \"binomial\", alpha=1, lambda=lambda_eps[i] )\n    \n    #get nonzero coefs\n    coefs = predict(mod,type=\"coefficients\") \n    nonzero_indices = which(coefs!=0)\n    selected_vars <- colnames(train_X)[nonzero_indices-1] \n    \n    #put num-coefs in vector\n    ncoef_vec[i] = length(nonzero_indices)\n    \n    #put variables selected in list\n    coef_list[[i]] = selected_vars\n    \n    #put model in list\n    model_list[[i]] = mod\n        \n}\n\nt2 <- Sys.time()\n\nprint('Time to fit several lasso models')\ndifftime(t2,t1, units=\"secs\")\n\n\n\n#viz num coeffs for diff lambda vals\ntibble(lambda = lambda_eps, ncoef_vec) %>%\n    ggplot(aes(x=lambda, y=ncoef_vec))+\n    geom_point()+\n    geom_line()+\n    geom_text(label=ncoef_vec, vjust=-1)+\n    theme(text = element_text(size=25))","metadata":{"execution":{"iopub.status.busy":"2022-07-14T00:20:06.963706Z","iopub.execute_input":"2022-07-14T00:20:06.965830Z","iopub.status.idle":"2022-07-14T00:21:27.071447Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gc()","metadata":{"execution":{"iopub.status.busy":"2022-07-14T00:23:44.452600Z","iopub.execute_input":"2022-07-14T00:23:44.454996Z","iopub.status.idle":"2022-07-14T00:23:44.784227Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# \n#\n#\n# PLACEHOLDER!\n#\n#\n#\n\n### check on outputs ###\n\n#class(ncoef_vec)\n#class(model_list)\n#class(coef_list)\n\n#best_lambda+eps\n#ncoef_vec\n#coef_list\n#model_list","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#pull out models that add information\ncoef_list\n#coef_list = coef_list[c(1,2,3,4,5,6,7,10,14)]\n#coef_list","metadata":{"execution":{"iopub.status.busy":"2022-07-14T00:23:49.538899Z","iopub.execute_input":"2022-07-14T00:23:49.540518Z","iopub.status.idle":"2022-07-14T00:23:49.564784Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#set up train_df to estimate/validate on\nnew_df = tibble( target=as.numeric(train_Y)-1, as_tibble(train_X) ) %>% mutate( splitter = runif(n=n()) )\nhead(new_df)\n\n#clean up\nrm(train_X, train_Y)\ngc()\n\n#set up train_df to estimate on\ntrain_df = new_df %>% filter(splitter < 0.7)\nprint('nrow, ncol in training data')\nnrow(train_df)\nncol(train_df)\n\n#set up test_df to validate on\ntest_df = new_df %>% filter(splitter >= 0.7)\nprint('nrow, ncol in test data')\nnrow(test_df)\nncol(test_df)","metadata":{"execution":{"iopub.status.busy":"2022-07-14T00:40:16.630928Z","iopub.execute_input":"2022-07-14T00:40:16.633893Z","iopub.status.idle":"2022-07-14T00:40:18.114257Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# test metric -- thank you: https://www.kaggle.com/code/igjit1/r-package-for-amex-competition-metric\nremotes::install_github(\"igjit/amexmetric\", quiet = TRUE)","metadata":{"execution":{"iopub.status.busy":"2022-07-14T00:40:23.388849Z","iopub.execute_input":"2022-07-14T00:40:23.390459Z","iopub.status.idle":"2022-07-14T00:40:58.492845Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#loop through list and fit models\n\n\n#store in-sample loss for each model\nloss_vec = numeric()\ncoef_len_vec = numeric()\n\n\n# loop thru list of models' coefficients\nfor (m in 1:length(coef_list)){\n    \n    #estimate model\n    coefs_m <- coef_list[[m]]\n    model_m = glm(target ~ ., data=(train_df %>% select(target, coefs_m)), family=\"binomial\") #summary(model_m)\n    \n    #predict\n    predictions = predict(model_m, newdata=(test_df %>% select(target, coefs_m)), type=\"response\") %>% replace_na(0) #fallback should be a 0\n    #head(predictions)\n    \n    #calc loss -- this should be done with CV unless you think it's not overfit\n    loss_m <- amexmetric::amex_metric( pull(test_df,target), predictions)\n    \n    #output\n    print(paste(\"Model loss equal to:\"))\n    print(loss_m)\n    print(\"estimated with coefs\")\n    print(coefs_m)\n    print(\"-----------------\")\n    \n    #store results (this is slightly redundant)\n    loss_vec[m] = loss_m\n    coef_len_vec[m] = length(coefs_m)\n    \n}\n\n\n#viz num coeffs for diff lambda vals\ntibble(mod_size=coef_len_vec, loss = loss_vec) %>%\n    ggplot(aes(x=coef_len_vec, y=loss_vec))+\n    geom_point()+\n    geom_line()+\n    geom_text(label=round(loss_vec,2),vjust=1,hjust=-1)+\n    theme(text = element_text(size=25))","metadata":{"execution":{"iopub.status.busy":"2022-07-14T00:44:20.490000Z","iopub.execute_input":"2022-07-14T00:44:20.493058Z","iopub.status.idle":"2022-07-14T00:45:05.995232Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# \n#\n#\n# ### PLACEHOLDER  --- KICKING OUT TEST PREDICTIONS BELOW ### #\n#\n#\n#","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# pick a model\ntest_model_coefs <- coef_list[[1]]","metadata":{"execution":{"iopub.status.busy":"2022-07-14T01:07:36.694435Z","iopub.execute_input":"2022-07-14T01:07:36.697566Z","iopub.status.idle":"2022-07-14T01:07:36.719407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# fit the model on all the data\ntest_model = glm(target ~ ., data=(new_df %>% select(target, test_model_coefs)), family=\"binomial\") \nsummary(test_model)","metadata":{"execution":{"iopub.status.busy":"2022-07-14T01:07:40.038586Z","iopub.execute_input":"2022-07-14T01:07:40.040751Z","iopub.status.idle":"2022-07-14T01:08:14.599597Z"},"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( new_df, train_df, test_df, predictions)\ngc()","metadata":{"execution":{"iopub.status.busy":"2022-07-14T01:08:51.199790Z","iopub.execute_input":"2022-07-14T01:08:51.201539Z","iopub.status.idle":"2022-07-14T01:08:51.578722Z"},"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 = c(\"customer_ID\",\"S_2\",test_model_coefs)) %>% \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)\nncol(test_X)\ncolnames(test_X)\n#test_X","metadata":{"execution":{"iopub.status.busy":"2022-07-14T01:08:54.547308Z","iopub.execute_input":"2022-07-14T01:08:54.548915Z","iopub.status.idle":"2022-07-14T01:15:45.555967Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gc()","metadata":{"execution":{"iopub.status.busy":"2022-07-14T01:23:57.045451Z","iopub.execute_input":"2022-07-14T01:23:57.049677Z","iopub.status.idle":"2022-07-14T01:23:57.232790Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#replace all nas with [fill in blank]\n\nnames <- colnames(test_X) \n\nfor (j in 3:ncol(test_X)){ \n    \n    name <- names[j]    \n    #\n    col <- test_X %>% select(., !!name) %>% pull()    \n    impute_val <- median(col, na.rm=T)\n    new_col <- ifelse(is.na(col),impute_val,col)\n    #    \n    test_X <- test_X %>% mutate( !!name := new_col )\n    \n}\n\n#check\nprint('there should be no missings anymore')\nsum(is.na(test_X))","metadata":{"execution":{"iopub.status.busy":"2022-07-14T01:23:31.579302Z","iopub.execute_input":"2022-07-14T01:23:31.581439Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#make predictions on test data\n\n#predict\npredictions = predict(test_model, newdata=test_X, type=\"response\") %>% replace_na(0) #fallback should definitely be a 0\nhead(predictions)\n\n#check distribution\nhist(predictions)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#bind test_IDs and predictions\ntest_df <- test_X %>% mutate(prediction=predictions) %>% select(customer_ID,prediction)\nnrow(test_df)\nhead(test_df)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#write to submission.csv \nfwrite(x=test_df, file=\"submission.csv\")","metadata":{},"execution_count":null,"outputs":[]}]}