Si te preguntas si habrías sobrevivido al Titanic, este post te va interesar. Describiré los gráficos y las técnicas que utilicé para la descripición inicial. Este proyecto es parte de la competencia de kaggle.com disponible en https://github.com/progamandoconro/Titanic-survivors. Espero sea de tu agrado.

Primero que nada, mis respetos a los familiares y víctimas de la catástrofe. El objetivo de este estudio es generar predicciones de la sobrevivencia de los pasajeros en función ciertas variables sociales y económicas. La evaluación es realizada por kaggle.com a partir de un set de datos nuevo, desconocido para los competidores.

Para este post trabajé con Python 3.7 instalado con Anaconda en un Sistema Operativo Linux con Debian 9 (Deepin).

Cargamos las librerías de Python que vamos utilizar y creamos una función para limpiar la consola de trabajo con cls(). La data puedes encontrarla en: https://www.kaggle.com/c/titanic/data.

import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

#clean console
def cls():
    os.system('cls' if os.name=='nt' else 'clear')

rain = pd.read_csv('~/PycharmProjects/titanic/train.csv')

test = pd.read_csv('~/PycharmProjects/titanic/test.csv')

Leímos la data con la librería pandas.  Ahora vamos a hacer unas exploraciones generales para entender mejor la data y descifrar patrones interesantes que podemos usar más adelante en nuestras predicciones.

Las 5 primeras filas de  la data lucen así:

train.head()

Ahora hagamos unas tablas de frecuencias:

survived_sex = pd.crosstab(index=train["Survived"],
                           columns=train["Sex"])
survived_class = pd.crosstab(index=train["Survived"],
                           columns=train["Pclass"])
survived_em = pd.crosstab(index=train["Survived"],
                           columns=train["Embarked"])

survived_class.index= ["died","survived"]
survived_sex.index= ["died","survived"]
survived_em.index= ["died","survived"]

survived_class
survived_sex
survived_em

Parece que la película  tiene razón hasta ahora, los ricos (1 clase) tuvieron más probabilidad de sobrevivir. El sexo también es una variable importante que determinó la sobrevivencia, «mujeres y niños primero». La zona de embarque, sobretodo en la S, también jugó un papel importante en la sobrevivencia de los pasajeros.

Ahora veamos combinaciones de interés para encontrar patrones:

sns.catplot(y="Age", x="Survived", hue="Pclass", data=train,kind="swarm");
plt.show()

Trágico, los niños pobres (3ra clase), representan la vasta mayoría de niños que murieron en la catástrofe. Había pocos niños en primera clase (4), sólo uno murió. Los niños de segunda clase (17) , afortunadamente sobrevivieron todos. En cuanto a las personas mayores, los más ancianos murieron, salvo el más anciano de todos, que viajaba en primera. En la sección de adultos (entre 20 y 50 años), vemos que se cumple lo visto en la película. Los pobres representaron la mayor proporción de los muertos entre las clases.

sns.catplot(y="Age", x="Survived", hue="Sex", data=train,kind="swarm");
plt.show()

Podemos ver en azul al sexo masculino y en naranja al sexo femenino. El primer grupo a la izquierda son los fallecidos, a la derecha los sobrevivientes. En el eje y tenemos la edad del pasajero.

En conclusión de la exploración inicial,  los sobrevivientes presentaron una proporción de clases y sexo distinta a cuando abordaron el barco (Leonardo se enteró). La distribución de sexos, edad y clase económica entre los fallecidos y los sobreviviente fue muy desigual.

Veamos ahora otras variables menos intuitivas de la data.

survived_sib = pd.crosstab(index=train["Survived"],
                           columns=train["SibSp"])
survived_parch = pd.crosstab(index=train["Survived"],
                           columns=train["Parch"])

survived_sib.index= ["died","survived"]
survived_parch.index= ["died","survived"]

survived_sib
survived_parch

