1.1. Load and check data

# Load packages
library('ggplot2') # visualization
library('ggthemes') # visualization
library('scales') # visualization
library('dplyr') # data manipulation
library('mice') # imputation
library('randomForest') # classification algorithm

train <- read.csv('train.csv')
test  <- read.csv('test.csv')

full  <- bind_rows(train, test) # bind training & test data

# check data
str(full)

## '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" ...

Dapat memahami variabel-variabel, jenis kelasnya, dan beberapa pengamatan awal dari masing-masing. Dengan 1309 pengamatan dari 12 variabel. Untuk membuatnya lebih eksplisit :

2.1. What’s in a name?

Variabelnama penumpang dapat dipecah menjadi beberapa variabel tambahan yang bermakna yang dapat digunakan untuk membuat prediksi atau digunakan dalam pembuatan variabel baru. Sebagai contoh, gelar penumpang terkandung dalam variabel nama penumpang dan dapat menggunakan nama keluarga untuk mewakili keluarga.

# Grab title from passenger names
full$Title <- gsub('(.*, )|(\\..*)', '', full$Name)

# Show title counts by sex
table(full$Sex, full$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

# Titles with very low cell counts to be combined to "rare" level
rare_title <- c('Dona', 'Lady', 'the Countess','Capt', 'Col', 'Don', 
                'Dr', 'Major', 'Rev', 'Sir', 'Jonkheer')

# Also reassign mlle, ms, and mme accordingly
full$Title[full$Title == 'Mlle']        <- 'Miss' 
full$Title[full$Title == 'Ms']          <- 'Miss'
full$Title[full$Title == 'Mme']         <- 'Mrs' 
full$Title[full$Title %in% rare_title]  <- 'Rare Title'

# Show title counts by sex again
table(full$Sex, full$Title)

##         
##          Master Miss  Mr Mrs Rare Title
##   female      0  264   0 198          4
##   male       61    0 757   0         25

# Finally, grab surname from passenger name
full$Surname <- sapply(full$Name,  
                      function(x) strsplit(x, split = '[,.]')[[1]][1])

cat(paste('Kami memiliki <b>', nlevels(factor(full$Surname)), '</b> nama keluarga yang unik. Saya tertarik untuk menyimpulkan etnisitas berdasarkan nama keluarga --- di lain waktu.'))

## Kami memiliki <b> 875 </b> nama keluarga yang unik. Saya tertarik untuk menyimpulkan etnisitas berdasarkan nama keluarga --- di lain waktu.

2.2. Do families sink or swim together?

Selanjutnya adalah denganmembuat beberapa variabel keluarga baru dengan membuat variabel ukuran keluarga berdasarkan jumlah saudara kandung/pasangan (mungkin ada yang memiliki lebih dari satu pasangan?) dan jumlah anak/orang tua

# Create a family size variable including the passenger themselves
full$Fsize <- full$SibSp + full$Parch + 1

# Create a family variable 
full$Family <- paste(full$Surname, full$Fsize, sep='_')

Seperti apa variabel ukuran keluarga kita? Untuk membantu kita memahami bagaimana variabel ini berhubungan dengan kelangsungan hidup,

# Use ggplot2 to visualize the relationship between family size & survival
ggplot(full[1:891,], aes(x = Fsize, fill = factor(Survived))) +
  geom_bar(stat='count', position='dodge') +
  scale_x_continuous(breaks=c(1:11)) +
  labs(x = 'Family Size') +
  theme_few()

Kita dapat melihat bahwa ada penalti untuk bertahan hidup bagi lajang dan mereka yang memiliki jumlah anggota keluarga di atas 4 orang. Kita dapat memecah variabel ini menjadi tiga level yang akan sangat membantu karena jumlah keluarga besar relatif lebih sedikit. Mari kita buat variabel ukuran keluarga yang didiskritisasi.

# Discretize family size
full$FsizeD[full$Fsize == 1] <- 'singleton'
full$FsizeD[full$Fsize < 5 & full$Fsize > 1] <- 'small'
full$FsizeD[full$Fsize > 4] <- 'large'

# Show family size by survival using a mosaic plot
mosaicplot(table(full$FsizeD, full$Survived), main='Family Size by Survival', shade=TRUE)

Plot mosaik menunjukkan bahwa kita mempertahankan aturan kita bahwa ada penalti untuk bertahan hidup di antara lajang dan keluarga besar, tetapi ada keuntungan bagi penumpang dalam keluarga kecil. Saya ingin melakukan sesuatu yang lebih jauh dengan variabel usia kita, tetapi 263 baris memiliki nilai usia yang hilang, jadi kita harus menunggu sampai setelah kita mengatasi kehilangan tersebut.

2.3. Treat a few more variables

Apa yang tersisa? Mungkin ada beberapa informasi yang mungkin berguna dalam variabel kabin penumpang, termasuk tentang deckk mereka.

# This variable appears to have a lot of missing values
full$Cabin[1:28]

##  [1] ""            "C85"         ""            "C123"        ""           
##  [6] ""            "E46"         ""            ""            ""           
## [11] "G6"          "C103"        ""            ""            ""           
## [16] ""            ""            ""            ""            ""           
## [21] ""            "D56"         ""            "A6"          ""           
## [26] ""            ""            "C23 C25 C27"

# The first character is the deck. For example:
strsplit(full$Cabin[2], NULL)[[1]]

## [1] "C" "8" "5"

# Create a Deck variable. Get passenger deck A - F:
full$Deck<-factor(sapply(full$Cabin, function(x) strsplit(x, NULL)[[1]][1]))

Ada banyak hal yang mungkin dapat dilakukan di sini, termasuk mencari kabin dengan beberapa kamar yang terdaftar (misalnya, baris 28: “C23 C25 C27”), tetapi mengingat sedikitnya kolom, kami akan berhenti di sini.

3.1. Sensible value imputation

# Passengers 62 and 830 are missing Embarkment
full[c(62, 830), 'Embarked']

## [1] "" ""

cat(paste('We will infer their values for **embarkment** based on present data that we can imagine may be relevant: **passenger class** and **fare**. We see that they paid<b> $', full[c(62, 830), 'Fare'][[1]][1], '</b>and<b> $', full[c(62, 830), 'Fare'][[1]][2], '</b>respectively and their classes are<b>', full[c(62, 830), 'Pclass'][[1]][1], '</b>and<b>', full[c(62, 830), 'Pclass'][[1]][2], '</b>. So from where did they embark?'))

## We will infer their values for **embarkment** based on present data that we can imagine may be relevant: **passenger class** and **fare**. We see that they paid<b> $ 80 </b>and<b> $ NA </b>respectively and their classes are<b> 1 </b>and<b> NA </b>. So from where did they embark?

# Get rid of our missing passenger IDs
embark_fare <- full %>%
  filter(PassengerId != 62 & PassengerId != 830)

# Use ggplot2 to visualize embarkment, passenger class, & median fare
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()

Tarif rata-rata untuk penumpang kelas satu yang berangkat dari Charbourg (‘C’) sama dengan $80 yang dikeluarkan Saya rasa kita dapat dengan aman mengganti nilai NA dengan ‘C’.

# Since their fare was $80 for 1st class, they most likely embarked from 'C'
full$Embarked[c(62, 830)] <- 'C'

Kami hampir memperbaiki beberapa nilai NA di sana-sini. Penumpang di baris 1044 memiliki nilai Tarif NA.

# Show row 1044
full[1044, ]

Ini adalah penumpang kelas tiga yang berangkat dari Southampton (‘S’). Mari kita bayangkan tarif di antara penumpang lain yang berbagi kelas dan keberangkatan (n = 494).

ggplot(full[full$Pclass == '3' & full$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()

Dari visualisasi ini, tampaknya cukup masuk akal untuk mengganti nilai Tarif NA dengan nilai median untuk kelas dan keberangkatan mereka, yaitu $8.05.

# Replace missing fare value with median fare for class/embarkment
full$Fare[1044] <- median(full[full$Pclass == '3' & full$Embarked == 'S', ]$Fare, na.rm = TRUE)

3.2. Predictive imputation

Terakhir, seperti yang telah kita catat sebelumnya, ada beberapa nilai Usia yang hilang dalam data kita. Kita akan sedikit lebih canggih dalam mengimputasi nilai usia yang hilang. Mengapa? Karena kita bisa. Kita akan membuat model yang memprediksi usia berdasarkan variabel lain.

# Show number of missing Age values
sum(is.na(full$Age))

## [1] 263

Kita tentu saja dapat menggunakan rpart (partisi rekursif untuk regresi) untuk memprediksi usia yang hilang, tetapi saya akan menggunakan paket mice untuk tugas ini untuk sesuatu yang berbeda. Anda dapat membaca lebih lanjut tentang imputasi berganda menggunakan persamaan berantai dalam r di sini (PDF). Karena kita belum melakukannya, saya akan memfaktorkan variabel faktor terlebih dahulu dan kemudian melakukan imputasi tikus.

# Make variables factors into factors
factor_vars <- c('PassengerId','Pclass','Sex','Embarked',
                 'Title','Surname','Family','FsizeD')

full[factor_vars] <- lapply(full[factor_vars], function(x) as.factor(x))

# Set a random seed
set.seed(129)

# Perform mice imputation, excluding certain less-than-useful variables:
mice_mod <- mice(full[, !names(full) %in% c('PassengerId','Name','Ticket','Cabin','Family','Surname','Survived')], method='rf')

## 
##  iter imp variable
##   1   1  Age  Deck
##   1   2  Age  Deck
##   1   3  Age  Deck
##   1   4  Age  Deck
##   1   5  Age  Deck
##   2   1  Age  Deck
##   2   2  Age  Deck
##   2   3  Age  Deck
##   2   4  Age  Deck
##   2   5  Age  Deck
##   3   1  Age  Deck
##   3   2  Age  Deck
##   3   3  Age  Deck
##   3   4  Age  Deck
##   3   5  Age  Deck
##   4   1  Age  Deck
##   4   2  Age  Deck
##   4   3  Age  Deck
##   4   4  Age  Deck
##   4   5  Age  Deck
##   5   1  Age  Deck
##   5   2  Age  Deck
##   5   3  Age  Deck
##   5   4  Age  Deck
##   5   5  Age  Deck

# Save the complete output 
mice_output <- complete(mice_mod)

Mari bandingkan hasil yang kami peroleh dengan distribusi usia penumpang yang asli untuk memastikan tidak ada yang salah.

# Plot age distributions
par(mfrow=c(1,2))
hist(full$Age, freq=F, main='Age: Original Data', 
  col='darkgreen', ylim=c(0,0.04))
hist(mice_output$Age, freq=F, main='Age: MICE Output', 
  col='lightgreen', ylim=c(0,0.04))

Semuanya terlihat baik, jadi mari kita ganti vektor usia pada data asli dengan output dari model tikus.

# Replace Age variable from the mice model.
full$Age <- mice_output$Age

# Show new number of missing Age values
sum(is.na(full$Age))

## [1] 0

Kita sudah selesai memasukkan nilai untuk semua variabel yang kita pedulikan untuk saat ini! Sekarang kita telah memiliki variabel Usia yang lengkap, hanya ada beberapa sentuhan akhir yang perlu dilakukan. Kita dapat menggunakan Age untuk melakukan sedikit rekayasa fitur.

3.3. Feature Engineering: Round 2

Setelah kita mengetahui usia setiap orang, kita dapat membuat beberapa variabel baru yang bergantung pada usia: Anak dan Ibu. Seorang anak adalah seseorang yang berusia di bawah 18 tahun dan seorang ibu adalah penumpang yang 1) berjenis kelamin perempuan, 2) berusia di atas 18 tahun, 3) memiliki lebih dari 0 anak (bukan main-main!), dan 4) tidak memiliki gelar ‘Nona’.

