5  ClojisR example: Titanic #0

(ns clojisr.v1.tutorials.titanic0
  (:require [scicloj.kindly.v4.kind :as kind]
            [scicloj.kindly.v4.api :as kindly]))

This notebook is a variation of Pradeep Tripathi’s Titanic Kaggle solution in R. Instead of writing it in R as the original, we write it in Clojure, and call R from Clojure.

The goal is to study the Clojure-R interop, and expecially experiment with various ways to define Clojure functions corresponding to R functions.

In this first, naive version, the corresponding Clojure functions are rather simple. They expect a varying number of arguments, and pass those arguments to R function calls by a rather generic way (as defined by Clojisrs)

We do not try to replace the rather imperative style of the original tutorial. Rather, we try to write something that is as close as possible to the original.

We have leanred a lot from this use case. It did expose lots of issues and open questions about the Clojisr API and implementation. Note, however, that the piece of R code that we are mimicing here is not so typical to the current tidyverse trends – there is no heavy use of dplyr, tidy evaluation, etc. It may be a good idea to study other examples that have more of those.

daslu, Jan. 2020


5.1 Bringing the neecessary R functions

Here are most of the functions that we need, brought by the standard require-r mechanism, inspired by libpython-clj’s require-python (though not as sophisticated at the moment). In function names, dots are changed to hyphens.

(require
 '[clojisr.v1.r :as r
   :refer [r r->clj
           r== r!= r< r> r<= r>= r& r&& r| r||
           str-md
           r+
           bra bra<- brabra brabra<- colon
           require-r]]
 '[clojisr.v1.applications.plotting :refer [plot->svg]]
 '[clojure.string :as string]
 '[clojure.java.io :as io])