El parch se refiere al número de familia sanguínea a bordo.  Las personas con 1, 2 o 3 hijos o padres a bordo fueron los más afortunados.  Con 0 familia, sólo un tercio sobrevivió. El caso de personas con un hijo o padre a bordo fue el único dónde el número de sobrevivientes superó al de fallecidos, siendo poco más de la mitad de sobreviventes. En el caso de la variable SibSp (esposa o familiar no sanguíneo), se cumplen relaciones parecidas a las descritas para el Parch, teniendo las personas con  1 SibSp mayor proporción de sobrevivientes. Vamos graficar a estas dos variables junto con la sobrevivencia para visualizar mejor que sucedió.

sns.catplot(y="SibSp", x="Parch", hue="Survived", data=train,kind="swarm")
plt.show()

Efectivamente, podemos ver secciones dónde hubo mayor sobrevivencia. Entre los parches 1 y dos, con SibSp también cercana a 1 y 2, se presentaron las mayores proporciones de sobrevivientes. Esto nos hace pensar que las familias pequeñas tuvieron mejor chance de sobrevivir que las grandes o las personas solas.

Todos estos patrones que hemos encontrado, vislumbran que también los algoritmos harán un buen trabajo y lograremos predicciones con precisión elevada.

Ahora, un poco de Machine Learning para tratar de predecir la sobrevivencia basados en la información que tenemos. Esta vez en R y Python en conjunto gracias a la librería «reticulate». El boosting lo haremos al final enteramente en Python y, como proyecto futuro, una app dónde, según la data socio-económica que el usuario suministre, se estime la probabilidad de sobrevivencia hipotética en caso de haber estado en el naufragio.

A continuación la minería de la data suministrada en la competencia, con el fin de generar un set preparado para el entrenamiento de los algoritmos. Ahora en R.

setwd("~/PycharmProjects/titanic")

library(caret)
library(randomForest)
library(MLmetrics)
library(mlbench)
library(zoo)
library(ROSE)
library(reticulate)
library(MASS)
library(qdap)
library(stringr)

normalize < - function(x)
{
return ((x - min(x)) / (max(x) - min(x)))
}

t < - read.csv("train.csv")
target < - t$Survived

colnames(t)
exp < - t[, c(3:8, 10:12)]
colnames(exp)

x < - data.frame(target, exp)
# x<-x[sapply(x, function(x) !any(is.na(x)))]
# vamos a atender el asunto de los NA

str(x)

# tenemos tres variables con factores, vamos a cambiarlas a numericas

x$Sex < -ifelse(x$Sex == "female", 0, 1)

table(x$Embarked)
empty < - x$Embarked == ""
x$Embarked[empty] < - "S"
x$Embarked < -ifelse(x$Embarked == "C", 1, ifelse(x$Embarked == "Q", 2, 3))

CABIN < - lapply(x$Cabin, substring, 1, 1)
x$Cabin < -unlist(CABIN)
table(x$Cabin)
x$Cabin < -as.factor(x$Cabin)
table(x$Cabin)
x$Cabin < -ifelse(x$Cabin == "A",
1, ifelse(x$Cabin == "B",
            2, ifelse(x$Cabin == "C",
                        3, ifelse(x$Cabin == "D",
                                    4, ifelse(x$Cabin == "E",
                                                5, ifelse(x$Cabin == "F",
                                                            6, ifelse(x$Cabin == "G",
                                                                        7, ifelse(x$Cabin == "T", 8, 0))))))))

# table(x$Cabin)

# x$Age==ifelse(x$Age<15,1,ifelse(x$Age<30,2,ifelse(x$Age<45,3,4)))

# x$SibSp2=ifelse(x$SibSp+x$Parch<1,1,ifelse(x$SibSp+x$Parch<4,2,ifelse(x$SibSp+x$Parch<6,3,4)))



nombres = t$Name % > %
str_replace_all(c(".*Mrs.*" = "MRS.", 
              ".*Mr.*" = "MR.", ".*Dr.*" = "DR.",
                  ".*Capt.*" = "CAPT", ".*Col.*" = "COL",
                  ".*Miss.*" = "MISS", ".*Rev.*" = "REV",
                  ".*Master.*" = "MASTER", ".*Major.*" = "MAJOR",   
                  ".*Countess.*" = "COUNTESS",
                  ".*Jonkheer.*" = "JONKHEER", 
                  ".*Don.*" = "DON", ".*Mlle.*" = "MLLE",
                  ".*Ms.*" = "MS", ".*Mme.*" = "MME"))