# First we'll look at the relationship between age & survival
ggplot(full[1:891,], aes(Age, fill = factor(Survived))) + 
  geom_histogram() + 
  # I include Sex since we know (a priori) it's a significant predictor
  facet_grid(.~Sex) + 
  theme_few()

# Create the column child, and indicate whether child or adult
full$Child[full$Age < 18] <- 'Child'
full$Child[full$Age >= 18] <- 'Adult'

# Show counts
table(full$Child, full$Survived)

##        
##           0   1
##   Adult 481 268
##   Child  68  74

Sepertinya menjadi seorang anak tidak ada salahnya, tetapi itu juga tidak akan menyelamatkan Anda! Kita akan menyelesaikan rekayasa fitur kita dengan menciptakan variabel Ibu. Mungkin kita bisa berharap bahwa ibu lebih mungkin untuk selamat dari Titanic.

# Adding Mother variable
full$Mother <- 'Not Mother'
full$Mother[full$Sex == 'female' & full$Parch > 0 & full$Age > 18 & full$Title != 'Miss'] <- 'Mother'

# Show counts
table(full$Mother, full$Survived)

##             
##                0   1
##   Mother      16  39
##   Not Mother 533 303

# Finish by factorizing our two new factor variables
full$Child  <- factor(full$Child)
full$Mother <- factor(full$Mother)

