---
title: "Random Forest And Nearest Neighbors on a Few Blocks"
author: "Alexandru Papiu"
date: "May 18, 2016"
output:
  html_document:
    fig_height: 5
    fig_width: 7
    highlight: tango
    theme: readable
---

In this competition we're given around 30 million (simulated) check-ins on Facebook in a 10km by 10km grid. The goal is to build a model that predicts what business a user checks into based on spatial and temporal information. The tricky part here is that there are around 100k different classes(`place_id`'s) so most supervised learning techniques won't work on the entire dataset. However most classes are clustered in only certain parts of the grid so the idea I'll pursue here is to select a small-ish square within the grid and try to see if we can do better within the small square. First I'll do some exploratory data analysis in the smaller square then I'll use a random forest algorithm for prediction and finally, I'll analyze the results.

### Read and Clean:

Let's load the required packages and read in the data:

```{r, message = FALSE, warning = FALSE}
library(data.table) #reading in the data
library(dplyr) #dataframe manipulation
library(ggplot2) #viz
library(ranger) #the random forest implementation
library(plotly) #3D plotting
library(tidyr) #dataframe manipulation
library(FNN) #k nearest neighbors algorithm
library(xgboost)

fb <- fread("../input/train.csv", integer64 = "character", showProgress = FALSE)
```

Now we'll select a subset of the data - I'll just pick a random 250 meters by 250 meters square in our imaginary Facebook city.

```{r}
fb %>% filter(x >1, x <1.25, y >2.5, y < 2.75) -> fb
head(fb, 3)
```

Notice that time is given to us simply as a numeric value. There have been quite a few great scripts exploring the timestamp and the unit of `time` here is almost certainly minutes. Since businesses tend to run on daily cycles let's extract a new feature called `hour` that gives the hour in the day (from 0 to 24). Let's also extract (approximations) of other time units such as `weekday`, `month` and `year`.

```{r}
fb$hour = (fb$time/60) %% 24
fb$weekday = (fb$time/(60*24)) %% 7
fb$month = (fb$time/(60*24*30)) %% 12 #month-ish
fb$year = fb$time/(60*24*365)
fb$day = fb$time/(60*24) %% 365

```

We will split our dataset into a training and validation set so we can check the results. We choose the validation set to be the more recent check-ins so that our validation structure is similar to the one kaggle does behind the scenes on the test set.

```{r}
small_train = fb[fb$time < 7.3e5,]
small_val = fb[fb$time >= 7.3e5,] 
```


### Exploratory Analysis:

Let's take a look at our training set:
 
```{r, fig.height = 8, fig.width = 10}
ggplot(small_train, aes(x, y )) +
    geom_point(aes(color = place_id)) + 
    theme_minimal() +
    theme(legend.position = "none") +
    ggtitle("Check-ins colored by place_id")
```

Ok, so the clusters are pretty visible, however there seems to be quite a lot of overlap - the place_id's are definitely not separable. Let's try plotting them using the hour component as our third variable. We will just look at the most popular clusters outherwise it gets really messy:

```{r, fig.height = 8, fig.width = 8}
small_train %>% count(place_id) %>% filter(n > 500) -> ids
small_trainz = small_train[small_train$place_id %in% ids$place_id,]

plot_ly(data = small_trainz, x = x , y = y, z = hour, color = place_id,  type = "scatter3d", mode = "markers", marker=list(size= 5)) %>% layout(title = "Place_id's by position and Time of Day")

```

Whoa very cool! Adding the time dimension definitely helps. The daily cycles are clearly visible above - for certain places the check in's stop for a few hours and then start picking up again. Other businesses have quite a few peaks throughtout the day, and the peaks tend to be rather different for different businesses. Also keep in mind that the upper z-square (z = 24) and the lower z-square (z = 0) are really the same thing since time of day is, well, prediodic. So really this thing we're looking at is better viewed not as a cube but as a (flat) solid torus!

Let's also take a look at day of week:
```{r}
plot_ly(data = small_trainz, x = x , y = y, z = weekday, color = place_id,  type = "scatter3d", mode = "markers", marker=list(size= 5)) %>% layout(title = "Place_id's by position and Day of Week")

```

There is some variation by day of week (perhaps some businesses are busier on the weekend) but the most visibile trend is still the day cycles.

However we still might have too many classes for something like random forest to work at its best. Let's check it out:

```{r}
length(unique(small_train$place_id))
```

Let's for now remove the `place_id`'s that have only three or less occurences in the city are we picked. This will decrease the number of classes by a lot. Since we have a validation set we can always come back and change the filter level to see if we get better results.

```{r}
small_train %>% count(place_id) %>% filter(n > 3) -> ids
small_train = small_train[small_train$place_id %in% ids$place_id,]
```

Sweet, now we have **`r nrow(small_train)`** training examples and **`r nrow(ids)`** classes and we're ready to do some machine learning!


### K Nearest Neighbors

The simplest thing to try in this case is probably K Nearest Neighbors. I was a little skeptical about it at first but it works surprisingly well for this problem (I added this section inspired by Sandro's great script [here](https://www.kaggle.com/svpons/facebook-v-predicting-check-ins/grid-knn/code)). 

The tricky part is figuring out the optimal weights for scaling your variables since knn is sensitive to the magnitutde of variables. The weights below are a combination of using a validation set (not the one here) and eye-balling.



```{r}
s = 2
l = 125
w = 500
m = 80

create_matrix = function(train) {
    cbind(s*train$y,
          train$x,
          train$hour/l,
          train$weekday/w,
          train$year/w,
          train$month/w,
          train$time/(w*60*24*7))
    }
    
create_matrix_1 = function(train) {
    cbind(s*train$y,
          train$x,
          train$hour/l)
    }

create_matrix_2 = function(train) {
    cbind(s*train$y,
          train$x)
    }
    
create_matrix_3 = function(train) {
    cbind(s*train$y,
          train$x,
          train$hour/l,
          train$accuracy/m)
    }


X = create_matrix(small_train)
X_val = create_matrix(small_val)

X2 = create_matrix_1(small_train)
X2_val = create_matrix_1(small_val)

X3 = create_matrix_2(small_train)
X3_val = create_matrix_2(small_val)

X4 = create_matrix_3(small_train)
X4_val = create_matrix_3(small_val)

```

Now we're ready to do the knn:

```{r}
model_knn = FNN::knn(train = X, test = X_val, cl = small_train$place_id, k = 15)

preds <- as.character(model_knn)
truth <- as.character(small_val$place_id)
mean(truth == preds)

model_knn2 = FNN::knn(train = X2, test = X2_val, cl = small_train$place_id, k = 15)
preds <- as.character(model_knn2)
truth <- as.character(small_val$place_id)
mean(truth == preds)

model_knn3 = FNN::knn(train = X3, test = X3_val, cl = small_train$place_id, k = 15)
preds <- as.character(model_knn3)
truth <- as.character(small_val$place_id)
mean(truth == preds)

model_knn4 = FNN::knn(train = X4, test = X4_val, cl = small_train$place_id, k = 15)
preds <- as.character(model_knn4)
truth <- as.character(small_val$place_id)
mean(truth == preds)


```

That's a pretty good accuracy I think for a model as simple as knn. Let's see if we can do do any better with a fancier method.