18-titanic.Rmd

# Case Study - Predicting Survivalship on the Titanic {#titanic}

This chapter demonstrates another example of classification with machine learning.  Kaggle made this exercise quite popular.  

In this study, the training and test sets have already been defined, so we

## Import the data.  
We have put our data into our google drive [here](https://drive.google.com/open?id=0ByHtvgo2NGDMV2VBWklMNFpVaVE) and [here](https://drive.google.com/open?id=0ByHtvgo2NGDMaFByZWRxVEJSeDg).  You can find them on Kaggle if need be.    
```{r warning=FALSE, message=FALSE}
library(tidyverse)

train_set <- read_csv("dataset/Kaggle_Titanic_train.csv")
test_set <- read_csv("dataset/Kaggle_Titanic_test.csv")  

## Let's bind both set of data for our exploratory analysis.
df2 <- bind_rows(train_set, test_set)  

## Let's have a first glimpse to our data
glimpse(df2)
```

## Tidy the data  

One can already see that we should put `Survived`, `Sex` and `Embarked` as factor.    
```{r}
df2$Survived <- factor(df2$Survived)
df2$Sex <- factor(df2$Sex)
df2$Embarked <- factor(df2$Embarked)
```

## Understand the data  

This step consists in massaging our variables to see if we can construct new ones or create additional meaning from what we have.  This step require some additional knowledge related to the data and getting familiar with the topics at hand.  

### A. Transform the data  
The great thing about this data set is all the features engineering one can do to increase the predictibilty power of our model.  

#### Dealing with names.  
One of the thing one can notice is the title associated with the name.  The full names on their own might have little predictibility power, but the *title* in the name might have some value and can be used as an additional variables.    
```{r warning=FALSE}
glimpse(df2$Name)  

## gsub is never fun to use.  But we need to strip the cell up to the comma, 
## then everything after the point of the title.
df2$title <- gsub('(.*,)|(\\..*)', "", df2$Name)
table(df2$Sex,df2$title)
```

Some titles are just translations from other languages.  Let's regroup those.  Some other titles aren't occuring often and would not justify to have a category on their own.  We have regroup some titles under common category.  There is some arbitraire in here.    
```{r warning=FALSE}
df2$title <- gsub("Mlle", "Miss", df2$title)
df2$title <- gsub("Mme", "Mrs", df2$title)
df2$title <- gsub("Ms", "Miss", df2$title)
df2$title <- gsub("Jonkheer", "Mr", df2$title)
df2$title <- gsub("Capt|Col|Major", "Army", df2$title)
df2$title <- gsub("Don|Dona|Lady|Sir|the Countess", "Nobility", df2$title)
df2$title <- gsub("Dr|Rev", "Others", df2$title)
df2$title <- factor(df2$title)
df2$title <- factor(df2$title, 
                  levels(df2$title)[c(5, 3, 2, 4, 7, 1, 6)] )
table(df2$Sex, df2$title)
```

It would be also interesting in fact to check the proportion of survivors for each type of title.    
```{r}
round(prop.table(table(df2$Survived, df2$title), 2), 2)
```
We can notice that `Mrs` are more likely to survive than `Miss`.  As expected, our `Mr` have a very low likelyhood of success.  Our `Noble` title managed mostly to survive.  

Our next step is to create a `Last_Name` variable.  This could be helpful as the ways family have escaped the boat might hold some pattens.    
```{r warning=FALSE}
## To get the last name we strip everything after the first comma.    
df2$last_name <- gsub(",.*", "", df2$Name)  

## We can now put this as factor and check how many families.   
df2$last_name <- factor(df2$last_name)
```
So we have `r n_distinct(df2$last_name)` different families on board of the Titanic.  Of course, there might have different families with the same last name. If that's the case, we won't know.   

### A.  Vizualize with families.  
We could add a variable about the family size.    
```{r warning=FALSE}
df2$family_size <- df2$SibSp + df2$Parch + 1
```

If we plot that to check survivalship in function of family size, one can notice interesting patterns.    
```{r titanic_survical_family_pic1, warning=FALSE}
x <- df2[1:891,]
ggplot(x, aes(x = family_size, fill = factor(Survived))) + 
  geom_bar(stat = 'count', position = "dodge") + 
  scale_x_continuous(breaks = c(1:11)) + 
  labs(x = "Family Size", fill = "Survived", 
       title = "Survivalship by Family Size") +
  theme(legend.position = c(0.9, 0.8), panel.background = NULL)
```

Obviously, we only have the survivalship for the train set of data, as we have to guess the test set of data.  So from what we have, there is a clear advantage in being a family of 2, 3 or 4.  We could collapse the variable `Family_Size` into 3 levels.    
```{r warning=FALSE, message=FALSE}
df2$family_size_type[df2$family_size == 1] <- "Singleton"
df2$family_size_type[df2$family_size <= 4 & df2$family_size > 1] <- "Small"
df2$family_size_type[df2$family_size > 4] <- "Large"
df2$family_size_type <- factor(df2$family_size_type, levels = c("Singleton", "Small", "Large"))
```

We can see how many people in each category, then we plot the proportion of survivers in each category.    
```{r titanic_survival_family_pic2, warning=FALSE, message=FALSE}
df3 <- df2[1:891,]
table(df3$Survived, df3$family_size_type)
df3 <- as_tibble(df3)

library(ggmosaic)
ggplot(data = df3) + 
  geom_mosaic(aes(weight = 1, x = product(family_size_type), 
                  fill = factor(Survived), na.rm = TRUE)) + 
  labs(x = "Family Size", y = "Proportion") + 
  theme(panel.background = NULL)
```

Clearly, there is an advantage in being in a family of size 2, 3 or 4; while there is a disadventage in being part of of a bigger family.    

We can try to digg in a bit further with our new family size and titles.  For people who are part of a *Small* family size, which *title* are more likely to surived?    
```{r titanic_survical_title_pic1, warning=FALSE, message=FALSE}
df4 <- df3 %>% dplyr::filter(family_size_type == "Small")
table(df4$Survived, df4$title)

ggplot(data = df4) + 
  geom_mosaic(aes(x = product(title), fill = Survived)) + 
  labs(x = "Survivorship for Small Families in function of their title", 
       y = "Proportion") + 
  theme(panel.background = NULL, axis.text.x = element_text(angle=90, vjust=1))
```

All masters in small families have survived.  Miss & Mrs in small family size have also lots of chane of survival.    

Similarly, for people who embarked alone (*Singleton*), which *title* are more likely to surived?   
```{r titanic_survival_title_pic2, warning=FALSE, message=FALSE}
df4 <- df3 %>% filter(family_size_type == "Singleton")
table(df4$Survived, df4$title)

ggplot(data = df4) + geom_mosaic(aes(x = product(title), fill = Survived)) + 
  labs(x = "Survivorship for people who boarded alone in function of their title", 
       y = "Proportion") + 
  theme(panel.background = NULL, axis.text.x = element_text(angle=90, vjust=1))
```

It might not comes as clear, but we could do the same for title and gender.  Vertically the stacks are ordered as `Singleton` then `Small` then `Large`.    
```{r titanic_survival_title_pic3, message=FALSE, warning=FALSE}
ggplot(data = df3) + geom_mosaic(aes(x = product(family_size_type, title), fill = Survived)) + 
  labs(x = "Survivorship in function of family type and title summary", 
       y = "Proportion") + 
  theme(panel.background = NULL, axis.text.x = element_text(angle=90, vjust=1))
```

## A. Visualize with cabins. 
Although there are many missing data there, we can use the cabin number given to passengers.  The first letter of the cabin number correspond to the deck on the boat.  So let's strip that deck location from the cabin number.   
```{r titanic_survival_class_pic1}
df3$deck <- gsub("([A-Z]+).*", "\\1", df3$Cabin)
df4 <- df3 %>% filter(!is.na(deck))

table(df3$Survived, df3$deck)

ggplot(data = df4) + geom_mosaic(aes(x = product(deck), fill = Survived)) + 
  labs(x = "Survivorship in function of Deck Location", y = "Proportion") + 
  theme(panel.background = NULL, axis.text.x = element_text(angle=90, vjust=1))

detach("package:ggmosaic", unload=TRUE)
```

There is a bit of an anomaly here as it almost as if most people survived.  Now let's keep in mind, that this is only for people which we have their cabin data.    

Let's have a look at how the `Passenger Class` are distributed on the decks. As we are also finishing this first round of feature engineering, let's just mention also how the Passenger Class is affecting survivalship.   
```{r warning=FALSE}
table(df3$Pclass, df3$deck)
round(prop.table(table(df3$Survived, df3$Pclass), 2), 2)
```
More first class people have survived than other classes.  

## B. Transform  Dealing with missing data.  

### Overview.  
I found this very cool package called `visdat` based on `ggplot2` that help us visualize easily missing data.  
```{r titanic_missingdata_pic1, warning=FALSE}
visdat::vis_dat(df2)
```
Straight away one can see that the variables `cabin` and and `Age` have quite a lot of missing data.   
For more accuracy one could check  
```{r warning=FALSE}
fun1 <- function(x){sum(is.na(x))}
map_dbl(df2, fun1)
```

So we can see some missing data in `Fare` and in `Embarked` as well.   
Let's deal with these last 2 variables first.  

#### Basic Replacement.  
We first start with the dessert and the variables that have few missing data.  For those, one can take the median of similar data.  
```{r warning=FALSE}
y <- which(is.na(df2$Embarked))
glimpse(df2[y, ])
```
So the 2 passengers that have no data on the origin of their embarqument are 2 ladies that boarded alone and that shared the same room in first class and that paid $80.  

Let's see who might have paid $80 for a fare.  
```{r titanic_missingdata_pic2, warning=FALSE}
y <- df2 %>% filter(!is.na(Embarked))
ggplot(y, aes(x = Embarked, y = Fare, fill = factor(Pclass))) + 
  geom_boxplot() + 
  scale_y_continuous(labels = scales::dollar, limits = c(0, 250)) + 
  labs(fill = "Passenger \n Class") + 
  geom_hline(aes(yintercept = 80), color = "red", linetype = "dashed", lwd = 1) + 
  theme(legend.position = c(0.9, 0.8), panel.background = NULL)
```

Following this graph, the 2 passengers without origin of embarcation are most likely from "C".  That said, one can argue that the 2 ladies should have embarked from "S" as this is where most people embarked as shown in this table.  
```{r warning=FALSE}
table(df2$Embarked)
```

That said, if we filter our data for the demographics of these 2 ladies, the likelhood of coming from "S" decreased quite a bit. 
```{r}
x <- df2 %>% filter(Sex == "female", Pclass == 1, family_size == 1)
table(x$Embarked)
```

So if we go with median price and with the demographics of the ladies, it would be more likely that they come from "C".  So let's input that.  
```{r warning=FALSE}
df2$Embarked[c(62, 830)] <- "C"
```

Now onto that missing `Fare` data
```{r warning=FALSE}
y <- which(is.na(df2$Fare))
glimpse(df2[y, ])
```
That passenger is a male that boarded in Southampton in third class.  So let's take the median price for similar passagers.  
```{r warning=FALSE}
y <- df2 %>% filter(Embarked == "S" & Pclass == "3" & Sex == "male" & 
                         family_size == 1 & Age > 40)
median(y$Fare, na.rm = TRUE)
df2$Fare[1044] <- median(y$Fare, na.rm = TRUE)
```

#### Predictive modeling replacement.  
First, we'll focus on the `Age` variable.  
There are several methods to input missing data.  We'll try 2 different ones in here.  
But before we can go forward, we have to factorise some variables.  
Let's do the same with `Sibsp` and `Parch` 
```{r warning=FALSE}
df2$Pclass <- factor(df2$Pclass)
```

The first method we'll be using is with the `missForest` package.  
```{r warning=FALSE, message=FALSE, cache=TRUE, results="hide", eval=FALSE}
y <- df2 %>% select(Pclass, Sex, Fare, Embarked, title, family_size, SibSp, Parch, Age)
y <- data.frame(y)

library(missForest)
z1 <- missForest(y, maxiter = 50, ntree = 500)
z1 <- z1[[1]]

# To view the new ages
# View(z1[[1]])

detach("package:missForest", unload=TRUE)
```

The process is fairly rapid on my computer (around 10~15 seconds)  

Our second method takes slightly more time.  
This time we are using the `mice` package.  
```{r warning=FALSE, message=FALSE, cache=TRUE, results="hide", eval=FALSE}
y <- df2 %>% select(Pclass, Sex, Fare, Embarked, title, family_size, SibSp, Parch, Age)
y$Pclass <- factor(y$Pclass)
y$family_size <- factor(y$family_size)
y <- data.frame(y)

library(mice)
mice_mod <- mice(y, method = 'rf')
z2 <- complete(mice_mod)

# To view the new ages
#View(z2[[1]])

detach("package:mice", unload=TRUE)
```

let's compare both type of imputations. 

```{r warning=FALSE, message=FALSE, echo=FALSE}
# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}
```


```{r titanic_missingdata_pic3, warning=FALSE, cache=TRUE}
p1 <- ggplot(df2, aes(x = df2$Age)) + 
        geom_histogram(aes(y = ..density.., fill = ..count..),binwidth = 5) + 
        labs(x = "Age", y = "Frequency", fil = "Survived") + 
        theme(legend.position = "none")
p1
````


```{r eval=FALSE}

p2 <- ggplot(z1, aes(x = z1$Age)) + 
        geom_histogram(aes(y = ..density.., fill = ..count..),binwidth = 5) + 
        labs(x = "Age", y = "Frequency", fil = "Survived") + 
        theme(legend.position = "none")

p3 <- ggplot(z2, aes(x = z2$Age)) + 
        geom_histogram(aes(y = ..density.., fill = ..count..),binwidth = 5) + 
        labs(x = "Age", y = "Frequency", fil = "Survived") + 
        theme(legend.position = "none")

multiplot(p1, p2, p3, cols = 3)
```

It does seem like our second method for imputation follow better our first graph.  So let's use that one and input our predicted age into our main dataframe.  
```{r warning=FALSE}
# df2$Age <- z2$Age
```


### C. Transform More feature engineering with the ages and others.  
Now that we have filled the `NA` for the age variable.  we can massage a bit more that variable.  
We can create 3 more variables: Infant from 0 to 5 years old.  Child from 5 to 15 years old.  Mothers if it is a woman with the variable `Parch` which is greater than one. 
```{r}
df2$infant <- factor(if_else(df2$Age <= 5, 1, 0))
df2$child <- factor(if_else((df2$Age > 5 & df2$Age < 15), 1, 0))
```

```{r}
df2$mother <- factor(if_else((df2$Sex == "female" & df2$Parch != 0), 1, 0))
df2$single <- factor(if_else((df2$SibSp + df2$Parch + 1 == 1), 1, 0))
```


## References.  
* Exploring the titanic dataset from Megan Risdal.  [here](https://www.kaggle.com/mrisdal/titanic/exploring-survival-on-the-titanic)  
* The `visdat` package.  [here](https://github.com/njtierney/visdat)  
* The `ggmosaic` package.  [here](https://github.com/haleyjeppson/ggmosaic)