Semua variabel yang kita pedulikan harus diperhatikan dan tidak boleh ada data yang hilang.Untuk memastikannya dapat diperiksa ulang

md.pattern(full)

##     PassengerId Pclass Name Sex Age SibSp Parch Ticket Fare Cabin Embarked
## 204           1      1    1   1   1     1     1      1    1     1        1
## 687           1      1    1   1   1     1     1      1    1     1        1
## 91            1      1    1   1   1     1     1      1    1     1        1
## 327           1      1    1   1   1     1     1      1    1     1        1
##               0      0    0   0   0     0     0      0    0     0        0
##     Title Surname Fsize Family FsizeD Child Mother Survived Deck     
## 204     1       1     1      1      1     1      1        1    1    0
## 687     1       1     1      1      1     1      1        1    0    1
## 91      1       1     1      1      1     1      1        0    1    1
## 327     1       1     1      1      1     1      1        0    0    2
##         0       0     0      0      0     0      0      418 1014 1432

4.2. Building the model

Kemudian dapat membangun model menggunakan randomForest pada set pelatihan.

# Set a random seed
set.seed(754)

# Build the model (note: not all possible variables are used)
rf_model <- randomForest(factor(Survived) ~ Pclass + Sex + Age + SibSp + Parch + 
                                            Fare + Embarked + Title + 
                                            FsizeD + Child + Mother,
                                            data = train)