table(nombres)

x$Name = ifelse(nombres == "MR.", 1, ifelse(nombres == "MRS.", 2, ifelse(nombres == "MISS", 3, 4)))

x < -x
colnames(x)

x < -na.aggregate(x)

str(x)

x < -lapply(x, normalize)

x < -as.data.frame(x)

# TT<- lapply(t$Ticket, substring, 1, 1)
# x$Ticket<-as.factor(unlist(TT))
pos = x[x$target == 1,];neg = x[x$target == 0,]
set.seed(7)
epos = sample(1:nrow(pos), floor(0.7 * nrow(pos)))
post = pos[epos,]
poscv = pos[-epos,]

set.seed(7)
eneg = sample(1:nrow(neg), floor(0.7 * nrow(neg)))
negt = neg[eneg,]
negcv = neg[-eneg,]
train = rbind(post, negt);
cv = rbind(poscv, negcv)

A continuación enlacemos R y Python para obtener un algoritmo más potente presente en librerías Python. Afortunadamente, es pósible directamente desde R lo que facilita la creación del código.

Necesitas instalar la librería scikit-learn, por ejemplo, con

conda install scikit-learn . En mi caso, con Anaconda 3.7.

# py_config
# use_python("/home/ro/anaconda3/envs/r-tensorflow/bin/python")

k =
import

("sklearn.ensemble")

rf = k$RandomForestClassifier()

apply(train, 2, function(x)
any( is.na(x)))

rf$fit(train[, -1], train[, 1])

k$RandomForestClassifier(bootstrap=TRUE, class_weight=NULL, criterion='gini',
max_depth = NULL, max_features = 'auto', max_leaf_nodes = NULL,
min_impurity_split = 1e-07, min_samples_leaf = 1,
min_samples_split = 2, min_weight_fraction_leaf = 0.0,
n_estimators = 10, n_jobs = 1, oob_score = FALSE, random_state = NULL,
verbose = 0, warm_start = FALSE)
y_pred = rf$predict(cv[, -1])

roc = roc.curve(y_pred, cv[, 1])

roc$auc

#############


k2 =
import

("sklearn.ensemble")

rf2 = k2$RandomForestClassifier()

rf2$fit(train[, -1], train[, 1])

k2$RandomForestClassifier(bootstrap=TRUE, class_weight=NULL, criterion='gini',
max_depth = NULL, max_features = 'auto', max_leaf_nodes = NULL,
min_impurity_split = 1e-07, min_samples_leaf = 1,
min_samples_split = 2, min_weight_fraction_leaf = 0.0,
n_estimators = 10, n_jobs = 1, oob_score = FALSE, random_state = NULL,
verbose = 0, warm_start = FALSE)
y_pred = rf2$predict(cv[, -1])
#######
k3 =
import

("sklearn.ensemble")

rf3 = k3$RandomForestClassifier()

rf3$fit(train[, -1], train[, 1])

k3$RandomForestClassifier(bootstrap=TRUE, class_weight=NULL, criterion='gini',
max_depth = NULL, max_features = 'auto', max_leaf_nodes = NULL,
min_impurity_split = 1e-07, min_samples_leaf = 1,
min_samples_split = 2, min_weight_fraction_leaf = 0.0,
n_estimators = 10, n_jobs = 1, oob_score = FALSE, random_state = NULL,
verbose = 0, warm_start = FALSE)
y_pred = rf3$predict(cv[, -1])

Hemos entrenado 3 algoritmos para obtener un promedio de sus predicciones. Leamos la data test proporcionada en la competición y realicemos las mismas operaciones que aplicamos a la data de entrenamiento.

test < -read.csv("test.csv")
colnames(test)
colnames(x)
exp < - test[, c(-1, -8)]
colnames(exp)

x < - exp

str(x)
x$Sex < -ifelse(x$Sex == "female", 0, 1)

