GStore: Crafting Reasonable Models *** --- title: "" author: "Alejandro Jiménez Rico" output: html_document: fig_width: 10 fig_height: 7 toc: yes number_sections : yes code_folding: hide theme: cosmo highlight: tango --- In today's competition we a’re challenged to analyze a Google Merchandise Store (also known as GStore, where Google swag is sold) customer dataset to predict revenue per customer. Finally, we have a competition more R-oriented. Our moment to shine. In this notebook, I'll try to explore the given dataset and make some inferences along the way in order to get insights on the way to build a baseline model to get started with. ```{r Libraries, message = FALSE} options(repr.plot.width=7, repr.plot.height=4, fig.align = "center") library(tidyverse) library(data.table) library(jsonlite) library(lubridate) library(magrittr) library(moments) library(caret) library(xgboost) library(dataPreparation) library(mlr) library(stringr) library(rsample) library(rlist) library(harrypotter) ``` # Retrieving Data As always, we have to start by taking a look at the actual data. ```{r} test <- fread("../input/test_v2.csv") test %>% glimpse() ``` ```{r, echo = FALSE, message = FALSE} rm(test);gc() ``` Do you notice those weird patterns in the columns `device`, `geoNetwork`, `totals` and `trafficSource`? This isn't suppose to be a `.csv` file. It has tree-developed features that we need to flatten out. ```{r Flattening, message = FALSE} flatten <- function(x){ pre_flatten <- . %>% str_c(., collapse = ",") %>% str_c("[", ., "]") %>% fromJSON(flatten = T) output <- x %>% bind_cols(pre_flatten(.$device)) %>% bind_cols(pre_flatten(.$geoNetwork)) %>% bind_cols(pre_flatten(.$trafficSource)) %>% bind_cols(pre_flatten(.$totals)) %>% select(-device, -geoNetwork, -trafficSource, -totals) return(output) } train <- read_csv("../input/train_v2.csv") %>% flatten() %>% data.table() test <- read_csv("../input/test_v2.csv") %>% flatten() %>% data.table() ``` # Missing Values: Broad Treatment After flattening the columns, we can spot _a lot_ of hidden missing values. And I mean missing by those `not available in demo dataset` and similars. Make no mistake, those are missing values also and should be treated as that. ```{r Hidden NA, message = FALSE} hidden_na <- function(x) x %in% c("not available in demo dataset", "(not provided)", "(not set)", "", "unknown.unknown", "(none)", "Not Socially Engaged") train <- train %>% mutate_all(funs(ifelse(hidden_na(.), NA, .))) %>% data.table() test <- test %>% mutate_all(funs(ifelse(hidden_na(.), NA, .))) %>% data.table() ``` ```{r Visualisation of NA, echo = FALSE, message = FALSE} train %>% summarise_all(funs(sum(is.na(.))/n()*100)) %>% gather(key = "feature", value = "missing_pct") %>% filter(missing_pct > 0) %>% ggplot(aes(x = reorder(feature,missing_pct), y = missing_pct ) ) + geom_bar(stat = "identity", fill = hp(3)[[3]], colour = "black") + labs(y = "Missing Values (%)", x = "Features") + coord_flip() ``` I don't want to seem careless, but I am just going to remove all columns with 100% missing values right now. They ain't telling us anything and are just occupying memory. ```{r remove 100 NA, message = FALSE} train <- train[,which(unlist(lapply(train, function(x)!all(is.na(x))))),with=F] test <- test[,which(unlist(lapply(test, function(x)!all(is.na(x))))),with=F] ``` # Defining the Target In this competition, our task is - as awlays - to predict a target $T$. The target today is defined as the natural log of the sum of all transactions ($t$) per user ($u$). which can be written down as: $$T_u = ln\left(\sum_{i = 1}^N t_u\right)$$ We can visualize its distribution (from those who have spent some money). ```{r Target, echo = FALSE, message = FALSE} train %>% select(fullVisitorId,transactionRevenue) %>% na.omit() %>% mutate(logRevenue = transactionRevenue %>% as.numeric() %>% log1p()) %>% filter(logRevenue > 0) %>% ggplot(aes(x = logRevenue)) + geom_histogram(aes(y = 100*..count../sum(..count..)), colour = "black", fill = hp(4, house = "Slytherin")[[3]], bins = 50) + xlab("Target") + ylab("(%)") ``` It looks like a beatifully distributed variable to be predicted, but please remember that we are filtering out here all the users that have spent nothing. And assuming that the missing values in the `transactionRevenue` column are simply an absence of transactions. Which is $0$ revenue. ```{r, message = FALSE} train[, transactionRevenue := ifelse(is.na(transactionRevenue), 0, transactionRevenue)] ``` > The Pareto principle - also known as the 80/20 rule - states that, for many events, roughly 80% of the effects come from 20% of the causes. And it is a consequent axiom of business management that _80% of sales come from 20% of clients_. Reality is, most people we have recorded in the dataset don't end up buying things in the store. ```{r Pareto Plot, echo = FALSE, message = FALSE} train %>% select(fullVisitorId,transactionRevenue) %>% group_by(fullVisitorId) %>% summarise(logRevenue = (sum(as.numeric(transactionRevenue)))) %>% mutate(logRevenue = ifelse(logRevenue == 0, "No", "Yes")) %>% ggplot(aes(x = logRevenue, fill = logRevenue)) + geom_bar(aes(y = 100*..count../sum(..count..)), colour = "black" ) + geom_text(aes(label = scales::percent(..count../sum(..count..)), y= 100*..count../sum(..count..) ), stat= "count", vjust = -.5) + theme_minimal() + xlab("Did they spend any money whatsoever?") + ylab("(%)") + scale_fill_hp(discrete = TRUE, house = "Gryffindor", name = "") ``` What these numbers are suggesting is that, before even considering predicting _how much_ money a customer is going to spend, we should begin by thinking _whether_ a customer is going to spend _any money_ whatsoever. This detail is crucial because it forces us to not start by constructing a _regression model_, but a _classification_ one. By now, we should have noticed the situation that we are facing defining the `target` this way. Note that we have far more instances for `target == 0` than for `target == 1`. This is what we call _unbalanced data_. And that might be a problem, because some models that we might build will tend to label every output as `0`, in order to increase their accuracy. And that is actually going to get the model more accuracte, indeed; but is not going to be useful for us. How is that? Imagine that you are building an AI model in order to spot and diagnose illnesses in people. An extremely accurate model would be one that _always_ diagnoses no illness; a model that just claims that every person is healthy. That would be extremely accurate, wouldn't it? After all, most people are healthy the vast majority of their time. Yet being accurate, such a model turns out to be utterly pointless, because is not telling you anything that you didn't know. Anagolously, our model might tend to spot every possible customer as `non-payer` if we just aim for accuracy. Making the model accurate, but pointless. There are a myriad of techniques in order to avoid it, most of involve resampling the data. Our target variable is not that extremely unbalanced, though. I think we can be happy just aiming for another metric different from accuracy, such as [AUC](https://en.wikipedia.org/wiki/Receiver_operating_characteristic#Area_under_the_curve). # EDA This EDA (Exploratory Data Analysis) is not intended to be general purpose. Most EDA out there are made without a clear aim, just to get to know the data, know the company or their product or draw business insights out of them. That is pretty useful from a company's perspective and even interesting. It actually is one of the roles of being a Data Scientist that I enjoy the most. But in this kernel I'd like to focus the attention on prediction. We are building a prediction model and our EDA should be focus upon that direction. ## Time Dependence of the Target variable ```{r} train[, date := ymd(date)] test[, date := ymd(date)] ``` ```{r} train_class <- train %>% mutate(target_class = ifelse(as.numeric(transactionRevenue) > 0, 1, 0)) ``` ```{r} train_class %>% group_by(date) %>% summarise(payer_perc = 100*sum(target_class)/n()) %>% ggplot(aes(x = ymd(date), y = payer_perc)) + geom_line(colour = "black", size = 0.5) + geom_smooth(method = "loess", formula = y ~ x, level = 0.9, colour = hp(5, house = "Ravenclaw")[[1]], size = 1.5) + ylab("(%) of Customers that bought something.") + xlab("Date") ``` What insights can we draw from this plot? First and foremost, we can not spot an obvious trend. Yet we can highlight that the usual percentage bounces between $1\%$ and $2\%$, we might get suspicious if the result of our model outputs something away from that. It seems that the store was somehow inefficient at the end of 2016 and beginning of 2017, and it surely suffered major changes that improved conversion and the tendency dramatically. I'd wonder, however, what the hell happened at the end of April 2017. That is just impressive, though it does not seem to keep improving over time. So I wouldn't expect an increasing trend in the % of payers from now on. What I really find interesting are the periodicity of some of the bounces. If you take a closer look at the plot, you can see that it has periodic spikes. I'd say that those ups and downs might have something to do with the day of the week. Let's check it out. ```{r} train_class %>% mutate(date = date %>% wday(label = TRUE, abbr = FALSE)) %>% group_by(date) %>% summarise(payer_perc = 100*sum(target_class)/n()) %>% ggplot(aes(x = date, y = payer_perc)) + geom_col(colour = "black", fill = hp(5, house = "Ravenclaw")[[1]]) + xlab("") + ylab("(%)") ``` As we suspected, something happens during the weekends. My shot is that more "casual" people enter wandering at the store, just for browsing out of curiosity without any intention of buying anything in the first place. Yes, some people might end up buying something, but I expected that even though the gross number of customers increase on weekends, the percentage of them that buy something would drecrease. ### Weekends The first variable that I'd like to explore is the one that we have already spotted studying the time dependence of the target. A variable that tells us whether the customer came to the store during a weekend. ```{r} train_class %>% mutate(weekend = ifelse((date %>% wday(label = TRUE, abbr = FALSE)) %in% c("Saturday", "Sunday"),TRUE,FALSE)) %>% ggplot(aes(x = as.factor(weekend), fill = target_class == 1)) + geom_bar(position = "fill", colour = "black") + scale_fill_hp(discrete = TRUE, name = "Payed?", house = "Gryffindor") + xlab("Weekend?") + scale_x_discrete(limits = c("FALSE", "TRUE")) + ylab("") ``` We can see that our common sense was somehow right, and the proportion of people that _didn't_ spend any money increases during the weekends. ### Distance between connections A bit more exotic than usual, we can focus a bit on the distance between the connections of the customers to the store. What does that mean? You'd fairly think that a user that keeps coming to the store regularly and with high frequency, has more interest in its products, and hence is more likely to end up spending some money in the store. How do we measure that? With a made-up function that creates a `lag` or a `lead` in the time series ```{r} create_time_fea <- function(fun = lag, n = 1){ tr <- train %>% data.table() tr[, date := as_datetime(visitStartTime)] tr[, time_var := (date - fun(date, n)) %>% as.integer()/3600] output <- tr$time_var rm(tr) return(output) } for(i in 1:5){ train_class[str_c("next_sess", i)] <- create_time_fea(lag,i) train_class[str_c("prev_sess", i)] <- create_time_fea(lead,i) } ``` And what is the observable effect of such variables? ```{r} train_class %>% select(prev_sess1, target_class) %>% na.omit() %>% ggplot(aes(x = prev_sess1, fill = as.factor(target_class))) + geom_histogram(aes(y=..count../sum(..count..)), binwidth = 2) + scale_x_continuous(limits = c(-50,50)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Previous Session Distance") + ylab("") + facet_wrap(.~target_class, scales="free_y") ``` ```{r} train_class %>% select(next_sess1, target_class) %>% na.omit() %>% ggplot(aes(x = next_sess1, fill = as.factor(target_class))) + geom_histogram(aes(y=..count../sum(..count..)), binwidth = 2) + scale_x_continuous(limits = c(-50,50)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Next Session Distance") + ylab("") + facet_wrap(.~target_class, scales="free_y") ``` The difference is subtle. Can you notice whether they are the same? According to Kolmogorov's Theorem, if we have two distributions with cumulative distributions functions defined $F_x$ and $F_y$, from which we measure $D$ defined as $$D = max(F_x - F_y)$$, the resulted distribution of $D$ shall follow a Kolmogorov-Smirnov distribution, if - and only if - the two distributions mentioned above are equivalent. Knowing that, we simply measure that using a KS test. ```{r} np <- train_class %>% filter(target_class == 0) %>% .$prev_sess1 p <- train_class %>% filter(target_class == 1) %>% .$prev_sess1 ks.test(p,np) ``` Such a small pvalue is just saying _hell no_, they are not the same distribution. You can repeat this process for the rest of variables that we have created and see that we have similar results. Which is good enough for me. Though there are some strong [criticism](https://community.dur.ac.uk/r.j.coe/teaching/critsig.htm) on statistical testing, which I find very reasonable and I could argue that hypothesis testing is kind of outdated in an era of Big Data and Machine Learning; apart from the _huge_ problem that supposes to [publication bias](https://www.jstor.org/stable/pdf/2982993.pdf). But let's leave it here and move on for now. ## adContent The straightforward variable that most marketing teams rush to look at is the `adContent`. Which ad the user has been shown to? Did it affect the purchase? Are there some ads more efficient than others? ```{r} train_class %>% group_by(adContent, target_class) %>% summarise(value = n()) %>% ggplot(aes(y = value, x = reorder(adContent, -value))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Ad Content") + ylab("") ``` Firstly, the vast majority of sessions didn't show any ad. Or that is what I understand from the `NA` in this variable. If we remove the missing values, we end up with something a bit more clear, but not clear enough. ```{r} train_class %>% select(adContent, target_class) %>% na.omit() %>% group_by(adContent, target_class) %>% summarise(value = n()) %>% ggplot(aes(y = value, x = reorder(adContent, -value))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Ad Content") + ylab("") ``` Actually, does showing an ad have any effect whatsoever? ```{r} train_class %>% mutate(adContent = ifelse(is.na(adContent), "Yes", "No")) %>% select(adContent, target_class) %>% na.omit() %>% group_by(adContent, target_class) %>% summarise(value = n()) %>% ggplot(aes(y = value, x = reorder(adContent, -value))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Saw any ad?") + ylab("") train_class %>% mutate(adContent = ifelse(is.na(adContent), "Yes", "No")) %>% ggplot(aes(x = adContent, fill = as.factor(target_class))) + geom_bar(colour = "black", position = "fill") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Saw any ad?") + ylab("") ``` In general, nope. Showing whatever ad you have to random customers is not going to get any effect at all. That doesn't mean that some of them might have a positive effect whereas others negative. Maybe if we split this by categories we can understand it. But we have far too many categories in order to draw any useful conclusions. Why? Because the quantity for most of them is negligible. I'd rather wrap them up in broader categories in order to get anything a bit more interesting. ### Keyword Ads The first _keyword_ (yikes) that I stopped at it's been `KeyWord`. There are many ads that are keyword-oriented. Are they useful? ```{r} train_class %>% select(adContent, target_class) %>% na.omit() %>% mutate(keyword_ad = ifelse(grepl("KeyWord", adContent), "Yes", "No")) %>% select(keyword_ad, target_class) %>% group_by(keyword_ad, target_class) %>% summarise(value = n()) %>% ggplot(aes(y = value, x = reorder(keyword_ad, -value))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Keyword-oriented Ad?") + ylab("") train_class %>% select(adContent, target_class) %>% na.omit() %>% mutate(keyword_ad = ifelse(grepl("KeyWord", adContent), "Yes", "No")) %>% select(keyword_ad, target_class) %>% ggplot(aes(x = keyword_ad, fill = as.factor(target_class))) + geom_bar(colour = "black", position = "fill") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Keyword-oriented Ad?") + ylab("") ``` If anything, I'd say that showing keyword-oriented ads could be misproductive. Honestly, though, what we are really seeing is that the difference is negligible. I might keep looking, but I don't see right now any other cohort that could be interesting from this feature. I'd simply drop it. ## channelGrouping ```{r} train_class %>% group_by(channelGrouping, target_class) %>% summarise(value = n()) %>% ggplot(aes(y = value, x = reorder(channelGrouping, -value))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Channel Grouping") + ylab("") train_class %>% ggplot(aes(x = channelGrouping, fill = as.factor(target_class))) + geom_bar(colour = "black", position = "fill") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Channel Grouping") + ylab("") ``` The results are quite informative, don't you think? It seems clear that people comming from `Referral` have a greater tendency to end up paying something than the average guy. We could say exactly the opposite for `Social`, since their results are plainly aweful. The rest are confusing or non-informative. Either by a lack of data or because the distribution is evenly splitted. ## Browser ```{r} train_class %>% group_by(browser, target_class) %>% summarise(value = n()) %>% ggplot(aes(y = value, x = reorder(browser, -value))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Browser") + ylab("") train_class %>% ggplot(aes(x = browser, fill = as.factor(target_class))) + geom_bar(colour = "black", position = "fill") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Browser") + ylab("") ``` If I have to be honest, at the beginning I jsut wanted to drop this feature. Almost everyone is using `Chrome`, which is splitted in the middle and tells us nothing. And the rest are so rare that any meaningful difference could be labeled as an outlier. But it is always a good idea to be a bit creative. What happens if we simplify this feature? We could say that from all these browsers we can distinguish `Chrome`, `No-Chrome-but-still-mainstream` and the rest. Let's plot that and see if we can spot something useful. ```{r} train_class %>% mutate(browser = ifelse(browser %in% c("Safari", "Firefox"), "mainstream", ifelse(browser == "Chrome", 'Chrome', "Other"))) %>% group_by(browser, target_class) %>% summarise(value = n()) %>% ggplot(aes(y = value, x = reorder(browser, -value))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Browser") + ylab("") train_class %>% mutate(browser = ifelse(browser %in% c("Safari", "Firefox"), "mainstream", ifelse(browser == "Chrome", 'Chrome', "Other"))) %>% group_by(browser, target_class) %>% ggplot(aes(x = browser, fill = as.factor(target_class))) + geom_bar(colour = "black", position = "fill") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Browser") + ylab("") ``` We should remember that _Chrome_ is - hands down - the most used browser both in pc and smartphones. And that after all, we are talking about a Merchandising Google Store. Chances are, that if our customer is not using _Chrome_, she might be not a big fan of Google. For now, I would just build a feature of type boolean that tracks whether the customer is using `Chrome` or not, and leave it there. ## Source & Medium ### Medium ```{r} train_class %>% group_by(medium, target_class) %>% summarise(value = n()) %>% ggplot(aes(y = value, x = reorder(medium, -value))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Medium") + ylab("") train_class %>% ggplot(aes(x = medium, fill = as.factor(target_class))) + geom_bar(colour = "black", position = "fill") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Medium") + ylab("") ``` We could appreciate some slight diffferents based on the medium through which the custimers arrived to the store, and we can see that there are some mediums more effective than others. If I was working in the _Analytics_ team of the GStore, I would recommend them to keep the `referral` system. It leads a lot of traffic, and impressively effective, compared to the others. The _path_ of this referral is something that we are going to study in more detail later on. ### Source ```{r} train_class %>% group_by(source, target_class) %>% summarise(value = n()) %>% ungroup() %>% top_n(10, value) %>% ggplot(aes(y = value, x = reorder(source, -value))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Source") + ylab("") top_sources <- train_class %>% group_by(source, target_class) %>% summarise(value = n()) %>% ungroup() %>% top_n(10, value) %>% .$source train_class %>% filter(source %in% top_sources) %>% ggplot(aes(x = source, fill = as.factor(target_class))) + geom_bar(colour = "black", position = "fill") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Source") + ylab("") ``` On the other hand, we have the _Source_. Which tells from which place the traffic came from. Note that I limited the plot to the top 10 with more customers, because the list of sources is humongous; and most of them have so small numbers that we couldn't say anything useful from it. See that the source `mail.googleplex.com` is incredibly efficient. It streams a lot of traffic and almost everyone is paying something. Even though we have the data rebalanced, that is impressive. On the contrary, `youtubecom` seems like a poor source to drag customers from. ### Source/Medium Variable In [this](https://www.kaggle.com/erikbruin/google-analytics-eda-lightgbm-screenshots) fantastic Kernel, I saw that the author created a new variable combining both the `Source` and the `Medium`, which could be a significant feature for a model. ```{r} train_class %>% mutate(source_medium = paste(source, medium, sep="/")) %>% group_by(source_medium, target_class) %>% summarise(value = n()) %>% ungroup() %>% top_n(10, value) %>% ggplot(aes(y = value, x = reorder(source_medium, -value))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Source/Medium") + ylab("") top_sourcemediums <- train_class %>% mutate(source_medium = paste(source, medium, sep="/")) %>% group_by(source_medium, target_class) %>% summarise(value = n()) %>% ungroup() %>% top_n(10, value) %>% .$source_medium train_class %>% mutate(source_medium = paste(source, medium, sep="/")) %>% filter(source_medium %in% top_sourcemediums) %>% ggplot(aes(x = source_medium, fill = as.factor(target_class))) + geom_bar(colour = "black", position = "fill") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Source/Medium") + ylab("") ``` However, the results are not much different that those that could be extracted from just the `Source`. Actually, the main important factors that I see here are whether the customer has come from the googleplex mail service, or whether the customer has come from youtube. One is positive and the other is negative, and that is all - the rest might be just for overfitting. So I would rather not adding this new feature yet. ## Source & Referral Another approach that I see in too many models is just leaving the entire feature as categorical and forget about it. You convert all levels of the factor to a unique number and voilà. This solution makes sense sometimes, I don't think it makes it now. ```{r} train_class$source %>% unique() %>% length() train_class$source %>% unique() %>% sample(15) ``` Take an honest look at that and tell me whether you can spot some categorical pattern on it. Make no mistake, you leave these features as that and put them in an all-purpose algorithm and you are forcing your model to train on noise. You are _begging_ for overfitting. It is far too easy to forget that training a Machine Learning models requires a bit of qualitative analysis of the data; not just hocus pocus on a couple of algorithms. I am a fair advocate of black magic, but this time I'd like to go a bit more old-fashioned and build manual features out of this situations. Same thing happens with the `referralPath` variable: ```{r} train_class$referralPath %>% unique() %>% length() train_class$referralPath %>% unique() %>% sample(15) ``` ### Length That does not mean that these variables are useless, we can draw interesting insights out of them and build useful features. The first thing that caught my eye was how variable the length of this values seem to be. Look at them, you have some of them with the length of a simple word, whereas others are as big as whole sentences. Let's take a look at that. ```{r} train_class %>% mutate(source_length = nchar(source)) %>% group_by(source_length, target_class) %>% summarise(value = n()) %>% ggplot(aes(y = value, x = as.numeric(source_length))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Source Length") + ylab("") ``` Weirdly enough, there are some interesting spikes that we could look at. ```{r} top_s_lengths <- train_class %>% mutate(source_length = nchar(source)) %>% group_by(source_length, target_class) %>% summarise(value = n()) %>% arrange(-value) %>% ungroup() %>% top_n(3) %>% .$source_length train_class %>% mutate(source_length = nchar(source)) %>% mutate(source_length_cat = ifelse(source_length %in% top_s_lengths, "Common", "Uncommon")) %>% group_by(source_length_cat, target_class) %>% summarise(value = n()) %>% ggplot(aes(y = value, x = as.factor(source_length_cat))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Source Length Category") + ylab("") train_class %>% mutate(source_length = nchar(source)) %>% mutate(source_length_cat = ifelse(source_length %in% top_s_lengths, "Common", "Uncommon")) %>% ggplot(aes(x = source_length_cat, fill = as.factor(target_class))) + geom_bar(colour = "black", position = "fill") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Source Length Category") + ylab("") ``` It exist a slight difference, indeed. It would be difficult to say whether it is representative enough to consider it a solid feature. Let's take glance at the same measure from `referralPath`. ```{r} train_class %>% mutate(path_length = nchar(referralPath)) %>% group_by(path_length, target_class) %>% summarise(value = n()) %>% ggplot(aes(y = value, x = as.numeric(path_length))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Path Length") + ylab("") ``` ```{r} top_s_lengths <- train_class %>% mutate(path_length = nchar(referralPath)) %>% group_by(path_length, target_class) %>% summarise(value = n()) %>% arrange(-value) %>% ungroup() %>% top_n(3) %>% .$path_length train_class %>% mutate(path_length = nchar(referralPath)) %>% mutate(path_length_cat = ifelse(path_length %in% top_s_lengths, "Common", "Uncommon")) %>% group_by(path_length_cat, target_class) %>% summarise(value = n()) %>% ggplot(aes(y = value, x = as.factor(path_length_cat))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("path Length Category") + ylab("") train_class %>% mutate(path_length = nchar(referralPath)) %>% mutate(path_length_cat = ifelse(path_length %in% top_s_lengths, "Common", "Uncommon")) %>% ggplot(aes(x = path_length_cat, fill = as.factor(target_class))) + geom_bar(colour = "black", position = "fill") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("path Length Category") + ylab("") ``` Again, same. Now I want to try something different and a bit less obvious. Observing closely the nature of these two variables, we could spot something interesting that has to be with the _depth_ of that value. By _depth_ I mean how deep in some site you are before coming to the store. Let's say that you got this value for the path: `/od/Things-To-Do-in-Silicon-Valley/fl/How-To-Visit-the-Googleplex-the-Google-Head-Office-in-Mountain-View.htm"`. This is far more deep than just `/analytics/web/`, right? How do we know that? By the number of `/` present. Same scenario with `source` but with dots `.` instead of `/`. Hence I suggest building a feature that comprises this sense of _depth_: ```{r} train_class %>% mutate(path_depth = str_count(referralPath, '/')) %>% group_by(path_depth, target_class) %>% summarise(value = n()) %>% ggplot(aes(y = value, x = as.numeric(path_depth))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Path Depth") + ylab("") train_class %>% mutate(source_depth = str_count(source, '.')) %>% group_by(source_depth, target_class) %>% summarise(value = n()) %>% ggplot(aes(y = value, x = as.numeric(source_depth))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Source Depth") + ylab("") ``` Instead of digging deeper into catogires within these values, I'd like to spot something against these features. It is obvious that the length and depth of both `source` and `referralPath` is going to correlate. I mean, given our definitions the deeper the value, the longer is going to be. So I thought of creating another new variable called _complexity_ (lacking a better word) that encapsules the proportion between depth and length. ```{r} train_class %>% mutate(path_depth = str_count(referralPath, '\\/')) %>% mutate(path_length = nchar(referralPath)) %>% mutate(path_complexity = path_depth/path_length) %>% ggplot(aes(x = as.numeric(path_complexity))) + geom_histogram(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Path Complexity") + ylab("") ``` ```{r} train_class %>% mutate(source_depth = str_count(source, '\\.')) %>% mutate(source_length = nchar(source)) %>% mutate(source_complexity = source_depth/source_length) %>% ggplot(aes(x = as.numeric(source_complexity))) + geom_histogram(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Source Complexity") + ylab("") ``` See what happens when the value of `path_complexity` is 1. What would that mean? That the `referralPath` is exactly `/`. Is this a direct link or something like that? I'm not sure, but it makes sense that it is something more _direct_ to the store. What could define a more interested customer, right? Let's check: ```{r} train_class %>% mutate(direct_referral = ifelse(referralPath == "/", "Yes", "No")) %>% ggplot(aes(x = as.factor(direct_referral))) + geom_bar(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Path Complexity") + ylab("") train_class %>% mutate(direct_referral = ifelse(referralPath == "/", "Yes", "No")) %>% ggplot(aes(x = as.factor(direct_referral))) + geom_bar(colour = "black", aes(fill = as.factor(target_class)), position = "fill") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Direct Referral") + ylab("") ``` Indeed, this can be used as variable. We got it. Now let's move to a bit more specific cases. Looking closely at the `source` variable, you can see some interesting kind of sources. For example, those who say `share`. There are many of them, and their source looks kind of weird compared to the rest. ```{r} train_class %>% mutate(shared_source = ifelse(grepl("shared", source), "Yes", "No")) %>% ggplot(aes(x = as.factor(shared_source))) + geom_bar(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Shared Source") + ylab("") train_class %>% mutate(shared_source = ifelse(grepl("login", source), "Yes", "No")) %>% ggplot(aes(x = as.factor(shared_source))) + geom_bar(colour = "black", aes(fill = as.factor(target_class)), position = "fill") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Shared Source") + ylab("") ``` There are only a few of them, and it does not look promising. Now I'd like to look at something similar. Some of the sources say `login` in it. Is this an indicator of a different pattern of coming customers? ```{r} train_class %>% mutate(login_source = ifelse(grepl("login", source), "Yes", "No")) %>% ggplot(aes(x = as.factor(login_source))) + geom_bar(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Login Source") + ylab("") train_class %>% mutate(login_source = ifelse(grepl("login", source), "Yes", "No")) %>% ggplot(aes(x = as.factor(login_source))) + geom_bar(colour = "black", aes(fill = as.factor(target_class)), position = "fill") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Login Source") + ylab("") ``` Kind of the same, there are only a few of them and they don't seem to behave any differently. Sadly, it is usual to spend most of the time looking for something that ends up being not that useful. This is science, folks. ## Device ```{r} train_class %>% group_by(deviceCategory, target_class) %>% summarise(value = n()) %>% ggplot(aes(y = value, x = reorder(deviceCategory, -value))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Device") + ylab("") train_class %>% ggplot(aes(x = deviceCategory, fill = as.factor(target_class))) + geom_bar(colour = "black", position = "fill") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Device") + ylab("") ``` As you might have noted by now, I like to simplify. In this feature I see a main driver: Whether a user is using the desktop or not. Mobile and Tablt interactions seem a bit more casual, less oriented to end up buying something. So I would suggest to build another boolean feature that jusks asks whether a user is in their desktop or not. ## Operating System As it has happened before, we have too many options, and most of them have simply not enough instances. So we are going to filter the top most important ones. ```{r} train_class %>% group_by(operatingSystem, target_class) %>% summarise(value = n()) %>% ungroup() %>% top_n(10, value) %>% ggplot(aes(y = value, x = reorder(operatingSystem, -value))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Operating System") + ylab("") top_os <- train_class %>% group_by(operatingSystem, target_class) %>% summarise(value = n()) %>% ungroup() %>% top_n(10, value) %>% .$operatingSystem train_class %>% filter(operatingSystem %in% top_os) %>% ggplot(aes(x = operatingSystem, fill = as.factor(target_class))) + geom_bar(colour = "black", position = "fill") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Operating System") + ylab("") ``` Now be careful because it seems easier than it is. Note that before this, we have stated that desktop users have more odds to end up buying something. And now we can see that users from `Chrome OS`, `Linux` and `Macintosh` also have higher chances. See the pattern? We are just seeing the same as before. If we decided to say that this is another useful variable, we would be fooling ourselves. This is just redundant information. But we can do a bit of magic, and filter in/out the desktop users and see how this changes. ```{r} train_class %>% group_by(operatingSystem, target_class, deviceCategory) %>% summarise(value = n()) %>% ungroup() %>% top_n(10, value) %>% ggplot(aes(y = value, x = reorder(operatingSystem, -value))) + geom_col(colour = "black", aes(fill = as.factor(target_class))) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Operating System") + ylab("") + facet_grid(ifelse(deviceCategory == 'desktop', "Desktop", "Non-Desktop") ~.) top_os <- train_class %>% filter(deviceCategory == "desktop") %>% group_by(operatingSystem, target_class) %>% summarise(value = n()) %>% ungroup() %>% top_n(10, value) %>% .$operatingSystem train_class %>% filter(operatingSystem %in% top_os) %>% ggplot(aes(x = operatingSystem, fill = as.factor(target_class))) + geom_bar(colour = "black", position = "fill") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + scale_fill_hp(name = "Payed?", discrete = TRUE, house = "Gryffindor") + xlab("Operating System") + ylab("") + facet_grid(ifelse(deviceCategory == 'desktop', "Desktop", "Non-Desktop") ~.) ``` There are some interesting insights that we can extract from here. Windows Phone users seem to have outsanding chances of buying something in the store. But don't be fooled by small numbers. They are simply not enough people to make that plot significative enough. People that is interested in actually buying goodies in the GStore aren't casual users, are mostly geeks. We can see that in the astonishing proportion appearing in `Chrome OS` and `Linux`, and quite impressive numbers in the `Macintosh`. The more _casual_ the operating system is - and less developer preferred -, the less proportion of customers end up spending money. The interactions between `Operating System` and whether the device is a `Dekstop` are interesting and informative, we should use them. # Feature Engineering Now we are going to put all this amazing stuff that we have found out about the nature of our data, and build interesting features out of them. In order to do that, we are going to build simple elegant functions, so we'll only have to call them in order to performe the Feature Engineering that we have designed. Firstly, we need a function that transforms categorical data into numeric. Most algorithms are not able to process non-numeric variables, and the ones who do they are doing this behind the scences; so better do it yourself and be sure that it is done the way you want. ```{r Custom Functions, message = FALSE, echo = FALSE} numerise_data <- function(data, numeric_columns){ features <- colnames(data) data[is.na(data)] <- 0 for(f in features) { if ((class(data[[f]])=="factor") || (class(data[[f]])=="character")) { levels <- unique(data[[f]]) %>% sort() data[[f]] <- (factor(data[[f]], levels=levels)) %>% as.numeric() }else{ data[[f]] <- data[[f]] %>% as.numeric() } } data[is.na(data)] <- 0 return(data) } ``` But we aren't encoding all categorical variables simply into unordered numbers, aren't we? We have discussed plenty of possibilities for some features in the previous section that we should address. In order to do that in an ordered manner and with elegant code, I like to define a function simply called `feature_engineering` where we do all the necessary mungings in order to create our features, and call it. ```{r} feature_engineering <- function(tr_te, models = NULL, tri = 0){ tr_te <- data.table(tr_te) # Building Simple Features cat("... Building Simple Features ... \n") tr_te[, weekend := ifelse((date %>% wday(label = TRUE, abbr = FALSE)) %in% c("Saturday", "Sunday"),1,0)] tr_te[, wday := date %>% wday(label = TRUE, abbr = FALSE)] tr_te[, year := date %>% year()] tr_te[, month := date %>% month()] tr_te[, hour := as.numeric(hour(as_datetime(visitStartTime)))] tr_te[, channel_referral := ifelse(channelGrouping == "Referral", 1, 0)] tr_te[, channel_social := ifelse(channelGrouping == "Social", 1, 0)] tr_te[, browser := ifelse(browser %in% c("Safari", "Firefox"), "mainstream", ifelse(browser == "Chrome", 'Chrome', "Other"))] tr_te[, is_chrome_the_browser := ifelse(browser == "Chrome", 1, 0) %>% as.numeric()] tr_te[, source_from_googleplex := ifelse(source == 'mail.googleplex.com', 1, 0) %>% as.numeric()] tr_te[, source_from_youtube := ifelse(source == 'youtube.com', 1, 0) %>% as.numeric()] tr_te[, is_medium_referral := ifelse(medium == 'referral', 1, 0) %>% as.numeric()] tr_te[, is_device_desktop := ifelse(deviceCategory == 'desktop', 1, 0) %>% as.numeric()] tr_te[, is_device_macbook := is_device_desktop*ifelse(operatingSystem == "Macintosh", 1, 0)] tr_te[, windows_desktop := is_device_desktop*ifelse(operatingSystem == 'Windows', 1, 0)] tr_te[, is_device_chromebook := ifelse(operatingSystem == "Chrome OS", 1, 0)] tr_te[, is_device_linux := ifelse(operatingSystem == "Linux", 1, 0)] tr_te[, is_phone_ios := ifelse(operatingSystem == "iOS", 1, 0)] tr_te[, is_phone_android := ifelse(operatingSystem == "Android", 1, 0)] tr_te[, is_phone_windows := ifelse(operatingSystem == "Windows Phone", 1, 0)] tr_te[, single_visit := ifelse(visitNumber == 1,1,0) ] tr_te[, hits_ratio := as.numeric(hits)/as.numeric(pageviews)] tr_te[, domain_site := gsub("^.*\\.","", networkDomain)] tr_te[, adwordsClickInfo.isVideoAd := ifelse(is.na(adwordsClickInfo.isVideoAd), 1, 0)] tr_te[, source_shared := ifelse(grepl("shared", source), 1, 0)] tr_te[, source_platform := ifelse(grepl("login", source), 1, 0)] tr_te[, source_depth := str_count(source, '.')] tr_te[, source_length := nchar(source)] tr_te[, source_complexity := source_depth/source_length] tr_te[, path_depth := str_count(referralPath, '/')] tr_te[, path_length := nchar(referralPath)] tr_te[, path_complexity := path_depth/path_length] # tr_te[,which(unlist(lapply(tr_te, function(x)!all(is.na(x))))),with=F] # Manual Combinations cat("... Manual Combinations ... \n") tr_te[, browser_dev := str_c(browser, "_", deviceCategory)] tr_te[, browser_os := str_c(browser, "_", operatingSystem)] tr_te[, browser_chan := str_c(browser, "_", channelGrouping)] tr_te[, campaign_medium := str_c(campaign, "_", medium)] tr_te[, chan_os := str_c(operatingSystem, "_", channelGrouping)] tr_te[, country_adcontent := str_c(country, "_", adContent)] tr_te[, country_medium := str_c(country, "_", medium)] tr_te[, country_source := str_c(country, "_", source_depth)] tr_te[, dev_chan := str_c(deviceCategory, "_", channelGrouping)] tr_te %<>% as_tibble() # Dummy Variables cat("... Dummy Variables ... \n") small_features <- c('channelGrouping','deviceCategory','adwordsClickInfo.slot','adwordsClickInfo.adNetworkType','medium','continent') dummies <- caret::dummyVars( ~ ., data = tr_te[small_features], fullRank=T) tr_te %<>% cbind(predict(dummies, tr_te)) # Automatic Combinations cat("... Automatic Combinations ... \n") for (i in c("city", "continent", "country", "metro", "domain_site", "region", "subContinent")) for (j in c("browser", "deviceCategory", "operatingSystem", "source_depth")) tr_te[str_c(i, "_", j)] <- str_c(tr_te[[i]], tr_te[[j]], sep = "_") # Time Features cat("... Time Features ... \n") create_time_fea <- function(fun = lag, n = 1){ tr_tmp <- tr_te %>% dplyr::slice(tri) %>% data.table() te_tmp <- tr_te %>% dplyr::slice(-tri) %>% data.table() tr_tmp[, time_var := (date - fun(date, n)) %>% as.integer()/3600] te_tmp[, time_var := (date - fun(date, n)) %>% as.integer()/3600] output <- c(tr_tmp$time_var, te_tmp$time_var) rm(tr_tmp, te_tmp) return(output) } for (i in 1:5) { tr_te[str_c("next_sess", i)] <- create_time_fea(lag, i) tr_te[str_c("prev_sess", i)] <- create_time_fea(lead, i) } tr_te %<>% data.table() all_ids <- tr_te$fullVisitorId tr_te <- tr_te %>% select(-visitId, -sessionId, -fullVisitorId) %>% numerise_data() %>% as_tibble() tr_te %<>% data.table() # Clean cat("... Cleaning up ... \n") tr_te[, date := NULL] tr_te[, visitStartTime := NULL] tr_te[, source := NULL] tr_te[, medium := NULL] tr_te[, browser := NULL] tr_te[, operatingSystem := NULL] tr_te[, deviceCategory := NULL] tr_te[, visitNumber := NULL] tr_te[, networkDomain := NULL] tr_te <- tr_te[,which(unlist(lapply(tr_te, function(x)!all(is.na(x))))),with=F] tr_te[is.na(tr_te)] <- 0 # K Means Features cat("... K-Means Features ... \n") tr_te$kmeans_cluster100 <- kmeans(tr_te, centers = 100) %>% .$cluster %>% as.numeric() tr_te$kmeans_cluster10 <- kmeans(tr_te, centers = 10) %>% .$cluster %>% as.numeric() tr_te$kmeans_cluster2 <- kmeans(tr_te, centers = 2) %>% .$cluster %>% as.numeric() return(tr_te) } # Remove Outliers train <- train %>% filter(country != 'Anguilla') %>% data.table() # Define Targets y <- log1p(as.numeric(train$transactionRevenue)) y[is.na(y)] <- 0 y_class <- ifelse(y>0,1,0) tri <- 1:nrow(train) tr_id <- train$fullVisitorId te_id <- test$fullVisitorId train[, transactionRevenue := NULL] train[, campaignCode := NULL] train %<>% as_tibble() has_many_values <- function(x) n_distinct(x) > 1 # Build tr_te tr_te <- train %>% bind_rows(test) %>% select_if(has_many_values) %>% data.table() # Basic Feature Engineering tr_te <- feature_engineering(tr_te) %>% as_tibble() ``` Most of the stuff that you see in there is usual. Typical automatic combinations between categorical variables, creating dummy variables out of small categorical variables and all the insights that we have extracted from the EDA. How do we check that every crazy variable and combination of them that we have constructed is actually useful? Well, there are many ways to that. One approach that I like to check two simple things: Constant columns and repeated columns. And when I say _repeated_ column, doesn't mean to be literally, it suffices to have two columns with the same shape; in that case one of them is redundant. ```{r} constant_columns <- whichAreConstant(tr_te, verbose=FALSE) redundant_columns <- whichAreBijection(tr_te, verbose=FALSE) print(paste0("Columns that are constant: ", colnames(tr_te)[constant_columns])) print(paste0("Columns that are redundant: ", colnames(tr_te)[redundant_columns])) # if(length(constant_columns > 0)) tr_te <- subset(tr_te,select = -c(constant_columns)) %>% as_tibble() ``` Which means that these columns are not providing _any_ information whatsoever to our model. They are just a waste of memory. So let's get rid of them: ```{r} if(length(constant_columns > 0)) tr_te <- subset(tr_te,select = -c(constant_columns)) %>% as_tibble() if(length(redundant_columns > 0)) tr_te <- subset(tr_te,select = -c(redundant_columns)) %>% as_tibble() ``` # Folding Now we are going to get Folds. What does that mean? Well, in order to prevent [overfitting](https://en.wikipedia.org/wiki/Overfitting) it is common to do what we call [cross validation](https://en.wikipedia.org/wiki/Cross-validation_(statistics)) or _out-of-sample testing_. The thing is to train the model on one data sets and making predictions on another. Trying to predict some data using a model trained on that same data is like asking your grandma who is the most beautiful boy in the world. You'll be told what you want to hear, but deep in your heart you know that the answer is not quite representative. Don't fool yourself, don't overfit. We are going to do this 5 times and average the results. ```{r} get_folds <- function(data, group, v = 5) { group_folds <- group_vfold_cv(data[group], group, v = v) list.zip(tr = tr_ind <- map(group_folds$splits, ~ .x$in_id), val = val_ind <- map(group_folds$splits, ~ setdiff(1:nrow(data), .x$in_id))) } # Get Folds tr_val_ind <- get_folds(train, "fullVisitorId", 5) ``` I'm getting folds now before mentioning any model because now I can do this: ```{r} rm(train,test);gc() ``` And spare some memory before going deeper. # XGB Models Now we are getting to the fancy stuff, Machine Learning! As you'd probably have noticed in the code above, I had saved for later a couple of targets. One continous and one categorical. We have discussed during the introduction of this long-ass kernel that we would better go first for a classification model. But wait, this is a regression problem, we need to minimize `RMSE`. And more problematic, we have instances for every connection of every user that has ever wandered into the GStore, but we need to predict the money spent by user, not by session. How do we commbine all of that? With *iterations*. We are going to train three chained models: * *Classification Model on Sessions*. This model will simply try to predict _whether_ a user has spent any money at each session. Simple, beautiful. * *Regression Model on Sessions*. This guy here is going to get the outcome of the previous model and use it as a feature in order to train another similar model that tries, this time, to predict _how much_ money a user has spent at each session. * *Regression Model on Users*. Finally, our last iteration is going to train a model using the outcome of the previous model as a feature, in which it will try to predict _how much_ money a user has spent _over all his sessions_. Note that from now on I'm going to get a bit technical on XGB models and I assume that you have read (or are able to do it) the documentation and have played with the models a bit. If not, well, good luck. ## Tuning Hyperparameters If you google a bit about XGB, you'll find that one of the trickiest aspects of training them is that you need to choose carefully and very wisely the values for each of the hyperparameters. Some of them lead the model to be prone to overfitting if too high, or too underfitting if too low. So we need to put some values. Of course, there are some standard values that everyone uses, but the rest might change largely from one problem to another. Choosing them is a form of art by itself, but reality is that most people end up choosing them by trial and error. We are going to do that, but more efficiently. We'll craft a function that iterates over random subsamples of the data using different values for those hyperparameters. We are going to evaluate how those hyperparameters effect the performance of our model using Cross Validation. And then we'll chose. ```{r} tune_xgb <- function(train_data, target_label, ntrees = 100, objective = "binary:logistic", eval_metric = "error", fast = TRUE){ train_data <- as.data.frame(train_data) # Count Event Rate if(objective == "binary:logistic") event_rate <- ceiling(1/(sum(train_data$target == 1)/length(train_data$target))) if(!(objective == "binary:logistic")) event_rate <- 10 if(fast){ parameterList <- expand.grid(subsample = seq(from = 1, to = 1, by = 1), colsample_bytree = seq(from = 0.5, to = 1, by = 0.5), lr = seq(from = 2, to = 10, by = 2), mtd = seq(from = 4, to = 10, by = 2), mcw = seq(from = event_rate, to = event_rate, by = event_rate)) }else{ parameterList <- expand.grid(subsample = seq(from = 0.5, to = 1, by = 0.5), colsample_bytree = seq(from = 0.4, to = 1, by = 0.2), lr = seq(from = 1, to = 15, by = 1), mtd = seq(from = 2, to = 16, by = 2), mcw = seq(from = floor(event_rate/2), to = event_rate*10, by = floor(event_rate*2))) } scores <- c() pb <- progress_bar$new( format = " Tuning Hyperparameters [:bar] :percent eta: :eta", total = nrow(parameterList), clear = FALSE, width= 60) for(i in 1:nrow(parameterList)){ pb$tick() # Define Subsample of Training Data sample_size <- floor(nrow(train_data)/100) sample_size <- max(c(sample_size,1e4)) if(nrow(train_data) <= 1e4) sample_size <- nrow(train_data) train_params <- train_data %>% sample_n(sample_size) y_params <- train_params[[target_label]] train_xgb_params <- xgb.DMatrix(data = train_params[,-which(names(train_params) %in% target_label)] %>% as.matrix(), label = y_params) #Extract Parameters to test currentSubSample <- parameterList[["subsample"]][[i]] currentColsampleRate <- parameterList[["colsample_bytree"]][[i]] lr <- parameterList[["lr"]][[i]] mtd <- parameterList[["mtd"]][[i]] mcw <- parameterList[["mcw"]][[i]] p <- list(objective = objective, booster = "gbtree", eval_metric = eval_metric, nthread = 4, eta = lr/ntrees, max_depth = mtd, min_child_weight = mcw, gamma = 0, subsample = currentSubSample, colsample_bytree = currentColsampleRate, colsample_bylevel = 0.632, alpha = 0, lambda = 0, nrounds = ntrees) xgb_cv <- xgb.cv(p, train_xgb_params, p$nrounds, print_every_n = 5, early_stopping_rounds = 25, nfold = 5, verbose = 0) if(eval_metric == "auc") scores[i] <- xgb_cv$evaluation_log$test_auc_mean %>% max() if(eval_metric == "error") scores[i] <- xgb_cv$evaluation_log$test_error_mean %>% min() if(eval_metric == "rmse") scores[i] <- xgb_cv$evaluation_log$test_rmse_mean %>% min() } parameterList$scores <- scores return(parameterList) } ``` See that the function is build such that you can decide which evaluation metric you want. ## First Iteration: Classificadion Model on Sessions ```{r Tuning Hyperparameters} dtest <- xgb.DMatrix(data = data.matrix(tr_te[-tri, ])) pred_tr <- rep(0, length(tri)) pred_te <- 0 train_tune <- tr_te[tri,] train_tune$target <- y_class tuning_scores <- tune_xgb(train_data = train_tune, target_label = "target", ntrees = 100, objective = "binary:logistic", eval_metric = "auc", fast = TRUE) rm(train_tune);gc() tuning_scores %>% melt(id.vars = "scores") %>% ggplot(aes(y = scores, x = as.factor(value))) + geom_boxplot(aes(fill = variable)) + facet_wrap(variable~.) + scale_fill_hp(discrete = TRUE, house = "Ravenclaw") + theme(legend.position ="none") m <- which.max(tuning_scores$scores) currentSubsampleRate <- tuning_scores[["subsample"]][[m]] currentColsampleRate <- tuning_scores[["colsample_bytree"]][[m]] lr <- tuning_scores[["lr"]][[m]] mtd <- tuning_scores[["mtd"]][[m]] mcw <- tuning_scores[["mcw"]][[m]] ``` Once we've got the parameters, we just train the model. ```{r Training First Iteration} ntrees <- 1e3 p <- list(objective = "binary:logistic", booster = "gbtree", eval_metric = "auc", nthread = 4, eta = lr/ntrees, max_depth = mtd, min_child_weight = 30, gamma = 0, subsample = currentSubsampleRate, colsample_bytree = currentColsampleRate, colsample_bylevel = 0.632, alpha = 0, lambda = 0, nrounds = ntrees) for (i in seq_along(tr_val_ind)) { cat("Group fold:", i, "\n") tr_ind <- tr_val_ind[[i]]$tr val_ind <- tr_val_ind[[i]]$val dtrain <- xgb.DMatrix(data = data.matrix(tr_te[tr_ind, ]), label = y_class[tr_ind]) dval <- xgb.DMatrix(data = data.matrix(tr_te[val_ind, ]), label = y_class[val_ind]) cv <- xgb.train(p, dtrain, p$nrounds, list(val = dval), print_every_n = 50, early_stopping_rounds = 300) pred_tr[val_ind] <- (predict(cv, dval, type = "prob")) pred_te <- pred_te + (predict(cv, dtest, type = "prob")) rm(dtrain, dval, tr_ind, val_ind) gc() } pred_tr <- ifelse(pred_tr < 0, 0, pred_tr) pred_te <- ifelse(pred_te < 0, 0, pred_te / length(tr_val_ind)) pred_class <- pred_te xgb.importance(model = cv) %>% as_tibble() %>% top_n(25, Gain) %>% ggplot(aes(x = reorder(Feature, Gain), y = Gain)) + geom_col(aes(fill = log(Gain))) + coord_flip() + scale_fill_hp(house = "Slytherin", name = "Gain") + xlab("Feature") + ylab("Gain") + ggtitle("Features Importance") + theme(legend.position = "none") rm(dtest, cv); gc() ``` So now we just save up this results for later. ```{r} class_pred <- c(pred_tr,pred_te) class_pred_te <- pred_te class_pred_tr <- pred_tr ``` ## Second Iteration: Regression Model on Sessions Now we are getting to the funny stuff. We are building the model that everyone's building. We need to predict a contiuous quantity, we do regression. Is that simple, isn't it? Let's do regression. But wait a minute. For complexity's sake, I thought that could be interesting to divide this iteration in two layers. One layer is a simple model that is just trained used payer's data. Come on. We've built a classification model already, we should be very proficient at knowing which users are going to be payers. We can get far more precise if we build a model that studies those who actually spent something and try to regress that quantity out of them. Then, using the result of our classification model, we can convert that number to $0$ right away for those who we predict that they won't spend anything. So we'll build two models. Payers and everyone. ### Payers Model How do we define the payers training data set? Easy, we just filter those who have paid. ```{r} tr_small <- which(y != 0) y_small <- y[tr_small] train_small <- tr_te %>% dplyr::slice(tri) %>% dplyr::slice(tr_small) ``` Actually is not that easy. Because I spent a lot of time thinking how to divide the training set into validation and actual training with this filtration. What I basically did is creating another index for this intentionally biased training set, and then modify everything that comes after accordingly so we don't screw up the process, either leaving some of the train rows without prediction or by predicting them using a cv model that has been trained with them. It is ugly as hell, so if you think of a smarter way, let me know please. And so we repeat the process above to tune the model, note that this time is a regression model and that we'v added tweaks in the index. ```{r Tuning Second Iteration Small, echo = TRUE, message = FALSE} train_small$target <- y_small tuning_scores <- tune_xgb(train_data = train_small, target_label = "target", ntrees = 100, objective = "reg:linear", eval_metric = "rmse", fast = TRUE) tuning_scores %>% melt(id.vars = "scores") %>% ggplot(aes(y = scores, x = as.factor(value))) + geom_boxplot(aes(fill = variable)) + facet_wrap(variable~.) + scale_fill_hp(discrete = TRUE, house = "Ravenclaw") + theme(legend.position ="none") ``` And also repeat the process of actually training the model. ```{r Training Second Iteration Small} m <- which.min(tuning_scores$scores) currentSubsampleRate <- tuning_scores[["subsample"]][[m]] currentColsampleRate <- tuning_scores[["colsample_bytree"]][[m]] lr <- tuning_scores[["lr"]][[m]] mtd <- tuning_scores[["mtd"]][[m]] mcw <- tuning_scores[["mcw"]][[m]] ntrees <- 1e3 p <- list(objective = "reg:linear", booster = "gbtree", eval_metric = "rmse", nthread = 4, eta = lr/ntrees, max_depth = mtd, min_child_weight = 30, gamma = 0, subsample = currentSubsampleRate, colsample_bytree = currentColsampleRate, colsample_bylevel = 0.632, alpha = 0, lambda = 0, nrounds = ntrees) train_small$fullVisitorId <- tr_id[tr_small] tr_val_ind_small <- get_folds(train_small, "fullVisitorId", 5) dtest <- xgb.DMatrix(data = data.matrix(tr_te[-tri, ])) train_small$fullVisitorId <- NULL train_small$target <- NULL pred_tr <- rep(0, length(tri)) pred_te <- 0 for (i in seq_along(tr_val_ind_small)) { cat("Group fold:", i, "\n") tr_ind <- tr_val_ind_small[[i]]$tr val_ind <- tr_val_ind_small[[i]]$val dtrain <- xgb.DMatrix(data = data.matrix(train_small[tr_ind, ]), label = y_small[tr_ind]) dval <- xgb.DMatrix(data = data.matrix(train_small[val_ind, ]), label = y_small[val_ind]) cv <- xgb.train(p, dtrain, p$nrounds, list(val = dval), print_every_n = 50, early_stopping_rounds = 300) pred_tr[val_ind] <- (predict(cv, dval)) pred_te <- pred_te + (predict(cv, dtest)) rm(dtrain, dval, tr_ind, val_ind); gc() } pred_tr <- ifelse(pred_tr < 0, 0, pred_tr) pred_te <- ifelse(pred_te < 0, 0, pred_te / length(tr_val_ind)) rm(dtest, cv); gc() ``` And with this very simple model and a headache, we've got another feature. ```{r} payers_pred_te <- pred_te*ifelse(class_pred_te >= 0.5,1,0) payers_pred <- c(pred_tr,payers_pred_te) ``` To be continued...