# Show model error
plot(rf_model, ylim=c(0,0.36))
legend('topright', colnames(rf_model$err.rate), col=1:3, fill=1:3)

Garis hitam menunjukkan tingkat kesalahan secara keseluruhan yang berada di bawah 20%. Garis merah dan hijau menunjukkan tingkat kesalahan untuk ‘meninggal’ dan ‘Survive’.

4.3. Variable importance

Selanjutnya dengan melihat tingkat kepentingan variabel relatif dengan memplotkan rata-rata penurunan Gini yang dihitung di semua pohon.

# Get importance
importance    <- importance(rf_model)
varImportance <- data.frame(Variables = row.names(importance), 
                            Importance = round(importance[ ,'MeanDecreaseGini'],2))

# Create a rank variable based on importance
rankImportance <- varImportance %>%
  mutate(Rank = paste0('#',dense_rank(desc(Importance))))

# Use ggplot2 to visualize the relative importance of variables
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()

4.4 Prediction!

langkah terakhir adalah dengan membuat prediksi namun kita mengulangi langkah-langkah sebelumnya dengan melakukan penyesuaian atau menyesuaikan data menggunakan model yang berbeda atau menggunakan kombinasi variabel yang berbeda untuk mencapai prediksi yang lebih baik.

# Predict using the test set
prediction <- predict(rf_model, test)

# Save the solution to a dataframe with two columns: PassengerId and Survived (prediction)
solution <- data.frame(PassengerID = test$PassengerId, Survived = prediction)

# Write the solution to file
write.csv(solution, file = 'rf_mod_Solution.csv', row.names = F)