Gara kaggle Titanic

Questo è il mio primo script Kaggle. Questo script sarà utile per muovere i primi passi all’interno di R e delle sue funzioni, allo stesso tempo però verranno esaminate anche tecniche più avanzate di machine learning in modo da stimolare la curiosità e dare l’opportunità di avere una panoramica più ampia delle potenzialità di R. Quindi dopo una parte iniziale di manipolazione e visualizzazione del dataset ho proceduto sino all’utilizzo di una tecnica di Stacking classifiers o un classificatore di classificatori per creare un modello di predizione dei sopravvissuti nel Titanic. Ho scelto di lavorare con il dataset di Titanic dopo aver preso spunto da alcuni scripts fatti da atri Kagglers come ispirazione. Le finalità di questo script sono didattiche, sono quindi ben accetti eventuali feedback. Questo script infatti è stato scritto in italiano con l’intento di favorire l’apprendimento e la diffusione della conoscenza di temi inerenti la data science, pensando perciò di fare una cosa gradita agli utenti interessati.

Lo script è composto di tre parti:

In R è necessario caricare i pacchetti che ci servono nella nostra elaborazione. Per caricare i pacchetti dovrei averli prima installati, scaricare continuamente pacchetti può sembrare un esercizio ripetitivo, perciò posso utilizzare una formula del genere per rendere il tutto meno noioso:

if (! ("ggplot2" %in% rownames(installed.packages()))) { install.packages("ggplot2") }
if (! ("dplyr" %in% rownames(installed.packages()))) { install.packages("dplyr") }
if (! ("ggthemes" %in% rownames(installed.packages()))) { install.packages("ggthemes") }
if (! ("scales" %in% rownames(installed.packages()))) { install.packages("scales") }
if (! ("mice" %in% rownames(installed.packages()))) { install.packages("mice") }
if (! ("randomForest" %in% rownames(installed.packages()))) { install.packages("randomForest") }
if (! ("caret" %in% rownames(installed.packages()))) { install.packages("caret") }
if (! ("e1071" %in% rownames(installed.packages()))) { install.packages("e1071") }
if (! ("tree" %in% rownames(installed.packages()))) { install.packages("tree") }

Carico e controllo i dati

# Carico i pacchetti
library('ggplot2') # visualizazione
library('ggthemes') # visualizazione
library('scales') # visualizazione
library('dplyr') #manipolaizone
library('mice') # modellizzazione
library('randomForest') # algoritmo di classificazione
library('caret') # algoritmo di classificazione e altro
library(e1071)# algoritmo di classificazione
library(tree)# algoritmo di classificazione

Ora che abbiamo caricato le librerie, carichiamo ed osserviamo i dati

train <- read.csv('E:/gara kaggle titanic/train.csv', stringsAsFactors = F)
test  <- read.csv('E:/gara kaggle titanic/test.csv', stringsAsFactors = F)

full  <- bind_rows(train, test) # metto insieme train e test, tanto devo pulire i dati

# osserviamo i dati
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" ...

Dopo aver osservato le mie variabili, la loro classe e tipo, le prime righe di ogni caratteristca.Stiamo lavorando con 1309 osservazioni e 12 variabili. Per rendere le cose più esplicite riportiamo le varibili con una breve descrizione:

Nome Variabil Description
Survived Supravvisuto (1) o morti (0)
Pclass Classe dei passeggeri
Name Nome del passeggero
Sex Sesso del passeggero
Age Età del passeggero
SibSp Numero di fratelli/spose a bordo
Parch Numero di genitori/figli a bordo
Ticket Numero del biglietto
Fare Tariffa
Cabin Cabina
Embarked Porto di imbarco

Manipolazione delle caratteristiche del dataset

Cosa può significare un nome?

La prima colonna che osservo è la nome del passeggero questo perchè è una variabile che può essere divisa in più parti significative che possono essere utili per la costruzione di ulteriori elementi utili o utilizzate direttamente come variabili. Per esempio , titolo del passeggero è contenuta con la variabile nome del passeggero inoltre si può utilizzare la caratteristica cognome per rappresentare le famiglie. Cominciamo ora con il manipolare caratteristiche!

Breve digressione su gsub e le regular expression

Il nostro scopo è prendere il cognome da un nome esteso scritto in questo modo: Holthen, Mr. Johan Martin
Quindi useremo la funzione gsub e ci aiuteremo con alcuni esempi per capire meglio come funziona, intanto è utile sapere che se vogliamo utilizzare l’help ci basta far girare ? ed il nome della funzione che ci interessa (nel nostro caso ?gsub).

