{"cells":[{"metadata":{"_uuid":"4e5a463b4a259b096e605bcfc6fcac870714ccfb"},"cell_type":"markdown","source":"**Introduction**\n\nIn R, converting the session based data to visitor based aggregation is no problem, but it takes forever without special tricks. This kernel mainly demonstrates how that can be done (see part 3).\n\nThis kernel follows the line of approach of Olivier's https://www.kaggle.com/ogrellier/teach-lightgbm-to-sum-predictions .  It is not a 1-1 translation from Python to R but rather adapts mostly years old code to this situation. \n\nFinally it also shows some very nasty pitfalls that can be caused by using tidyverse and data.table.\n\n\n**Part 0: Libraries and other preliminaries**\n"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"## Importing packages silently\nsuppressPackageStartupMessages({\n    # Data Manipulation: on data. -frames ( or -tables) because standard R is too slow\n    library(tidyverse) # I don't like it, but...\n    library(data.table) # I don't like it either, but...\n    library(splitstackshape) # obscure packge to reshape nested frames\n    library(lubridate) # Because standard date manipulation is never good enough\n\n    # JSON from Jeroen Ooms\n    library(jsonlite)\n\n    # ML from everywhere\n    library(xgboost)\n    library(lightgbm)\n    library(catboost)\n})\n## -----------------------------------------------------\n## options - to prevent some unwanted side effects\noptions(stringsAsFactors= FALSE)\noptions(digits= 4)\n\nuntibble <- function(tibble) {\n  if(\"tbl\" %in%  class(tibble)) {\n    return(as.data.frame(unclass(tibble)))\n  } else return(tibble)\n}\n\nfast_csv <- function(...){\n    require(readr)\n    tib <- suppressMessages(read_csv(...))\n    return(untibble(tib))\n}\n\nrmse <- Metrics::rmse\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5b45e008cc1511a3bb11fe758da8e647df7cbcd4"},"cell_type":"markdown","source":"**Part 1: make the data ready for modeling**\n\nMake a dataset. Output is train, test and id that maps the original  visitor key to a unique number."},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"# make a data-set like Olivier's \n\nn.time <- Sys.time()\ncat(\"\\n### Start...\")\n\ntrain <- fast_csv(\"../input/train.csv\")\ntest <- fast_csv(\"../input/test.csv\")\nboth <- rbind(train, test)\nboth <- cbind(id=1: nrow(both), both)\nl.train <- nrow(train)\n\nb_device <- paste(\"[\", paste(both$device, collapse = \",\"), \"]\") %>% fromJSON(flatten = TRUE)\nb_geoNetwork <- paste(\"[\", paste(both$geoNetwork, collapse = \",\"), \"]\") %>% fromJSON(flatten = TRUE)\nb_totals <- paste(\"[\", paste(both$totals, collapse = \",\"), \"]\") %>% fromJSON(flatten = TRUE)\nb_trafficSource <- paste(\"[\", paste(both$trafficSource, collapse = \",\"), \"]\") %>% fromJSON(flatten = TRUE)\nboth <- both %>%\n  cbind(b_device, b_geoNetwork, b_totals, b_trafficSource) %>%\n  select(-device, -geoNetwork, -totals, -trafficSource)\nrm(train); rm(test); rm(b_device); rm(b_geoNetwork); rm(b_totals); rm(b_trafficSource)\ncat(\"\\n### JSON conversion in...\", Sys.time() - n.time)\n# --------------------------------------\n# now use data table functions\n\nsetDT(both)\nboth[, campaignCode := NULL]\n# delete some columns\nboth[is.na(both$adwordsClickInfo.isVideoAd), \n     adwordsClickInfo.isVideoAd := FALSE]\nfeats<- copy(names(both)) ## very nasty:  for (nm in names(both)) does not work\nfor (nm in feats){\n  v <-length(unique(both[[nm]]))\n  if (v == 1) {\n    both[, (nm):= NULL]\n    ## cat(\"\\n### Deleted column ...\", nm)\n  } ## else cat(\"\\n### Skip column ...\", nm, v)\n}\n#-----------------------------------------------------------\n# Columns that Olivier ddidnt use:\nfeats <- c(\"adwordsClickInfo.page\", \"adwordsClickInfo.slot\"\n           , \"adwordsClickInfo.gclId\", \"adwordsClickInfo.adNetworkType\")\nfor (nm in feats) both[, (nm):= NULL]\ncat(\"\\n### Dimensions of data ...\", dim(both)) # should be 31\n# time feats\nv <- as.POSIXlt(both[[\"visitStartTime\" ]], origin=\"1970-01-01 00:00.00 UTC\")\nboth[, v_dow :=  as.integer(wday(v))]\nboth[, v_hours := as.integer(hours(v))]\nboth[, v_dom := as.integer(day(v))]\nexcluded_features = c(\n    'date', 'fullVisitorId', 'sessionId', 'transactionRevenue', \n    'visitId', 'visitStartTime', 'id'\n)\n\ncategory_feats <- setdiff(names(both)[sapply(names(both), function(nm){\n    class(both[[nm]]) == \"character\"})],excluded_features) \n# don't clean_up yet, just factorize\nfor(f in category_feats) both[, (f) := as.integer(as.factor(both[[f]]))]\n\n# ------------------------------------------------------------\nu <- unique(both$fullVisitorId); l<- length(u)\nid <- data.frame(id=l, Vid=u)\nboth[, fullVisitorId := match(fullVisitorId, u)]\nsetDF(both)\ntrain <- both[1:l.train,]\ntest  <- both[-(1:l.train), ]\ntrain$transactionRevenue <- ifelse(is.na(train$transactionRevenue), \n                                0, as.numeric(train$transactionRevenue))\n#### train, test and id to next stage","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6436531dd758fff8168fd4831312655ae4f54282"},"cell_type":"markdown","source":"**Part 2: stage one model on session data**\n\nMake a model at session level. Okay, lets try a lightgbm on 5 subsets"},{"metadata":{"trusted":true,"_uuid":"3be74521f01267b484f51e82aaee4d0e9fe24b03"},"cell_type":"code","source":"# prepare data\nforbidden = c(\n    'date', 'fullVisitorId', 'sessionId', 'transactionRevenue', \n    'visitId', 'visitStartTime', 'id'\n)\ny<- log1p(train$transactionRevenue)\nallowed<- setdiff(names(train), forbidden)\nd_train<- train[, allowed]\nd_test <- test[, allowed]\n# ---------------------------------------------------------\n# make model\n\nNROUNDS=5\nfolds <- 1+ (1:nrow(train) %%NROUNDS) # stratify on Vid is better\ndev.result <- rep(0, nrow(train))\ntst.result <- matrix(0, nrow=NROUNDS, ncol=nrow(test))\nparam= list( ## far from optimal, but quick\n  boosting = \"gbdt\",\n  num_leaves = 31, \n  learning_rate = 0.03,\n  subsample= 0.9,\n  colsample_bytree= 0.9\n)\ncat( \"\\n=== Params\", unlist(param), \"===\\n\")\nfor (this.round in 1:NROUNDS){\n      cat(\"\\n=== Making Model\", this.round, \":\",  format(Sys.time()), \"===\\n\")      \n      valid <- c(1:length(y)) [folds == this.round]\n      dev <- c(1:length(y)) [folds != this.round]\n\n      dtrain<- lgb.Dataset(data= as.matrix(d_train[dev,]), \n                           label= y[dev])\n      dvalid <- lgb.Dataset(data= as.matrix(d_train[valid,]) , \n                            label= y[valid])\n      valids <- list(val = dvalid)\n\n      Starttime <- Sys.time()\n      cat(\"\\n=== Start run\", format(Starttime), \"===\\n\")\n      model<- lgb.train(data = dtrain,\n                params= param, \n                nrounds=10000,\n                valids = valids,\n                metric=\"rmse\",\n                obj = \"regression\",\n                eval_freq = 250, \n                early_stopping_rounds=50,\n                verbose= 1\n      )\n      pred <- predict(model,as.matrix(d_train[valid,]))\n      dev.result[valid] <- pred\n      pred_t  <- predict(model, as.matrix(test[, allowed]))\n      tst.result[this.round, ]<- pred_t\n      if (this.round == 1 ) tst <- tst.result[1,] else\n        tst <- log1p(colMeans(expm1(tst.result[1:this.round,])))\n      tst[tst < 0] <- 0\n}\ncat(\"\\n### Overall RMSE:\", rmse(dev.result, y))\n# ---------------------------------------------------------\n#  model output is dev.result, dev.test[1:NROUNDS, ]\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9ff087b5ed1b6044e9af9e851d763b8da0dad83c"},"cell_type":"markdown","source":"**Part 3: from session data to visitor data**\n\nFrom session level to user level is the most time-critical section of this kernel. In stndard R it is almost impossible and it takes hours. Here I go to an obscure package and tricky DT syntax."},{"metadata":{"trusted":true,"_uuid":"17be4df9f0be5faaf31ca62254b648e26b4e066a"},"cell_type":"code","source":"# first prepare the data again\nforbidden = c(\n    'date', 'sessionId', \n    'visitId', 'visitStartTime', 'id'\n)\nallowed <- setdiff(names(train), forbidden)\nd_train <- train[, c(allowed, 'transactionRevenue')]\nd_test  <- test[, allowed]\n\nst.time<- Sys.time()\n\nl.train=nrow(d_train)\ndev.result <- pmax(0, expm1(dev.result))\ntst <- pmax(0, colMeans(expm1(tst.result)))\n\ndf_tr <- data.frame(fullVisitorId= d_train$fullVisitorId, \n          pred= dev.result)\ndf_te <- data.frame(fullVisitorId= d_test$fullVisitorId, \n          pred= tst)\nrequire(splitstackshape)\ndf_tr <- aggregate(pred ~ fullVisitorId, data=df_tr, \n         FUN= function(r) return(as.list(r)))\ndf_tr_plus <- t( sapply(df_tr[,2], \n       function(li){ li<- expm1(unlist(li))\n       c(mean_pred=log1p(mean(li)), max_pred=log1p(max(li)), \n         median_pred= log1p(median(li)), \n         nr_pred = length(li),\n         nr_nb_pred = sum(li != 0))\n         })\n)\ndf_tr<- cbind(df_tr, df_tr_plus)\ndf_tr<- listCol_w(df_tr, \"pred\", fill= NA) \ndf_te <- aggregate(pred ~ fullVisitorId, data=df_te, \n         FUN= function(r) return(as.list(r)))\ndf_te_plus <- t( sapply(df_te[,2], \n      function(li){ li<- expm1(unlist(li))\n      c(mean_pred=log1p(mean(li)), max_pred=log1p(max(li)), \n        median_pred= log1p(median(li)), \n        nr_pred = length(li),\n        nr_nb_pred = sum(li != 0))\n      })\n)\ndf_te<- cbind(df_te, df_te_plus)\ndf_te<- listCol_w(df_te, \"pred\", fill= NA) \n                   \nsetDT(d_train)\nsetDT(d_test)\nfeats <- setdiff(names(d_train), \n       c(\"fullVisitorId\" , \"transactionRevenue\"))\nd <- d_train[, lapply(.SD, mean), by =fullVisitorId, .SDcols = feats]\ndf_tr <- merge(d, df_tr, by=\"fullVisitorId\", all.x=TRUE)\nfeats <- setdiff(names(d_test), \n       c(\"fullVisitorId\" , \"transactionRevenue\"))\nd <- d_test[, lapply(.SD, mean), by =fullVisitorId, .SDcols = feats]\ndf_te <- merge(d, df_te, by=\"fullVisitorId\", all.x=TRUE)\n\ntarget <- d_train[, lapply(.SD, sum), \n           by =fullVisitorId, .SDcols = \"transactionRevenue\"]\nnames(target) <- c(\"fullVisitorId\", \"y\")\nsetDF(df_tr); setDF(df_te); setDF(target)\n\ncat(\"\\n=== Processing took\", Sys.time()-st.time, \"mins\\n\")\n# ---------------------------------------------------------\n## the above three plus id go to the next Phase\nall_data<- ls(); to_remain= c(\"df_tr\", \"df_te\", \"target\", \"id\")\nrm(list=setdiff(all_data, to_remain))\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0092b89de2f1c386d446a8b1ada0299db3befb6b"},"cell_type":"markdown","source":"**Part 4: predict from visitor data**\n\nThe second model - let's take xgboost"},{"metadata":{"trusted":true,"_uuid":"ea38371a2618823d5e9043dee7e01d5ff1811fe1"},"cell_type":"code","source":"# some local parameters\nversion=1 \nTESTFLAG= FALSE\nNROUNDS = 5\nset.seed(1948-02-25)\nsubmit_file <- function( nround) sprintf(\"mrooijer01_%02d_%s_%d.csv\"\n                    , version, ifelse(TESTFLAG, 'T', 'P'), nround)\n# ---- data ------\nforbidden <- c(\"Vid\", \"transactionRevenue\", \"sessionId\" \n               , \"fullVisitorId\"\n               , paste0(\"pred_fl_\", 100:300)\n)\nallowed <- setdiff(names(df_te), forbidden)\ny<-log1p(target$y)\nd_train <- df_tr[, allowed]\nd_test <- df_te[, allowed]\nNAvar= NA\n\ndev.result <- rep(0, length(y))\ntst.result <- matrix(0, nrow=NROUNDS, ncol=nrow(d_test))\ncvs= 1+((1:length(y)) %% NROUNDS)\n\nparam = list( # sub optimal\n  objective=\"reg:linear\",\n  eval_metric= \"rmse\",\n  max_depth = 9,\n  min_child_weight = 6,\n  subsample = 0.7,\n  colsample_bytree = 1,\n  colsample_bylevel = 0.7,\n  base_score= mean(y), \n  tree_method= \"exact\", \n  eta                 = ifelse(TESTFLAG, 0.05 , 0.02), \n  \"my_n_rounds\"       = ifelse(TESTFLAG, 5000, 20000),\n  \"my_delta\"          = ifelse(TESTFLAG, 50, 50), \n  \"my_print_n\"        = ifelse(TESTFLAG, 10, 250),\n  nthreads = 4, \n  \"my_ncol\" = ncol(d_train)\n) \ncat( \"\\n=== Params\", unlist(param), \"===\\n\")\ndtest <- as.matrix(d_test)\ndtest <- xgb.DMatrix(dtest, missing=  NAvar)\n\nlimit_loop <- ifelse(TESTFLAG, 1, NROUNDS)\nfor (this.round in 1:limit_loop){\n  cat(\"\\n=== Making Model\", this.round, \":\",  format(Sys.time()), \"===\\n\")\n  valid <- c(1:length(y)) [cvs == this.round]\n  dev <- c(1:length(y)) [cvs != this.round]\n  dmodel <- xgb.DMatrix(\n    as.matrix(d_train[dev,]),\n    label = y[dev], \n    missing= NAvar)\n  dvalid <- xgb.DMatrix(\n    as.matrix(d_train[valid,]), \n    label = y[valid], missing= NAvar)\n  \n  model <- xgb.train(\n    data = dmodel,\n    nrounds = param[[\"my_n_rounds\"]], \n    params = param,  \n    maximize= FALSE,\n    watchlist = list( val = dvalid),\n    print_every_n = param[[\"my_print_n\"]],\n    early_stopping_rounds = param[['my_delta']]\n  )\n  cat(\"\\n=== Model best Score:\", model$best_score\n      , \"at step\", model$best_iteration)\n  pred <- predict(model, dvalid, ntreelimit = model$best_iteration)\n  dev.result[valid] <- pred\n  \n  pred_t  <- predict(model, dtest, ntreelimit = model$best_iteration)\n  tst.result[this.round, ]<- pred_t\n  if (this.round == 1 ) tst <- tst.result[1,] else\n    tst <- log1p(colMeans(expm1(tst.result[1:this.round,])))\n}\n\nif (! TESTFLAG) {\n  cat(\"\\n\\n=== RMSE:\", Metrics::rmse(dev.result, y), \"===\\n\")\n}\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9e7eb815f0370b6361dbca7f8b8a359472412718"},"cell_type":"code","source":"# write the prediction and look at importance\nsub <- read_csv(\"../input/sample_submission.csv\")\nsub$fullVisitorId  <- id$Vid[df_te$fullVisitorId]\nsub$PredictedLogRevenue <- tst\nwrite.csv(sub, submit_file(this.round), row.names = F)\n\nimp <-xgb.importance(model, feature_names=colnames(d_train))\nimp$rank<- 1:nrow(imp)\n## plot this \n# -------------------\nxgb.plot.importance(imp[1:20,])","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"R","language":"R","name":"ir"}},"nbformat":4,"nbformat_minor":1}