---
title: "Titanic 0"
author: ''
output:
  html_document:
    toc: yes
  pdf_document: null
  pdf_presentation: null
  slidy_presentation:
    css: slidy.css
    widescreen: yes
  word_document: null
---

# Objectives

 + Learn how to fit a random forest model
 
 + We will hindcast who survived the Titanic disaster
 
 + Understand OOB and variable importance in practice
 

# Surviving the titanic?

### The context

<img width=400 src="https://www.selection.ca/wp-content/uploads/sites/15/2012/04/titanic-bateau-1000x675.jpg">

<img width=400 src="https://ds1.static.rtbf.be/article/image/1240x800/8/5/4/652e3bdbd9109dd5a4ed9efebd1b1ecb-1553536091.jpg">

<img width=400 src="https://media2.ledevoir.com/images_galerie/nwdp_587443_439454/image.jpg">


<img width=400 src="https://a57.foxnews.com/media2.foxnews.com/BrightCove/694940094001/2019/09/24/931/524/694940094001_6089174013001_6089176547001-vs.jpg?ve=1&tl=1">

<img width=400 src="https://www.rcinet.ca/fr/wp-content/uploads/sites/2/2019/04/45116283_2280871128622014_1465511320353767424_n.jpg">


<img width=400 src="https://d6jf304m27oxw.cloudfront.net/titanic-conspiracy-the-ship-that-never-sank/titanic-newspaper-headline.jpg">


<img width=400 src="https://cdn.vox-cdn.com/thumbor/Z3M6rKDTsQO3qjwu_a4v-J5rQKg=/0x10:2000x1510/1200x800/filters:focal(0x10:2000x1510)/cdn.vox-cdn.com/uploads/chorus_image/image/46125356/titanicillustration.0.0.jpg">

<img width=400 src="https://www.ilmetapontino.it/wp-content/uploads/2019/07/Titanic.jpg">


<img width=400 src="https://images-na.ssl-images-amazon.com/images/I/91skX9RJkFL._AC_SL1500_.jpg">

<img width=400 src="https://www.digitalfireflymarketing.com/wp-content/uploads/drupal/Titanic.Everyone.Safe__1.jpeg">


<img width=400 src="https://anti-dialectics.co.uk/Titanic-sinking.jpg">


<img width=400 src="https://www.titanicuniverse.com/wp-content/uploads/2016/02/Titanicdrawing8.jpeg">


### The cARPATHIA


<img width=400 src="https://i.ytimg.com/vi/NkG-z9pVJTg/maxresdefault.jpg">

<img width=400 src="https://farm8.staticflickr.com/7185/7049314345_0d9857d73b_b.jpg">


<img width=400 src="https://i.pinimg.com/originals/49/5a/d1/495ad194e31c856823abada2053d7057.jpg">


<a href="https://nationalpost.com/news/canada/charles-joughin-titanic-anniversary-april-15-drunk">
<img width=400 src="https://nationalpost.com/news/canada/charles-joughin-titanic-anniversary-april-15-drunk"></a>

### The deads




# Analysis


### Load the librairies

```{r}
library(Hmisc)
library("caret")
library("rpart")
library("tree")
library("e1071")
library(ggplot2) # Data visualization
library(readr) # CSV file I/O, e.g. the read_csv function
library(randomForest)
set.seed(1)
```

### Load the data

This is a considerable shortcut.

This specific dataset has been prepared by <a href="">Thilaksha Silva</a>

Please see <a href="./TitanicExploratoryAnalysis.html">THIS LINK</a> for an explanatory analysis of the Titanic data.

```{r}
mytrain = read.csv("titanictrain.csv")
mytest  = read.csv("titanictest.csv")
mytitanic = rbind(mytest,mytrain)
```

### Fit a Random Forest model

```{r}
# Fitting Random Forest Classification to the Training set
set.seed(432)
myclassifier = randomForest(as.factor(Survived) ~ ., data = mytrain, importance=TRUE)
```


### How many trees in the forest?

```{r}
# Choosing the number of trees
plot(myclassifier,main="Prediction errors: OOB (black), Death (red), Survival (green) ")
```
The figure shows how the errors get to a statistical equilibrium for a large enough forest.

The OOB error rate converges to around 17%. 

Prediction error of death are lower than prediction error of survival. 

The default value of 500 trees in the `randomForest` function seemd to be ok here.


### Calculate predictions on test


And also the misclassification error.


```{r}
# Predicting the Validation set results
y_pred = predict(myclassifier, newdata = mytest[,-which(names(mytest)=="Survived")])

# Checking the prediction accuracy
table(mytest$Survived, y_pred) # Confusion matrix
error <- mean(mytest$Survived != y_pred) # Misclassification error
paste('Accuracy',round(1-error,4))
```


y_pred
      0   1
  0 103   7
  1  18  50
[1] "Accuracy 0.8596"

# Analysis of predictors

### Variable importance measures