primo esempio prendiamo more e lo facciamo diventare cuore, sostituiamo m con cu :

x <- "more"
gsub("m","cu",x)
## [1] "cuore"

con questo comando ’\..*’ sostituiamo tutto ciò che c’è dopo del punto compreso il punto:

pippo <- 'Holthen, Mr. Johan Martin'

prova1 <- gsub('(\\..*)', '', pippo)
prova1
## [1] "Holthen, Mr"

con questo comando ’\,.*’ sostituiamo tutto ciò che c’è dopo la virgola compresa la virgola:

prova2 <- gsub('(\\,.*)', '', pippo)
prova2
## [1] "Holthen"

con questo comando ’.*, ’ sostituiamo tutto ciò che c’è prima dello spazio compreso lo spazio:

prova3 <- gsub('(\\s.*)', '', pippo)
prova3
## [1] "Holthen,"

con questo comando ’.*, ’ sostituiamo tutto ciò che c’è prima della virgola compresa la virgola:

prova4 <- gsub('(.*, )', '', pippo)
prova4
## [1] "Mr. Johan Martin"

Adesso che abbiamo un poco più chiaro il funzionamento di gsub possiamo proseguire

# recupero del titolo dal nome del passeggero
full$Title <- gsub('(.*, )|(\\..*)', '', full$Name)

Poniamo un attimo attenzione a quello che abbiamo appena fatto. full è il nostro dataset, con $ selezioniamo la colonna Title al suo interno. La questione è che nel nostro dataset non c’è una colonna che si chiama Title. Ossia con questo comando R aggiunge direttamente questa nuova colonna al dataset, quindi non abbiamo bisogno di creare una colonna vuota che poi andiamo a riempire con dei dati, possimo direttamte dirgli metti questo in questa colonna ed R se non trova una colonna con quel nominativo ne aggiunge una nuova in cui scrive i dati che ci interessano.

#Con questo grafico il conteggio dei titoli per sesso
table(full$Sex, full$Title)
##         
##          Capt Col Don Dona  Dr Jonkheer Lady Major Master Miss Mlle Mme
##   female    0   0   0    1   1        0    1     0      0  260    2   1
##   male      1   4   1    0   7        1    0     2     61    0    0   0
##         
##           Mr Mrs  Ms Rev Sir the Countess
##   female   0 197   2   0   0            1
##   male   757   0   0   8   1            0
# I titoli che compaiono poche volte li raccolgo insieme nel gruppo rare
rare_title <- c('Dona', 'Lady', 'the Countess','Capt', 'Col', 'Don', 'Dr', 'Major', 'Rev', 'Sir', 'Jonkheer')

# Riassegno anche i titoli con refuso tipo mlle, ms, e mme
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'

# Rivediamo i titoli ordinati per sex
table(full$Sex, full$Title)
##         
##          Master Miss  Mr Mrs Rare Title
##   female      0  264   0 198          4
##   male       61    0 757   0         25
# Finalmente ricaviamo il cognome 
# per sapply, per ora, diciamo solo che è un modo per fare cicli for in R
# ossia prende una funziona e la fa girare sugli elementi di un dataset 
full$Surname <- sapply(full$Name,  
                      function(x) strsplit(x, split = '[,.]')[[1]][1])
cat(paste('We have <b>', nlevels(factor(full$Surname)), '</b> unique surnames. I would be interested to infer ethnicity based on surname --- another time.'))

We have 875 unique surnames. I would be interested to infer ethnicity based on surname — another time.

I componenti della stessa famiglia si aiutano tra loro?

Questa può essere una buona intuizione, ed apre una nuova via di intepretazione dei dati. Ossia i dati rappresentano anche cose della nostra esistenza, quindi possiamo armarci di interpretazione per comprendere più affondo cosa li muove o le possibili relazioni che ci sono. In altre parole trarre valore dai dati significa anche questo, porsi domande semplici: I componenti della stessa famiglia si aiutano tra loro?

# Creiamo la variabole dimensione famiglia + la persona stessa che stiamo analizzando
# ricordo che SibSp sarebbe numero di fratelli + marito o moglie mentre Prach sta per 
# genitori e figli.
 
full$Fsize <- full$SibSp + full$Parch + 1

