---
title: "Digit Recognition"
author: "Clement WALTER"
output: html_document
---

## Introduction
This is my very first Kaggle exercise. I am familiar with `R` and `Rmarkdown` but never really did any _Data Science_ missions.

The goal of this task is to create a metamodel predicting the digit written in a `28 x 28` pixels image with a handwritten digit. Each image is black and white and each pixel has thus an integer value between $0$ and $255$.

## Inputs
Let us load the data.
```{r}
list.files("../input")
train <- read.csv('../input/train.csv')
test <- read.csv('../input/test.csv')
```

The first column of the `train` set is supposed to be the true label of each image:

```{r}
head(names(train))
```

Let us see how many of each sample do we have in hand:

```{r}
prop.table(table(train[,'label']))
```
So it seems that all digits are approximately evenly sampled. This minimises the risk of mispredicting one given digit.

## Exploration
Before going with a prediction model, I am intersted in digging a little bit into conditional statistics. Precisely if there are significant difference between digits regarding the _darkness_ of a picture.

### Picture taint
Since each pixel gets a number between $0$ and $255$, we can first have a look at the average value of the pixels for each digit.
```{r}
train$Darkness <- rowMeans(train[,-1])
train$label <- factor(train$label)
library(ggplot2)
ggplot(train, aes(x=label, y=Darkness)) + geom_boxplot()
```

From this plot it seems that this value is not a significant classifier, apart for the digit $1$ with lots more _white_ pixels. This information may always be useful for final prediction so we add it to `test`.
```{r}
test$Darkness <- rowMeans(test)
```

### Principal Component Analysis
We can also have a look at the distribution of the darkness amongst all the pixels of each digit. In a way, this means doing a Principal Component Analysis.
```{r}
# Remember that the first column is label and that we have added a column Darkness
# Here I do a manual centering because I want black and white to be of equal importance
# Scaling is useless 
pca <- prcomp(train[,2:785] - 255/2, center = FALSE, scale. = FALSE)
```

Let us make some visualisation now. First, we look at the total explained variance against the principal components
```{r}
ss <- summary(pca)
ggplot(data.frame(x=1:784, y = ss$importance[3,]), aes(x,y)) +
  geom_line() +
  geom_vline(xintercept = c(22, 50), linetype = 2)
```

We can see that $90\%$ of the variance is explained by only $20$ components. It seems reasonable to retain only those ones for further exploration. We can also use biplots to look at the repartition of the digits amongst the first components.
```{r}
# # This block does not work on kaggle, try it on your computer
# devtools::install_github("vqv/ggbiplot")
# library(ggbiplot) # install from github
# ggbiplot::ggbiplot(pca, groups = train[,1], var.axes = FALSE)
# # one does not see much
# ggbiplot::ggbiplot(pca, groups = train[,1], var.axes = FALSE, choices = 2:3)
# ggbiplot::ggbiplot(pca, groups = train[,1], var.axes = FALSE, choices = 3:4)
# # looks like a heart !!
```

Even though we can guess that some components will help determine the digits, I don't like that much these pictures. Let us look at the loadings insteads, ie. the distribution of each digits on the principal components.
```{r, fig.show="asis"}
ld <- as.data.frame(as.matrix(train[,2:785])%*%pca$rotation)
ld$label <- train[,'label']
p <- lapply(paste0('PC', 1:9), function(pc){
  eval(parse(text=paste0('ggplot(ld, aes(x=label, y=',pc,')) + geom_boxplot()')))
})
gridExtra::marrangeGrob(p, ncol = 3, nrow = 3)
```

It is now clear that the components will help determine the digits. Let us add to the `train` and `test` data frames the loadings on the $20$ first components.
```{r}
train <- data.frame(train, ld[,1:20])
test <- data.frame(test, as.data.frame(as.matrix(test[,1:784])%*%pca$rotation[,1:20]))
```

## Prediction

### Build the model
We now want to use a `randomForest` for our model.
```{r}
nTrain <- 1e4
nTrees <- 100
sel <- 1:nrow(train) %in% sample(1:dim(train)[1], nTrain, replace = FALSE)
rf <- randomForest::randomForest(label ~ ., data = train[sel,], ntree = nTrees)
```
We can check the importance of each variable:
```{r}
ind <- sort(rf$importance, decreasing = TRUE, index.return = TRUE)
as.matrix(rf$importance[ind$ix,][1:30])
```
Here it appears that our auxilliary variables (Darkness and PCA) play a key role in prediction ! Also let us have a look at the importance of all the variables with a plot:
```{r}
ggplot(data.frame(Gini = c(rf$importance),
                  Variable = 1:805,
                  type = c(rep('Pixel', 784), 'Darkness', rep('PC', 20))),
       aes(y=Gini, x=Variable, color = type)) +
  geom_point()
```

The colors let see what the variables are. It is interesting to notice the behaviour of the importance of the pixels. Let us redo this plot with  `geom_line` instead:
```{r}
ggplot(data.frame(Gini = c(rf$importance),
                  Variable = 1:805,
                  type = c(rep('Pixel', 784), 'Darkness', rep('PC', 20))),
       aes(y=Gini, x=Variable, color = type)) +
  geom_line()
```

Here it is, we can see clear peaks. How many peaks they are? On the plot I can count $25$ of them. But it seams like there is space for $28$ (flat line at the beginning). My guess is that there is one peak per line of pixels. In other words, one peak each 28 pixels. Let us see that:
```{r}
ggplot(data.frame(Gini = c(rf$importance),
                  Variable = 1:805,
                  type = c(rep('Pixel', 784), 'Darkness', rep('PC', 20))),
       aes(y=Gini, x=Variable, color = type)) +
  geom_line() +
  geom_vline(xintercept = 1 + c(0:27)*28, linetype = 2, alpha = 0.5)
```

Bingo! So if we were interested in lightening the model, we could retain _only one pixel per line_. Which one would it be? The peaks seem to be more or less centred (no interesting pixel in the edges...) but this decision is not straight forward to us.

In the end since we do not need this kind of simplification we keep our model with all the pixels. Nevertheless was this study quite surprising to me!

### Predict used train data
```{r}
train_test <- predict(rf, train[!sel,-1])
score <- sum(train_test==train[!sel, 1])/length(train_test)
```
Hence on our training set we are able to reach a score of `r score` good predictions.

### Write output
We now predict the `test` data. We have completed this data frame along the way so it is ready to use.
```{r}
results <- predict(rf, test)
results <- data.frame(ImageId = 1:length(results), Label = results)
write.csv(results, file = "results.csv", row.names = FALSE)
```