(r/set-default-session-type! :rserve)
{:session-type :rserve}
(r/discard-all-sessions)
{}
(require-r
 '[base :refer [round names ! set-seed sum which rnorm lapply sapply %in% table list-files c paste colnames row-names cbind gsub <- $ $<- as-data-frame data-frame nlevels factor expression is-na strsplit as-character summary table]]
 '[stats :refer [median predict]]
 '[ggplot2 :refer [ggsave qplot ggplot aes facet_grid geom_density geom_text geom_histogram geom_bar scale_x_continuous scale_y_continuous labs coord_flip geom_vline geom_hline geom_boxplot]]
 '[ggthemes :refer [theme_few]]
 '[scales :refer [dollar_format]]
 '[graphics :refer [par plot hist dev-off legend]]
 '[dplyr :refer [mutate bind_rows summarise group_by]]
 '[utils :refer [read-csv write-csv head]]
 '[mice :refer [mice complete]]
 '[randomForest :refer [randomForest importance]])
nil

5.2 Introduction – Prediction of Titanic Survival using Random Forest

Pradeep Tripathi’s solution will use randomForest to create a model predicting survival on the Titanic.

5.3 Reading test and train data

This step assumes that the Titanic data lies under the resources/data/ path under your Clojure project.

(def data-path
  (-> (io/resource "data/titanic")
      (io/file)
      (.getAbsolutePath)
      (str "/")))
# Original code:
list.files('../input')
(list-files data-path)
[1] "test.csv.gz"  "train.csv.gz"
# Original code:
train <-read.csv('../input/train.csv', stringsAsFactors = F)
test  <-read.csv('../input/test.csv', stringsAsFactors = F)
(def train (read-csv
            (str data-path "train.csv.gz")
            :stringsAsFactors false))
(def test (read-csv
            (str data-path "test.csv.gz")
            :stringsAsFactors false))

5.4 Combining test and train data

As explained by Thripathi, the Random Forest algorithm will use the Bagging method to create multiple random samples with replacement from the dataset, that will be treated as training data, while the out of bag samples will be treated as test data.

# Original code:
titanic<-bind_rows(train,test)
(def titanic
  (bind_rows train test))

5.5 Data check

# Original code:
str(titanic)
summary(titanic)
head(titanic)
(kind/md
 (str-md titanic))
'data.frame':   1309 obs. of  12 variables:
 $ PassengerId: int  1 2 3 4 5 6 7 8 9 10 ...
 $ Survived   : int  0 1 1 1 0 0 0 0 1 1 ...
 $ Pclass     : int  3 1 3 1 3 3 1 3 3 2 ...
 $ Name       : chr  "Braund, Mr. Owen Harris" "Cumings, Mrs. John Bradley (Florence Briggs Thayer)" "Heikkinen, Miss. Laina" "Futrelle, Mrs. Jacques Heath (Lily May Peel)" ...
 $ Sex        : chr  "male" "female" "female" "female" ...
 $ Age        : num  22 38 26 35 35 NA 54 2 27 14 ...
 $ SibSp      : int  1 1 0 1 0 0 0 3 0 1 ...
 $ Parch      : int  0 0 0 0 0 0 0 1 2 0 ...
 $ Ticket     : chr  "A/5 21171" "PC 17599" "STON/O2. 3101282" "113803" ...
 $ Fare       : num  7.25 71.28 7.92 53.1 8.05 ...
 $ Cabin      : chr  "" "C85" "" "C123" ...
 $ Embarked   : chr  "S" "C" "S" "S" ...
(summary titanic)
  PassengerId      Survived          Pclass          Name          
 Min.   :   1   Min.   :0.0000   Min.   :1.000   Length:1309       
 1st Qu.: 328   1st Qu.:0.0000   1st Qu.:2.000   Class :character  
 Median : 655   Median :0.0000   Median :3.000   Mode  :character  
 Mean   : 655   Mean   :0.3838   Mean   :2.295                     
 3rd Qu.: 982   3rd Qu.:1.0000   3rd Qu.:3.000                     
 Max.   :1309   Max.   :1.0000   Max.   :3.000                     
                NA's   :418                                        
     Sex                 Age            SibSp            Parch      
 Length:1309        Min.   : 0.17   Min.   :0.0000   Min.   :0.000  
 Class :character   1st Qu.:21.00   1st Qu.:0.0000   1st Qu.:0.000  
 Mode  :character   Median :28.00   Median :0.0000   Median :0.000  
                    Mean   :29.88   Mean   :0.4989   Mean   :0.385  
                    3rd Qu.:39.00   3rd Qu.:1.0000   3rd Qu.:0.000  
                    Max.   :80.00   Max.   :8.0000   Max.   :9.000  
                    NA's   :263                                     
    Ticket               Fare            Cabin             Embarked        
 Length:1309        Min.   :  0.000   Length:1309        Length:1309       
 Class :character   1st Qu.:  7.896   Class :character   Class :character  
 Mode  :character   Median : 14.454   Mode  :character   Mode  :character  
                    Mean   : 33.295                                        
                    3rd Qu.: 31.275                                        
                    Max.   :512.329                                        
                    NA's   :1                                              
(head titanic)
  PassengerId Survived Pclass
1           1        0      3
2           2        1      1
3           3        1      3
4           4        1      1
5           5        0      3
6           6        0      3
                                                 Name    Sex Age SibSp Parch
1                             Braund, Mr. Owen Harris   male  22     1     0
2 Cumings, Mrs. John Bradley (Florence Briggs Thayer) female  38     1     0
3                              Heikkinen, Miss. Laina female  26     0     0
4        Futrelle, Mrs. Jacques Heath (Lily May Peel) female  35     1     0
5                            Allen, Mr. William Henry   male  35     0     0
6                                    Moran, Mr. James   male  NA     0     0
            Ticket    Fare Cabin Embarked
1        A/5 21171  7.2500              S
2         PC 17599 71.2833   C85        C
3 STON/O2. 3101282  7.9250              S
4           113803 53.1000  C123        S
5           373450  8.0500              S
6           330877  8.4583              Q

Tripathi: We’ve got a sense of our variables, their class type, 3 and the first few observations of each. We know we’re working with 1309 observations of 12 variables.

5.6 Feature engineering

Thripathi’s explanation: We can break down Passenger name into additional meaningful variables which can feed predictions or be used in the creation of additional new variables. For instance, passenger title is contained within the passenger name variable and we can use surname to represent families.

# Original code:
colnames(titanic)
(colnames titanic)
 [1] "PassengerId" "Survived"    "Pclass"      "Name"        "Sex"        
 [6] "Age"         "SibSp"       "Parch"       "Ticket"      "Fare"       
[11] "Cabin"       "Embarked"   

Retrieve title from passenger names

# Original code:
titanic$title<-gsub('(.*, )|(\..*)', '', titanic$Name)
(def titanic
  ($<- titanic 'title
       (gsub "(.*, )|(\\\\..*)"
             ""
             ($ titanic 'Name))))

Show title counts by sex

# Original code:
table(titanic$Sex, titanic$title)

Clojisr can covert an R frequency table to a Clojure data structure:

(-> (table ($ titanic 'Sex)
           ($ titanic 'title))
    r->clj)

_unnamed [36 3]:

0 1 :$value
female Capt 0
female Col 1
female Don 0
female Dona 4
female Dr 0
female Jonkheer 1
female Lady 1
female Major 0
female Master 1
female Miss 7
male Major 757
male Master 197
male Miss 0
male Mlle 2
male Mme 0
male Mr 0
male Mrs 8
male Ms 0
male Rev 1
male Sir 1
male the Countess 0

Sometimes, it is convenient to first convert it to an R data frame:

(as-data-frame
 (table ($ titanic 'Sex)
        ($ titanic 'title)))
     Var1         Var2 Freq
1  female         Capt    0
2    male         Capt    1
3  female          Col    0
4    male          Col    4
5  female          Don    0
6    male          Don    1
7  female         Dona    1
8    male         Dona    0
9  female           Dr    1
10   male           Dr    7
11 female     Jonkheer    0
12   male     Jonkheer    1
13 female         Lady    1
14   male         Lady    0
15 female        Major    0
16   male        Major    2
17 female       Master    0
18   male       Master   61
19 female         Miss  260
20   male         Miss    0
21 female         Mlle    2
22   male         Mlle    0
23 female          Mme    1
24   male          Mme    0
25 female           Mr    0
26   male           Mr  757
27 female          Mrs  197
28   male          Mrs    0
29 female           Ms    2
30   male           Ms    0
31 female          Rev    0
32   male          Rev    8
33 female          Sir    0
34   male          Sir    1
35 female the Countess    1
36   male the Countess    0

Sometimes, it is convenient to use the way R prints a frequency table.

(table ($ titanic 'Sex)
       ($ titanic 'title))
        
         Capt Col Don Dona  Dr Jonkheer Lady Major Master Miss Mlle Mme  Mr Mrs
  female    0   0   0    1   1        0    1     0      0  260    2   1   0 197
  male      1   4   1    0   7        1    0     2     61    0    0   0 757   0
        
          Ms Rev Sir the Countess
  female   2   0   0            1
  male     0   8   1            0

Convert titles with low count into a new title, and rename/reassign Mlle, Ms and Mme.

# Original code:
unusual_title<-c('Dona', 'Lady', 'the Countess','Capt', 'Col', 'Don',
                 'Dr', 'Major', 'Rev', 'Sir', 'Jonkheer')
(def unusual-title
  ["Dona", "Lady", "the Countess","Capt", "Col", "Don",
   "Dr", "Major", "Rev", "Sir", "Jonkheer"])
# Original code:
titanic$title[titanic$title=='Mlle']<-'Miss'
titanic$title[titanic$title=='Ms']<-'Miss'
titanic$title[titanic$title=='Mme']<-'Mrs'
titanic$title[titanic$title %in% unusual_title]<-'Unusual Title'
(def titanic
  (-> titanic
      (bra<- (r== ($ titanic 'title) "Mlle")
             "title"
             "Miss")
      (bra<- (r== ($ titanic 'title) "Ms")
             "title"
             "Miss")
      (bra<- (r== ($ titanic 'title) "Mme")
             "title"
             "Mrs")
      (bra<- (%in% ($ titanic 'title) unusual-title)
             "title"
             "Mrs")))

Check the title count again:

# Original code:
table(titanic$Sex, titanic$title)

trying again:

(table ($ titanic 'Sex)
       ($ titanic 'title))
        
         Master Miss  Mr Mrs
  female      0  264   0 202
  male       61    0 757  25

Create a variable which contain the surnames of passengers.

# Original code:
titanic$surname<-sapply(titanic$Name, function(x) strsplit(x,split='[,.]')[[1]][1])
nlevels(factor(titanic$surname)) ## 875 unique surnames
(def titanic
  ($<- titanic 'surname
       (sapply ($ titanic 'Name)
               (r '(function [x]
                             (bra (brabra (strsplit x :split "[,.]") 1) 1))))))
(-> titanic
    ($ 'surname)
    factor
    nlevels)
[1] 875

Tripathi: Family size variable: We are going to create a variable “famsize” to know the number of family members. It includes number of sibling/number of parents and children+ passenger themselves

# Original code:
titanic$famsize <- titanic$SibSp + titanic$Parch + 1
(def titanic
  ($<- titanic 'famsize
       (r+ ($ titanic 'SibSp)
           ($ titanic 'Parch)
           1)))

Create a family variable:

# Original code:
titanic$family <- paste(titanic$surname, titanic$famsize, sep='_')
(def titanic
  ($<- titanic 'family
       (paste ($ titanic 'surname)
              ($ titanic 'famsize)
              :sep "_")))

Visualize the relationship between family size & survival:

 ggplot(titanic[1:891,], aes(x = famsize, fill = factor(Survived))) +
   geom_bar(stat='count', position='dodge') +
   scale_x_continuous(breaks=c(1:11)) +
   labs(x = 'Family Size') +
   theme_few()
(-> titanic
    (bra (colon 1 891)
         nil)
    (ggplot (aes :x 'famsize
                 :fill '(factor Survived)))
    (r+ (geom_bar :stat "count"
                  :position "dodge")
        (scale_x_continuous :breaks (colon 1 11))
        (labs :x "Family Size")
        (theme_few))
    plot->svg)
0 100 200 300 1 2 3 4 5 6 7 8 9 10 11 Family Size count factor(Survived) 0 1

Tripathi: Explanation: We can see that there’s a survival penalty to single/alone, and those with family sizes above 4. We can collapse this variable into three levels which will be helpful since there are comparatively fewer large families.

Discretize family size:

# Original code:
titanic$fsizeD[titanic$famsize == 1] <- 'single'
titanic$fsizeD[titanic$famsize < 5 & titanic$famsize> 1] <- 'small'
titanic$fsizeD[titanic$famsize> 4] <- 'large'
(def titanic
  (-> titanic
      (bra<- (r== ($ titanic 'famsize) 1)
             "fsizeD"
             "single")
      (bra<- (r& (r< ($ titanic 'famsize) 5)
                 (r> ($ titanic 'famsize) 1))
             "fsizeD"
             "small")
      (bra<- (r> ($ titanic 'famsize) 4)
             "fsizeD" "large")))

Let us check if it makes sense:

(-> titanic
    ($ 'fsizeD)
    table)

 large single  small 
    82    790    437 

And let us make sure there are no missing values:

(-> titanic
    ($ 'fsizeD)
    is-na
    table)

FALSE 
 1309 

Tripathi: There’s could be some useful information in the passenger cabin variable including about their deck, so Retrieve deck from Cabin variable.

# Original code:
titanic$Cabin[1:28]
(-> titanic
    (bra (colon 1 28)
         "Cabin"))
 [1] ""            "C85"         ""            "C123"        ""           
 [6] ""            "E46"         ""            ""            ""           
[11] "G6"          "C103"        ""            ""            ""           
[16] ""            ""            ""            ""            ""           
[21] ""            "D56"         ""            "A6"          ""           
[26] ""            ""            "C23 C25 C27"

The first character is the deck:

# Original code:
strsplit(titanic$Cabin[2], NULL) [[1]]
(-> titanic
    ($ 'Cabin)
    (bra 2)
    (strsplit nil)
    (brabra 1))
[1] "C" "8" "5"

Deck variable:

# Original R code:
titanic$deck<-factor(sapply(titanic$Cabin, function(x) strsplit(x, NULL)[[1]][1]))
(def titanic
  ($<- titanic 'deck
       (factor (sapply ($ titanic 'Cabin)
                       '(function [x]
                                  (bra (brabra (strsplit x nil) 1) 1))))))

Let us check:

(-> titanic
    ($ 'deck)
    table)

 A  B  C  D  E  F  G  T 
22 65 94 46 41 21  5  1 

5.7 Missing values

updated summary

(summary titanic)
  PassengerId      Survived          Pclass          Name          
 Min.   :   1   Min.   :0.0000   Min.   :1.000   Length:1309       
 1st Qu.: 328   1st Qu.:0.0000   1st Qu.:2.000   Class :character  
 Median : 655   Median :0.0000   Median :3.000   Mode  :character  
 Mean   : 655   Mean   :0.3838   Mean   :2.295                     
 3rd Qu.: 982   3rd Qu.:1.0000   3rd Qu.:3.000                     
 Max.   :1309   Max.   :1.0000   Max.   :3.000                     
                NA's   :418                                        
     Sex                 Age            SibSp            Parch      
 Length:1309        Min.   : 0.17   Min.   :0.0000   Min.   :0.000  
 Class :character   1st Qu.:21.00   1st Qu.:0.0000   1st Qu.:0.000  
 Mode  :character   Median :28.00   Median :0.0000   Median :0.000  
                    Mean   :29.88   Mean   :0.4989   Mean   :0.385  
                    3rd Qu.:39.00   3rd Qu.:1.0000   3rd Qu.:0.000  
                    Max.   :80.00   Max.   :8.0000   Max.   :9.000  
                    NA's   :263                                     
    Ticket               Fare            Cabin             Embarked        
 Length:1309        Min.   :  0.000   Length:1309        Length:1309       
 Class :character   1st Qu.:  7.896   Class :character   Class :character  
 Mode  :character   Median : 14.454   Mode  :character   Mode  :character  
                    Mean   : 33.295                                        
                    3rd Qu.: 31.275                                        
                    Max.   :512.329                                        
                    NA's   :1                                              
    title             surname             famsize          family         
 Length:1309        Length:1309        Min.   : 1.000   Length:1309       
 Class :character   Class :character   1st Qu.: 1.000   Class :character  
 Mode  :character   Mode  :character   Median : 1.000   Mode  :character  
                                       Mean   : 1.884                     
                                       3rd Qu.: 2.000                     
                                       Max.   :11.000                     
                                                                          
    fsizeD               deck     
 Length:1309        C      :  94  
 Class :character   B      :  65  
 Mode  :character   D      :  46  
                    E      :  41  
                    A      :  22  
                    (Other):  27  
                    NA's   :1014  

Thripathi’s explanation, following the summary: - Age : 263 missing values - Fare : 1 missing values - Embarked : 2 missing values - survived:too many - Cabin : too many

Missing value in Embarkment – Tripathi: Now we will explore missing values and rectify it through imputation. There are a number of different ways we could go about doing this. Given the small size of the dataset, we probably should not opt for deleting either entire observations (rows) or variables (columns) containing missing values. We’re left with the option of replacing missing values with sensible values given the distribution of the data, e.g., the mean, median or mode.

To know which passengers have no listed embarkment port:

# Original code:
titanic$Embarked[titanic$Embarked == ""] <- NA
titanic[(which(is.na(titanic$Embarked))), 1]

Marking as missing:

(def titanic
  (bra<- titanic
         (r== ($ titanic 'Embarked) "")
         "Embarked"
         'NA))

Checking which has missing port:

(-> titanic
    (bra (-> titanic
             ($ 'Embarked)
             is-na
             which)
         1))
[1]  62 830

Tripathi: Passengers 62 and 830 are missing Embarkment.

# Original code:
titanic[c(62, 830), 'Embarked']
(-> titanic
    (bra [62 830]
         "Embarked"))
[1] NA NA

Tripathi: So Passenger numbers 62 and 830 are each missing their embarkment ports. Let’s look at their class of ticket and their fare.

# Original code:
titanic[c(62, 830), c(1,3,10)]
(-> titanic
    (bra [62 830]
         [1 3 10]))
    PassengerId Pclass Fare
62           62      1   80
830         830      1   80

Alternatively:

(-> titanic
    (bra [62 830]
         ["PassengerId" "Pclass" "Fare"]))
    PassengerId Pclass Fare
62           62      1   80
830         830      1   80

Thripathi’s explanation: Both passengers had first class tickets that they spent 80 (pounds?) on. Let’s see the embarkment ports of others who bought similar kinds of tickets.

First way of handling missing value in Embarked:

# Original code:
titanic%>%
  group_by(Embarked, Pclass) %>%
  filter(Pclass == "1") %>%
  filter(Pclass == "1") %>%
  filter(Pclass == "1") %>%
  summarise(mfare = median(Fare),n = n())
(-> titanic
    (group_by 'Embarked 'Pclass)
    (r.dplyr/filter '(== Pclass "1"))
    (summarise :mfare '(median Fare)
               :n '(n)))
# A tibble: 4 × 4
# Groups:   Embarked [4]
  Embarked Pclass mfare     n
  <chr>     <int> <dbl> <int>
1 C             1  76.7   141
2 Q             1  90       3
3 S             1  52     177
4 <NA>          1  80       2

Tripathi: Looks like the median price for a first class ticket departing from ‘C’ (Charbourg) was 77 (in comparison to our 80). While first class tickets departing from ‘Q’ were only slightly more expensiive (median price 90), only 3 first class passengers departed from that port. It seems far more likely that passengers 62 and 830 departed with the other 141 first-class passengers from Charbourg.

Second Way of handling missing value in Embarked:

# Original code:
embark_fare <- titanic %>%
  filter(PassengerId != 62 & PassengerId != 830)
embark_fare
(def embark_fare
  (-> titanic
      (r.dplyr/filter '(& (!= PassengerId 62)
                          (!= PassengerId 830)))))

Use ggplot2 to visualize embarkment, passenger class, & median fare:

# Original code:
ggplot(embark_fare, aes(x = Embarked, y = Fare, fill = factor(Pclass))) +
geom_boxplot() +
geom_hline(aes(yintercept=80),
              colour='red', linetype='dashed', lwd=2) +
scale_y_continuous(labels=dollar_format()) +
theme_few()
(-> embark_fare
    (ggplot (aes :x 'Embarked
                 :y 'Fare
                 :fill '(factor Pclass)))
    (r+ (geom_boxplot)
        (geom_hline (aes :yintercept 80)
                    :colour "red"
                    :linetype "dashed"
                    :lwd 2)
        (scale_y_continuous :labels (dollar_format)))
    plot->svg)
$0 $100 $200 $300 $400 $500 C Q S Embarked Fare factor(Pclass) 1 2 3

Tripathi: From plot we can see that The median fare for a first class passenger departing from Charbourg (‘C’) coincides nicely with the $80 paid by our embarkment-deficient passengers. I think we can safely replace the NA values with ‘C’. Since their fare was $80 for 1st class, they most likely embarked from ‘C’.

# Original code:
titanic$Embarked[c(62, 830)] <- 'C'
(def titanic
  (bra<- titanic [62 830] "Embarked"
         "C"))

A missing value in fare. Thripathi’s explanation: To know Which passenger has no fare information:

# Original code:
titanic[(which(is.na(titanic$Fare))) , 1]
(-> titanic
    (bra (-> titanic
             ($ 'Fare)
             is-na
             which)
         1))
[1] 1044

Tripathi: Looks like Passenger number 1044 has no listed Fare

Where did this passenger leave from? What was their class?

# Original code:
 titanic[1044, c(3, 12)]
(-> titanic
    (bra 1044 [3 12]))
     Pclass Embarked
1044      3        S

Tripathi: Another way to know about passenger id 1044 :Show row 1044

# Original code:
titanic[1044, ]
(-> titanic
    (bra 1044 nil))
     PassengerId Survived Pclass               Name  Sex  Age SibSp Parch
1044        1044       NA      3 Storey, Mr. Thomas male 60.5     0     0
     Ticket Fare Cabin Embarked title surname famsize   family fsizeD deck
1044   3701   NA              S    Mr  Storey       1 Storey_1 single <NA>

Thripathi’s explanation: Looks like he left from ‘S’ (Southampton) as a 3rd class passenger. Let’s see what other people of the same class and embarkment port paid for their tickets.

6 First way:

titanic%>% filter(Pclass == ‘3’ & Embarked == ‘S’) %>% summarise(missing_fare = median(Fare, na.rm = TRUE))

(-> titanic
    (r.dplyr/filter '(& (== Pclass "3")
                        (== Embarked "S")))
    (summarise :missing_fare '(median Fare :na.rm true)))
  missing_fare
1         8.05

Tripathi: Looks like the median cost for a 3rd class passenger leaving out of Southampton was 8.05. That seems like a logical value for this passenger to have paid.

Second way:

# Original code:
ggplot(titanic[titanic$Pclass == '3' & titanic$Embarked == 'S', ],
       aes(x = Fare)) +
  geom_density(fill = '#99d6ff', alpha=0.4) +
  geom_vline(aes(xintercept=median(Fare, na.rm=T)),
             colour='red', linetype='dashed', lwd=1) +
  scale_x_continuous(labels=dollar_format()) +
  theme_few()
(-> titanic
    (bra (r& (r== ($ titanic 'Pclass) 3)
             (r== ($ titanic 'Embarked) "S"))
         nil)
    (ggplot (aes :x 'Fare))
    (r+ (geom_density :fill "#99d6ff"
                      :alpha 0.4)
        (geom_vline (aes :xintercept
                         '(median Fare :na.rm true))
                    :colour "red"
                    :linetype "dashed"
                    :lwd 1)
        (scale_x_continuous :labels (dollar_format)))
    plot->svg)
0.00 0.05 0.10 0.15 $0 $20 $40 $60 Fare density

Tripathi: From this visualization, it seems quite reasonable to replace the NA Fare value with median for their class and embarkment which is $8.05.

Replace that NA with 8.05

# Original code:
titanic$Fare[1044] <- 8.05
summary(titanic$Fare)
(def titanic
  (bra<- titanic 1044 "Fare"
         8.05))
(-> titanic
     ($ 'Fare)
     summary)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   7.896  14.454  33.276  31.275 512.329 

Tripathi: Another way of Replace missing fare value with median fare for class/embarkment:

# Original code:
titanic$Fare[1044] <- median(titanic[titanic$Pclass == '3' & titanic$Embarked == 'S', ]$Fare, na.rm = TRUE)
(def titanic
  (bra<- titanic 1044 "Fare"
         (-> titanic
             (bra (r& (r== ($ titanic 'Pclass) 3)
                      (r== ($ titanic 'Embarked) "S"))
                  "Fare")
             (median :na.rm true))))

Missing Value in Age.

Tripathi: Show number of missing Age values.

# Original code:
sum(is.na(titanic$Age)) ```

before


::: {.sourceClojure}
```clojure
(-> titanic
    ($ 'Age)
    is-na
    sum)

:::

[1] 263

Tripathi: 263 passengers have no age listed. Taking a median age of all passengers doesn’t seem like the best way to solve this problem, so it may be easiest to try to predict the passengers’ age based on other known information.

To predict missing ages, I’m going to use the mice package. To start with I will factorize the factor variables and then perform mice(multiple imputation using chained equations).

Set a random seed:

# Original code:
set.seed(129)
(set-seed 129)
NULL

Tripathi: Perform mice imputation, excluding certain less-than-useful variables:

# Original code:
mice_mod <- mice(titanic[, !names(titanic) %in% c('PassengerId','Name','Ticket','Cabin','Family','Surname','Survived')], method='rf')
(def mice-mod
  (-> titanic
      (bra nil
           (-> titanic
               names
               (%in% ["PassengerId","Name","Ticket","Cabin","Family","Surname","Survived"])
               !))
      (mice :method "rf")))

Save the complete output.

# Original code:
mice_output <- complete(mice_mod)
(def mice-output
  (complete mice-mod))

Tripathi: Let’s compare the results we get with the original distribution of passenger ages to ensure that nothing has gone completely awry.

Plot age distributions:

# Original code:
par(mfrow=c(1,2))
hist(titanic$Age, freq=F, main='Age: Original Data',
     col='darkred', ylim=c(0,0.04))
hist(mice_output$Age, freq=F, main='Age: MICE Output',
     col='lightgreen', ylim=c(0,0.04))
(plot->svg
 (fn []
   (par :mfrow [1 2])
   (-> titanic
       ($ 'Age)
       (hist :freq 'F
             :main "Age: Original Data"
             :col "darkred"
             :lim [0 0.04]
             :xlab "Age"))
   (-> mice-output
       ($ 'Age)
       (hist :freq 'F
             :main "Age: MICE Output"
             :col "lightgreen"
             :lim [0 0.04]
             :xlab "Age"))))
Age: Original Data Age Density 0 20 40 60 80 0.00 0.01 0.02 0.03 Age: MICE Output Age Density 0 20 40 60 80 0.00 0.01 0.02 0.03 0.04

Tripathi: Things look good, so let’s replace our age vector in the original data with the output from the mice model.

Replace Age variable from the mice model:

# Original code:
titanic$Age <- mice_output$Age
(def titanic
  ($<- titanic 'Age
       ($ mice-output 'Age)))

Show new number of missing Age values

# Original code:
sum(is.na(titanic$Age))

after

(-> titanic
    ($ 'Age)
    is-na
    sum)
[1] 0

6.1 Feature Enginnering: Part 2

Tripathi: I will create a couple of new age-dependent variables: Child and Mother. A child will simply be someone under 18 years of age and a mother is a passenger who is 1) female, 2) is over 18, 3) has more than 0 children and 4) does not have the title ‘Miss’.

Relationship between age & survival: I include Sex since we know it’s a significant predictor.

# Original code:
ggplot(titanic[1:891,], aes(Age, fill = factor(Survived))) +
  geom_histogram() + facet_grid(.~Sex) + theme_few()
(-> titanic
    (bra (colon 1 891) nil)
    (ggplot (aes 'Age :fill '(factor Survived)))
    (r+ (geom_histogram)
        (facet_grid '(tilde . Sex))
        (theme_few))
    plot->svg)
female male 0 20 40 60 80 0 20 40 60 80 0 20 40 60 Age count factor(Survived) 0 1

Tripathi: Create the column Child, and indicate whether child or adult:

# Original code:
 titanic$Child[titanic$Age < 18] <- 'Child'
 titanic$Child[titanic$Age >= 18] <- 'Adult'
(def titanic
  (-> titanic
      (bra<- (r< ($ titanic 'Age) 18) "Child"
             "Child")
      (bra<- (r>= ($ titanic 'Age) 18) "Child"
             "Adult")))

Show counts:

# Original code:
table(titanic$Child, titanic$Survived)
(table ($ titanic 'Child)
       ($ titanic 'Survived))
       
          0   1
  Adult 484 273
  Child  65  69

Adding Mother variable:

# Original code:
titanic$Mother <- 'Not Mother'
titanic$Mother[titanic$Sex == 'female' & titanic$Parch >0 & titanic$Age > 18 & titanic$title != 'Miss'] <- 'Mother'
(def titanic
  (-> titanic
      ($<- 'Mother
           "Not Mother")
      (bra<- (reduce r&
                     [(r== ($ titanic 'Sex) "female")
                      (r> ($ titanic 'Parch) 0)
                      (r> ($ titanic 'Age) 18)
                      (r!= ($ titanic 'title) "Miss")])
             "Mother"
             "Mother")))

Show counts:

# Original code:
table(titanic$Mother, titanic$Survived)
(table ($ titanic 'Mother)
       ($ titanic 'Survived))
            
               0   1
  Mother      15  38
  Not Mother 534 304

Factorizing variables:

# Original code:
titanic$Child <- factor(titanic$Child)
titanic$Mother <- factor(titanic$Mother)
titanic$Pclass <- factor(titanic$Pclass)
titanic$Sex <- factor(titanic$Sex)
titanic$Embarked <- factor(titanic$Embarked)
titanic$Survived <- factor(titanic$Survived)
titanic$title <- factor(titanic$title)
titanic$fsizeD <- factor(titanic$fsizeD)
(def titanic
  (reduce (fn [data symbol]
            ($<- data symbol
                 (factor ($ data symbol))))
          titanic
          '[Child Mother Pclass Sex Embarked Survived title fsizeD]))

Check classes of all columns:

(lapply titanic (r "class"))
$PassengerId
[1] "integer"

$Survived
[1] "factor"

$Pclass
[1] "factor"

$Name
[1] "character"

$Sex
[1] "factor"

$Age
[1] "numeric"

$SibSp
[1] "integer"

$Parch
[1] "integer"

$Ticket
[1] "character"

$Fare
[1] "numeric"

$Cabin
[1] "character"

$Embarked
[1] "factor"

$title
[1] "factor"

$surname
[1] "character"

$famsize
[1] "integer"

$family
[1] "character"

$fsizeD
[1] "factor"

$deck
[1] "factor"

$Child
[1] "factor"

$Mother
[1] "factor"

7 Prediction

Split into training & test sets:

# Original code:
train <- titanic[1:891,]
test <- titanic[892:1309,]
(def train
  (bra titanic (colon 1 891) nil))
(def test
   (bra titanic (colon 892 1309) nil))

Building the model:

Tripathi: We then build our model using randomForest on the training set.

Set a random seed:

# Original code:
set.seed(754)
(set-seed 754)
NULL

Tripathi: Build the model (note: not all possible variables are used):

# Original code:
titanic_model <- randomForest(Survived ~ Pclass + Sex + Age + SibSp + Parch +
                                Fare + Embarked + title +
                                fsizeD + Child + Mother,
                              data = train)
(def titanic-model
  (randomForest '(tilde Survived
                        (+ Pclass Sex Age SibSp Parch
                           Fare Embarked title
                           fsizeD Child Mother))
                :data train))

Show model error:

# Original code:
plot(titanic_model, ylim=c(0,0.36))
legend('topright', colnames(titanic_model$err.rate), col=1:3, fill=1:3)
(plot->svg
 (fn []
   (plot titanic-model :ylim [0 0.36]
         :main "Model Error")
   (legend "topright"
           (colnames ($ titanic-model 'err.rate))
           :col (colon 1 3)
           :fill (colon 1 3))))
0 100 200 300 400 500 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 Model Error trees Error OOB 0 1

Tripathi: The black line shows the overall error rate which falls below 20%. The red and green lines show the error rate for ‘died’ and ‘survived’, respectively. We can see that right now we’re much more successful predicting death than we are survival.

7.1 Variable Importance

Get importance:

# Original code:
importance    <- importance(titanic_model)
varImportance <- data.frame(Variables = row.names(importance),
                            Importance = round(importance[ ,'MeanDecreaseGini'],2))
(importance titanic-model)
         MeanDecreaseGini
Pclass          29.479780
Sex             61.243379
Age             46.395186
SibSp           11.921763
Parch            7.785026
Fare            59.766976
Embarked         9.710976
title           65.424749
fsizeD          17.363445
Child            4.306876
Mother           2.240933
(def importance-info
  (importance titanic-model))
(def var-importance
  (data-frame :Variables (row-names importance-info)
              :Importance (-> importance-info
                              (bra nil "MeanDecreaseGini")
                              round)))
importance-info
         MeanDecreaseGini
Pclass          29.479780
Sex             61.243379
Age             46.395186
SibSp           11.921763
Parch            7.785026
Fare            59.766976
Embarked         9.710976
title           65.424749
fsizeD          17.363445
Child            4.306876
Mother           2.240933
var-importance
         Variables Importance
Pclass      Pclass         29
Sex            Sex         61
Age            Age         46
SibSp        SibSp         12
Parch        Parch          8
Fare          Fare         60
Embarked  Embarked         10
title        title         65
fsizeD      fsizeD         17
Child        Child          4
Mother      Mother          2

7.2 Variable importance

Create a rank variable based on importance:

# Original code:
rankImportance <- varImportance %>%
   mutate(Rank = paste0('#',dense_rank(desc(Importance))))
(def rank-importance
  (-> var-importance
      (mutate :Rank '(paste0 "#" (dense_rank (desc Importance))))))
rank-importance
         Variables Importance Rank
Pclass      Pclass         29   #5
Sex            Sex         61   #2
Age            Age         46   #4
SibSp        SibSp         12   #7
Parch        Parch          8   #9
Fare          Fare         60   #3
Embarked  Embarked         10   #8
title        title         65   #1
fsizeD      fsizeD         17   #6
Child        Child          4  #10
Mother      Mother          2  #11

Tripathi: Use ggplot2 to visualize the relative importance of variables

# Original code:
ggplot(rankImportance, aes(x = reorder(Variables, Importance),
                           y = Importance, fill = Importance)) +
  geom_bar(stat='identity') +
  geom_text(aes(x = Variables, y = 0.5, label = Rank),
            hjust=0, vjust=0.55, size = 4, colour = 'red') +
  labs(x = 'Variables') +
  coord_flip() +
  theme_few()
(-> rank-importance
    (ggplot (aes :x '(reorder Variables Importance)
                 :y 'Importance
                 :fill 'Importance))
    (r+ (geom_bar :stat "Identity")
        (geom_text (aes :x 'Variables
                        :y 0.5
                        :label 'Rank)
                   :hjust 0
                   :vjust 0.55
                   :size 4
                   :colour "red")
        (labs :x "Variables")
        (coord_flip)
        (theme_few))
    plot->svg)
#5 #2 #4 #7 #9 #3 #8 #1 #6 #10 #11 Mother Child Parch Embarked SibSp fsizeD Pclass Age Fare Sex title 0 20 40 60 Importance Variables 20 40 60 Importance

Tripathi: From the plot we can see that the ‘title’ variable has the highest relative importance out of all of our predictor variables.

7.3 Final Prediction

Predict using the test set:

# Original code:
prediction <- predict(titanic_model, test)
prediction
(def prediction
  (predict titanic-model test))

Tripathi: Save the solution to a dataframe with two columns: PassengerId and Survived (prediction).

# Original code:
Output<- data.frame(PassengerID = test$PassengerId, Survived = prediction)
Output
(def output (data-frame :PassengerId ($ test 'PassengerId)
                        :Survived prediction))
(r->clj output)

_unnamed [418 3]:

:$row.names :PassengerId :Survived
892 892 :0
893 893 :0
894 894 :0
895 895 :0
896 896 :0
897 897 :0
898 898 :1
899 899 :0
900 900 :1
901 901 :0
1299 1299 :0
1300 1300 :1
1301 1301 :1
1302 1302 :1
1303 1303 :1
1304 1304 :0
1305 1305 :0
1306 1306 :1
1307 1307 :0
1308 1308 :0
1309 1309 :1

Write the Output to file:

# Original code:
write.csv(Output, file = 'pradeep_titanic_output.csv', row.names = F)
(write-csv output
           :file "/tmp/pradeep_titanic_output.csv"
           :row.names 'F)
NULL

7.4 Conclusion

Tripathi: Thank you for taking the time to read through my first exploration of a Titanic Kaggle dataset. Again, this newbie welcomes comments and suggestions!

source: notebooks/clojisr/v1/tutorials/titanic0.clj