# Creiamo una variabile con cognome e numero di componenti
full$Family <- paste(full$Surname, full$Fsize, sep='_')

A questo punto proviamo ad osservare tramite un grafico se effettivamente otteniamo una correlazione tra il numero di componenti di un nucleo famigliare ed i soppravvissuti.

# Adesso utiliziamo ggplot per fare il grafico
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()

Si può osservare che c’è una correlazione tra la dimensione delle famiglie e la probabilità di sopravvivere.

# Trasformiamo la variabile famiglia in una variabile discreta
full$FsizeD[full$Fsize == 1] <- 'singleton'
full$FsizeD[full$Fsize < 5 & full$Fsize > 1] <- 'small'
full$FsizeD[full$Fsize > 4] <- 'large'

Valori mancanti

Quando abbiamo dei dati mancanti come possiamo comportarci? Alcuni potrebbero dire di sotituirli con la media ed infatti se ho pochi valori mancanti posso farlo:

# recupero dei dati mancanti tramite la media
full$Fare[1044] <- median(full[full$Pclass == '3' & full$Embarked == 'S', ]$Fare, na.rm = TRUE)

Ora descriveremo una seconda possibilità, ossia usare i dati che abbimo a disposizione, le altre varibili disponibili o colonne, per costruire un modello di predizione in modo da ricavare i dati mancanti.

Analisi predittiva

Ora procediamo osservando il numero di valori mancanti. In particolare analizzeremo la variabile Age, come primo passo cerchiamo di capire quanti sono i valori che non hanno valore.

# Mostra il numero di valori mancanti
sum(is.na(full$Age))
## [1] 263

Non useremo rpart (recursive partitioning for regression) per predire i valori Age mancanti, ma useremo mice un pacchetto apposito giusto per provare una cosa divesa. Per approfondimenti clicca qui (PDF). Dato che non lo abbiamo fatto prima faremo diventare factor le variabili, questo si fa per ottimizioamo l’uso della memoria di R.

# Raccogliamo tutte le varibili insieme, e trasformiamole in factor
factor_vars <- c('PassengerId','Pclass','Sex',
                 'Title','Surname','Family','FsizeD')

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

# poniamo un valore random
set.seed(129)

# Adesso faremo la nostr predizione togliendo (abbimo messo il punto esclamativo che
# significa negazione) le varibili che hanno poco peso per determinare l'età, come nome cognome identificativo ecc. :

mice_mod <- mice(full[, !names(full) %in% c('PassengerId','Name','Ticket','Cabin','Family','Surname','Survived')], method='rf') 
## 
##  iter imp variable
##   1   1  Age
##   1   2  Age
##   1   3  Age
##   1   4  Age
##   1   5  Age
##   2   1  Age
##   2   2  Age
##   2   3  Age
##   2   4  Age
##   2   5  Age
##   3   1  Age
##   3   2  Age
##   3   3  Age
##   3   4  Age
##   3   5  Age
##   4   1  Age
##   4   2  Age
##   4   3  Age
##   4   4  Age
##   4   5  Age
##   5   1  Age
##   5   2  Age
##   5   3  Age
##   5   4  Age
##   5   5  Age
# Salviamo il risultato 
mice_output <- complete(mice_mod)

Andiamo a comparare le di stribuzioni di dati originali con quelli da noi creati

# Grafico della distribuzione di Age
par(mfrow=c(1,2))
hist(full$Age, freq=F, main='Age: Dati originali', 
  col='darkgreen', ylim=c(0,0.04))
hist(mice_output$Age, freq=F, main='Age: Dati modello', 
  col='lightgreen', ylim=c(0,0.04))

Sembra che le cose vadano bene. Adesso andiamo ad inserire i dati creati con il modello mice nella nostra colonna Age.

# Riscriviamo la varibile Age.
full$Age <- mice_output$Age

# Verifichiamo che la somma dei valori persi sia zero
sum(is.na(full$Age))
## [1] 0

Modello predittivo

Infine siamo in grado di predire chi sopravviverà tra i passeggeri del Titanic, ora utilizzeremo una randomForest classification algorithm.

Dividiamo i nostri dati

come primo passo dividiamo i nostri dati in train e test.

# dividiamo i dati come erano in origine
train <- full[1:891,]
test <- full[892:1309,]

Costruzione del modello

Dapprimo costruiamo la nostra randomForest sul training set.

# settiamo il valore random
set.seed(754)

