# A first approach using Regularization.

## The problem
The data provided contains thousands of different predictors, namely the *place_ids*. This produces an extremely sparse information distribution, making virtually impossible to tackle it using regular multiple regression and finding significant coefficients based upon ordinary least squares in an straightforward manner.

On the other hand I think we may not need to find the "right" answer, but more precisely a valid and meaningful answer for each situation. We may not need to provide an unique and deterministic answer, but a recommendation that is good enough in each situation.

## The strategy
After analyzing some generic aspects of the data, a first approach may consist on defining which are the relationships among distinct places within the grid. Due to the sparsity of the data *Regularization* seems to be an optimal first way to extract significant information from the data set.

The atomicity of the records, having few non-zero values for each time unit, among thousands of predictors (the place_ids), and the importance of each single point above average weights among places, due to the categorical nature of the variable to analyze, suggests the convenience of *Lasso Regularization* over Ridge. However I think using *ElasticNet* and let the model learn from the data is also a good approach.

It may be not possible to directly make variable selection or to assess the exact variable importance by using Lasso penalty, but I think its predicitive accuracy out-of-the-sample may result of great interest as a first step to build a predictive model for this specific issue. There are many sources that describe how to use Lasso and more generally EN to highlight singificative variables and then proceed with other methods.

I won't go into this, but sparse PCA is another technique suitable for this problem as well, that could help to fit each original record in a clustered space area.

```{r}
library(data.table)
#Get Data
pathToFiles <- '../input/'
numberOfRowsToRead <- -1
set.seed(55446)
dt.train <- read.table(paste0(pathToFiles, "train.csv"),  header=T, sep=",",nrows=numberOfRowsToRead)
```

## Transforming the data set

First, we don't need to study each milimeter of the grid and each second of activity. Let's group the data into 1Kmx1Km grid and expand time 120 times.

```{r}
#Reducing sample size and preparing data for analysis:
dt.train<-as.data.table(dt.train)
dt.train <- dt.train[,c("x","y","time","place_id","accuracy"),with=F]

#X and Y will be removed and moved to an unique XY column called "coordinates".
dt.train <- dt.train[,list(
   x = ifelse(x<10,paste0("0",round(x)),as.character(round(x)))
  ,y = ifelse(y<10,paste0("0",round(y)),as.character(round(y)))
  ,time = factor(round(time/120,0))
  ,place_id = factor(place_id)
  ,accuracy=accuracy
)]
dt.train[,coordinates:=as.factor(paste0(x,y))]
dt.train[,coordinates:=as.factor(coordinates)]

dt.train<-dt.train[,-c("x","y"),with=F]
```

Now the data is ready and we can focus on one of the cells in the grid. We will just select that cell in the grid with more events in this example and study only it.
```{r}
#Getting number of transactions per coordinate, place_id and time.
dt.train.summary <- dt.train[,list(transactionsNumber=length(accuracy)),by=list(coordinates,place_id,time)]

#From the previous dataset, now select which coordinate has more records than the rest.
dt.train.top.coordinates <-  dt.train.summary[,c("coordinates","transactionsNumber"),with=F]

dt.train.top.coordinates<-dt.train.top.coordinates[,list(transactionsNumber=sum(transactionsNumber)),by=coordinates]

dt.train.top.coordinates<-dt.train.top.coordinates[order(-transactionsNumber)]

#Selecting the coordinate with more transactions.
topCoordinate <- as.vector(dt.train.top.coordinates$coordinates[1])

print(topCoordinate)

#Extracting the subset with more transactions.
dt.train.summary.subset <- dt.train.summary[coordinates==topCoordinate]

```

We now want to understand how each place behaves at a certain ploint. There are many choices present at each location. How does each of these places contribute to the total number of events in absolute terms? That is the answer we want to find in this example.

In order to respond we are going to compare the total number of transactions with the total number of *place_ids* selected on each row of data. All the *place_ids* are transformed into column names and for each time unit we put the number of transactions corresponding to each.

Then, we add a new column that reflects the total number of places selected on that row to get an idea of which places are selected more than the rest. This new column is called *totalOperations*.
```{r}
#Converting place_ids in columns
library(reshape2)
dt.train.summary.subset.wide <- as.data.table(dcast(dt.train.summary.subset
                                                  ,coordinates + time  ~ place_id
                                                  ,fill=0))

#Getting number of operations per time unit 
dt.train.summary.subset.wide.operations <- dt.train.summary.subset.wide[,-c(1:2),with=F] 

dt.train.summary.subset.wide.operations[,totalOperations:=                                   apply(dt.train.summary.subset.wide.operations,1
               ,FUN = function(x) length(which(x>0)))
                                 ]
dt.train.summary.subset.wide[,totalOperations:=dt.train.summary.subset.wide.operations$totalOperations]
```
Now we have all the data set ready for analysis. As mentioned above Lasso regression seems a priori the best initial choice due to the atomicity of the data, but we will use ElasticNet and let the model to decide.
```{r}
#Getting test and train X and Y.
trainProportion <- 0.75
smp_size <- floor(trainProportion * nrow(dt.train.summary.subset.wide))

## set the seed to make your partition reproductible
train_ind <- sample(seq_len(nrow(dt.train.summary.subset.wide)), size = smp_size)

Y <- dt.train.summary.subset.wide$totalOperations
X <- dt.train.summary.subset.wide[,-c("totalOperations","coordinates","time"),with=F]

#Train and test subsets.
y_train <- Y[train_ind]
y_test <- Y[-train_ind]
x_train <- X[train_ind]
x_test <- X[-train_ind]
```

Now we work with the ElasticNet model.
```{r}
library(glmnet)

#ElasticNet Regularization.
mod.en.cv <- cv.glmnet(as.matrix(x_train),y_train,alpha=0.5)

#Mean Squared Error plot
plot(mod.en.cv)

#Use alpha with less squared error to predict:
en.pred=predict(mod.en.cv,as.matrix(x_test),s = "lambda.min")

en.mse=mean((en.pred-y_test)^2)
sprintf("MSE_ols absolute = %.3f",en.mse)  
en.mse = sqrt(en.mse) /median(y_test)
sprintf("MSE_ols relative = %.3f",en.mse)  

```

And finally the coefficients. We may be interested on those non-zero.

```{r}
#Got this useful function from stackOverflow to extract coefficients.
get_glmnet_coefs <- function(cvfit, s="lambda.min") {
  ind <- which(coef(cvfit, s=s) != 0)
  df <- data.table(
    feature=rownames(coef(cvfit, s=s))[ind],
    coeficient=coef(cvfit, s=s)[ind]
  )
  #kable(df)
  df <- df[order(-coeficient)]
  df
}

#A list with coefficients for each predictor (This just says that nonzero ones should be taken into account):
mod.en.coefficients <- get_glmnet_coefs(mod.en.cv) 

#Show 50 initial and 50 final
rbind(
  head(mod.en.coefficients,50)
  ,tail(mod.en.coefficients,50)
)
```
```