table(x$Embarked)
empty < - x$Embarked == ""
x$Embarked[empty] < - "S"
x$Embarked < -ifelse(x$Embarked == "C", 
1, ifelse(x$Embarked == "Q", 2, 3))

CABIN < - lapply(x$Cabin, substring, 1, 1)
x$Cabin < -unlist(CABIN)
table(x$Cabin)
x$Cabin < -as.factor(x$Cabin)
table(x$Cabin)
x$Cabin < -ifelse(x$Cabin == "A",
1, ifelse(x$Cabin == "B",
              2, ifelse(x$Cabin == "C",
              3, ifelse(x$Cabin == "D",
              4, ifelse(x$Cabin == "E",
              5, ifelse(x$Cabin == "F",
              6, ifelse(x$Cabin == "G",
                                                                        7, ifelse(x$Cabin == "T", 8, 0))))))))

table(x$Cabin)

nombres = test$Name % > %
str_replace_all(c(".*Mrs.*" = "MRS.",
                   ".*Mr.*" = "MR.",
                   ".*Dr.*" = "DR.",
                   ".*Capt.*" = "CAPT", ".*Col.*" = "COL",
                   ".*Miss.*" = "MISS", ".*Rev.*" = "REV",
                   ".*Master.*" = "MASTER", ".*Major.*" = "MAJOR",
                   ".*Countess.*" = "COUNTESS",
                   ".*Jonkheer.*" = "JONKHEER",
                   ".*Don.*" = "DON",
                   ".*Mlle.*" = "MLLE",
                   ".*Ms.*" = "MS", ".*Mme.*" = "MME"))

table(nombres)

x$Name = ifelse(nombres == "MR.", 1, 
         ifelse(nombres == "MRS.", 2, 
         ifelse(nombres == "MISS", 3, 4)))

colnames(x)

x < -na.aggregate(x)

str(x)

x < -lapply(x, normalize)

x < -as.data.frame(x)

xtest < - x

Ahora realizamos las predicciones con los 3 algoritmos creados previamente. Los tres algoritmos entrenados han funcionado cómo en la película Minority Report. Calculamos un promedio entre ellos y transformamos el resultado a binario nuevamente. En este punto se podría intentar ajustar el umbral de decisión, aunque 0.5 es lo acostumbrado.

Survived1 < -rf$predict(xtest) 
Survived2 < -rf2$predict(xtest) 
Survived3 < -rf3$predict(xtest) 

Survived < - (Survived1 + Survived2 + Survived3) / 3 

Survived < -ifelse(Survived < 0.5, 0, 1)

 Unos detalles más y tenemos nuestra data para ser enviada a ser evaluada en kaggle.

PassengerId < -test$PassengerId 

submi < -data.frame(PassengerId, Survived) 

View(submi) 

write.csv(submi, "SubmissionC6.csv", row.names = FALSE)

 Esto ha permitido una mejoría en las predicciones, hasta colocarnos en el centro de la tabla de posiciones, con un 78% de acierto.

Accuracy = 0.77990 

Ahora es tiempo de realizar unas revisiones de nuestras predicciones y modificar los valores que, por el conocimiento que tenemos de la data, no se ajustan a lo que pensamos esperar.

Modificadas las predicciones por el criterio de información que obtuvimos en el análisis inicial, sumamos nuestro conocimiento e intuición a la matriz de predicciones. El resultado: nos colocamos en el top 15% en el rating de posiciones de la competencia.

Accuracy = 0.79904 

(Top 15% con predicciones para el día 13 de Mayo)

He escrito otro artículo  para describir como mejorar estas predicciones utilizando ensamblaje de algoritmos Machine Learning. Te dejo el link para acceder. Allí se logra mejorar el código para alcanzar el 8% en el top de la tabla de posiciones.

Accuracy = 0.800861

TOP 8% Usuario DataScienceconRo

https://www.kaggle.com/c/titanic/leaderboard

(Predicciones para el día 21 de Mayo)

Hemos alcanzado colocarnos en el top 8% de la competición. Te invito a escribirime si tienes dudas o sugerencias. También puedes descargar el código actualizado en mi repo de github.

Hasta pronto. Gracias por leer.