* for each observation in training or test sets, bagging will give a prediction, and the predictions tend to be better than for individual trees
* one downside of the bagging procedure is that it doesn't give a single tree, so has less interprability than the single tree estimate.  **Bagging improves predictability at the expense of interpretation**
* We can record the total amount that the Gini index is decreased by splits over a particular
   predictor, averaged overall all $B$ trees. This gives an overall estimate of importance of
   each predictor variable.
*  We can do the same using prediction accuracy, measuring the average 
   decrease in prediction error resulting from splits over each predictor.
   
   
Two possible measures of variable importance are returned by randomforest:

 + Measure-1: Decrease in accuracy

Computed from permuting OOB data: For each tree, the prediction error on the OOB portion of the data is recorded (error rate for classification, MSE for regression). 

Then the same is done after permuting each predictor variable. The difference between the two are then averaged over all trees, and normalized by the standard deviation of the differences. If the standard deviation of the differences is equal to 0 for a variable, the division is not done (but the average is almost always equal to 0 in that case).

 + Measure-2: Decrease in node impurity
 
 This measure is the total decrease in node impurities from splitting on the variable, averaged over all trees. For classification, the node impurity is measured by the Gini index. For regression, it is measured by residual sum of squares.


```{r}
importance(myclassifier)
```

(measure 1): MeanDecreaseAccuracy 

(measure 2): MeanDecreaseGini

Plot of variable importance

```{r}
varImpPlot(myclassifier)
```

### Can we look at predictors more directly?


Why is title important?

 + I'd better be a Miss or Mrs

```{r}
table(mytest$Title,mytest$Survived)
```


oh well...

```{r}
table(mytest$Title,mytest$Sex)
```


 + I'd better not be a man (but Titanic is an exception!)

```{r}
table(mytest$Sex,mytest$Survived)
```



```{r}
hist(mytest$Fare)
median(mytest$Fare)
```


 + I'd better had money!

```{r}
table(mytest$Fare < 15,mytest$Survived)
```



### How about we change the validation set?

```{r}
rm(mytrain,mytest)
#mytrain
```


```{r}
dotitan <- function(ntree,mtry,iopt){
  
set.seed(432)
acc = NULL
for(i in 1:100){
rm(mytrain,mytest)
nrec=nrow(mytitanic);
ntrain=ifelse(nrec%%2==0,nrec/2,(nrec+1)/2);
#ntrain
train=sample(1:nrec,ntrain,replace=F)
mytrain = mytitanic[train,]
mytest  = mytitanic[-train,]
#train
# Fitting Random Forest Classification to the Training set
#if(iopt == 0){
myclassifier = randomForest(as.factor(Survived) ~ ., data = mytrain, ntree = 500)
#}
#if(iopt == 1){
#myclassifier = randomForest(as.factor(Survived) ~ ., #ntree=ntree, data = mytrain)
#}
#if(iopt == 2){
#myclassifier = randomForest(as.factor(Survived) ~ ., ntree #=ntree, mtry = mtry, data = mytrain)
#}

# Predict response on test covariates
y_pred = predict(myclassifier, newdata = mytest[,-which(names(mytest)=="Survived")])
# Check the prediction accuracy
table(mytest$Survived, y_pred) # Confusion matrix
error <- mean(mytest$Survived != y_pred) # Misclassification error
acc = c(acc,round(1-error,4))

}

return(acc)


}

```



```{r}
h <- dotitan(100,4,0)
hist(h)
mean(h)
```

Now a different exercise. Freeze the validation set


```{r}

iopt = 1
ntree=500
mtry=4
set.seed(432)
acc = NULL
rm(mytrain,mytest)
nrec=nrow(mytitanic);
ntrain=ifelse(nrec%%2==0,nrec/2,(nrec+1)/2);
#ntrain
train=sample(1:nrec,ntrain,replace=F)
mytrain = mytitanic[train,]
mytest  = mytitanic[-train,]

```

Run 50 replicates (for example)

```{r}
for(i in 1:50){
myclassifier = randomForest(as.factor(Survived) ~ ., ntree=ntree, data = mytrain)
# Predict response on test covariates
y_pred = predict(myclassifier, newdata = mytest[,-which(names(mytest)=="Survived")])
# Check the prediction accuracy
table(mytest$Survived, y_pred) # Confusion matrix
error <- mean(mytest$Survived != y_pred) # Misclassification error
acc = c(acc,round(1-error,4))
}
```



```{r}
hist(acc,main="Yap RF is not very variable")
mean(acc)
```


# Credits and links

https://nationalpost.com/news/canada/charles-joughin-titanic-anniversary-april-15-drunk

<ol>
<li><a href="">Credit: data preparation due to: </a>
<li><a href="RFOOB.html">OOB error rate explained</a>
<li><a href="https://nationalpost.com/news/canada/charles-joughin-titanic-anniversary-april-15-drunk">Sinking into alcohol</a>
<li>
</ol>