# Costruiamo il nostro modello
rf_model <- randomForest(factor(Survived) ~ Pclass + Sex + Age + SibSp + Parch + Title + 
                                            FsizeD,
                                            data = train)

# osserviamo l'errore nel modello
plot(rf_model, ylim=c(0,0.36))
legend('topright', colnames(rf_model$err.rate), col=1:3, fill=1:3)

La linea nera mostra l’errore complessivo che commettiamo, che cala intorno al 20%.La linea rossa e la linea verde mostrano l’errore per ‘died’ e ‘survived’. Quindi possiamo dire se siamo più bravi a predire ‘died’ oppure ‘survived’ a seconda della posizione che occupano nel grafico.

Stacking classifiers - utliziamo piu modellei insieme

La prima cosa che faremo è costruire il nostro modello con i dati del train

#ripartiamo dalla random forest che abbiamo già visto
rf1 <- randomForest(factor(Survived) ~ Pclass + Sex + Age + SibSp + Parch + Title + 
                      FsizeD,
                    data = train)
#utiliziamo un modello di regressione lineare
glm1    <-  glm(factor(Survived) ~ Pclass + Sex + Age + SibSp + Parch + Title + 
              FsizeD,
            data = train,   family = 'binomial')
#un modello di classificazione ad albero
ct1 <-  tree(factor(Survived) ~ Pclass + Sex + Age + SibSp + Parch + Title + 
              FsizeD,data = train)
#uno di classificazione naive baise
nb1 <-naiveBayes(factor(Survived) ~ Pclass + Sex + Age + SibSp + Parch + Title + 
                   FsizeD,data = train)

#adesso andiamo a fare le nostre predizioni sul train
pre.glm1 <-  predict(glm1,  train, type = 'response')
pre.nb1 <-  predict(nb1,    train)
pre.ct1 <-  predict(ct1,    train, type = 'class')
pre.rf1 <- predict(rf1, train, type = 'class')

#adesso possimo mettere tutte le nostre predizioni in un unico file
sample.predictions <- data.frame(pre.glm1  ,pre.rf1 ,pre.nb1 ,pre.ct1 )

#questa cosa ci serve perchè stiamo per costruire un modello e abbiamo bisogno dell'output
sample.predictions$Survived <- train$Survived

#a questo punto possiamo utilizzare ancora una random forest che fa una previsione sul #file costruito per mezzo di altri modelli predittivi
rf_stack <- randomForest(Survived ~.,
                         data = sample.predictions)

Abbiamo costruito il nostro modello di predizione basato su più modelli. Poniamo attenzione al fatto che abbiamo semplicemente inserito i dati di output di più modelli che poi abbiamo utlizzato come input di un ulteriore modello.

#a questo punto possiamo proseguire facendo la nostra predizione con i dati del test
pre.glm1 <-  predict(glm1,  test, type = 'response')
pre.nb1 <-  predict(nb1,    test)
pre.ct1 <-  predict(ct1,    test, type = 'class')
pre.rf1 <- predict(rf1, test, type = 'class')

# accorpiamo le preidizioni dei nostri modelli in un unico file
sample.predictions <- data.frame(pre.glm1  ,pre.rf1 ,pre.nb1 ,pre.ct1 )

Predizione

Infiene siamo pronti per procere con la nostra predizione. Una volta che lo abbiamo fatto possiamo riprendere i punti precendenti e aggiungere delle variabili nel modello o utilizzare calissificatori differenti, ossia questo potebbe essere un ponto di partenza e non l’ultimo step.

# procediamo andando a fare la predizione inserendo nel modello questo file di dati
prediction  <-  predict(rf_stack,   sample.predictions, type='class')
prediction[(prediction > 0) & (prediction <= 0.5)] <- 0
prediction[(prediction > 0.5) & (prediction <= 1)] <- 1
prediction

# Adesso salviamo i dati previsti in un file con le seguenti colonne: PassengerId and Survived (predetti)
solution <- data.frame(PassengerID = test$PassengerId, Survived = prediction)

# scriviamo la soluzione in un file csv
write.csv(solution, file = 'rf_mod_Solution.csv', row.names = F)

Conclusioni

Spero che questo scritto possa essere uno strumento per stimolare la curiosità intorno al mondo della data science e favorire l’apprendimento, anche tramite l’intuizione ed il piacere di acquisire nuove conoscenze, consapevoli che si sta lavorando con strumenti molto utili e trasversali. Sono benvenuti commenti e suggerimenti! Grazie per il vostro tempo.