Arturo S. Lewinger

Objetivo

El Titanic fue, para su época, el mayor transatlantico a nivel mundial. En su viaje inaugural chocó contra un iceberg que lo hundió el 15 de abril de 1912. Debido a las pobres normas de seguridad, 1514 de las 2223 personas que viajaban fallecieron en la tragedia. En el presente trabajo se realizarán modelos de regresión logística para clasificar, analizar y comparar las chances de supervivencia de personas con distintas características.

Preprocesamiento de datos

#Cargo los paquetes que utilizaremos en el TP
library(tidyverse)
library(dplyr)
library(broom)
library(ISLR)
library(ggplot2)
library(GGally)
library(modelr)
library(pROC)
library(cowplot)
library(OneR)
library(rlang)
library(readr)
library(caret)
library(funModeling)
set.seed(1234)
#Cargo los datos (dajamos los de test para el final)
titanic_complete_train <- read_csv("C:/Users/Artur/Desktop/Arturo/Maestria-DM/github/EEA2019/TP3/titanic_complete_train.csv")
titanic_complete_test <- read_csv("C:/Users/Artur/Desktop/Arturo/Maestria-DM/github/EEA2019/TP3/titanic_complete_test.csv")
#vemos las 12 variables y una muestra de los primero datos del universo de 891 datos
glimpse(titanic_complete_train)
## Observations: 891
## Variables: 12
## $ PassengerId <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,...
## $ Survived    <dbl> 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0,...
## $ Pclass      <dbl> 3, 1, 3, 1, 3, 3, 1, 3, 3, 2, 3, 1, 3, 3, 3, 2, 3,...
## $ Name        <chr> "Braund, Mr. Owen Harris", "Cumings, Mrs. John Bra...
## $ Sex         <chr> "male", "female", "female", "female", "male", "mal...
## $ Age         <dbl> 22.00000, 38.00000, 26.00000, 35.00000, 35.00000, ...
## $ SibSp       <dbl> 1, 1, 0, 1, 0, 0, 0, 3, 0, 1, 1, 0, 0, 1, 0, 0, 4,...
## $ Parch       <dbl> 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 1, 0, 0, 5, 0, 0, 1,...
## $ Ticket      <chr> "A/5 21171", "PC 17599", "STON/O2. 3101282", "1138...
## $ Fare        <dbl> 7.2500, 71.2833, 7.9250, 53.1000, 8.0500, 8.4583, ...
## $ Cabin       <chr> NA, "C85", NA, "C123", NA, NA, "E46", NA, NA, NA, ...
## $ Embarked    <chr> "S", "C", "S", "S", "S", "Q", "S", "S", "S", "C", ...
#Utilizo la libreria funModeling para revisar si existen valores faltanten, ceros y valores únicos
status(titanic_complete_train)

Como se observa, luego del recorte de variables no quedan registros faltantes en nuestro dataset. Parch y SibSP tienen una alta proporción de ceros lo que indica una mayoría de gente viajando sin familiares.

#Selecciono las variables que nos interesan. Luego cambio el tipo de variable a factor en tres casos.
titanic_complete_train <- titanic_complete_train %>% 
  select(PassengerId, Survived, Pclass, Sex, Age, SibSp,Parch, Fare, Embarked) %>% 
  mutate(Survived=as.factor(Survived),Pclass=as.factor(Pclass),Embarked=as.factor(Embarked))
titanic_complete_train

Detalle de las variables:

ggpairs(titanic_complete_train %>% 
          select(Survived, Pclass, Sex, Age, Fare),
        mapping = aes(colour= Survived),legend = c(1,1)) + 
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
  theme_bw()+
  theme(legend.position = "bottom")

Algunas primeras observaciones del ggpairs:

Adicionalmente, en el primer gráfico (arriba a la izquierda) se observa la distribución de los sobrevivientes. Para ver exactamente la cantidad de personas en cada clase hacemos lo siguiente:

#sumo los casos que sobrevivieron y que no en el dataset de training
titanic_complete_train %>% 
  group_by(Survived) %>% 
  summarise(numero_casos=n(),porcentaje=numero_casos/nrow(titanic_complete_train)*100)

Vemos que en nuestra dataset de training el 38% de los pasajeros sobrevivió.

Partición en entrenamiento y validación

#hacemos el corte del dataset en Train y test 70% y 30%
set.seed(19450)
train_test <- titanic_complete_train %>% 
  resample_partition(c(train=0.7,test=0.3))

titanic_train <- train_test$train %>% as_tibble()
titanic_validacion <- train_test$test %>% as_tibble()
#Chequeamos si quedaton balanceadas

rbind(titanic_train %>% 
  group_by(Survived) %>% 
  summarise(dataset="train",numero_casos=n(),porcentaje=numero_casos/nrow(titanic_train)*100),
  titanic_validacion %>% 
  group_by(Survived) %>% 
  summarise(dataset="validacion",numero_casos=n(),porcentaje=numero_casos/nrow(titanic_validacion)*100)
)

Apenas cambio del 38,4% al 38,5% en el “train” y 38,1% en la “validación”. Habiendo confirmado que se encuentran balanceados, procedenmos a los entrenamientos para realizar predicciones.

Primer modelo de predicción

Modelo con variables “Pclass”, “Sex” y “Age”

Buscamos calcular la probabilidad \(P(Sobrevivir=Yes|X)=P(X)\) , es decir, la probabilidad de sobrevivir dado que conocemos en que clase viajaba, su sexo y edad. Luego elegiremos un punto de corte que nos separe lo mejor posible los que sobreviven y los que no.

#realizamos el modelo de regresión logistica con las variables mencionadas
modelo1 <- glm(Survived ~ Pclass + Sex + Age, family = 'binomial', data = titanic_train) # GLM se refiere a un modelo lineal generalizado, que al especificarlo como "binomial", calcula la regresión logistica.

#Hacemos una salida para revisar los coeficientes y su significatividad
coef_modelo1 <- tidy(modelo1)
coef_modelo1

Del análisis de los coeficientes se obtiene que:

  • Tanto el intercepto como el resto de los coeficientes resultaron ser significativos de acuerdo al P-value.

  • Las clases tomadas como base (es decir, que son absorvidas por el intercepto) son: “1ra Clase” en cuanto a Pclass y “Mujer” en Sex.

  • Pclass : para la 2da, el coeficiente es de -1,4 lo que indica que si en lugar de viajar en primera clase se viaja en 2da las chances de sobrevivir disminuyen. En 3ra, pasa lo mismo pero las chances de no sobrevivir respecto a primera clase disminuyen aún mas.

-Sex: Coeficiente para “Male” -2,8. Con lo cual podemos afirmar que los hombres tienen menos probabilidades de sobrevivir que las mujeres (categoria basal elegida por R)

-age: Coeficiente de -0,04, es decir, que a medida que se incrementa las probabilidades de sobrevivir disminuyen ligeramente.

Predicción puntual

Para saber el caso puntual de quien tiene mayor probabilidad de sobrevivir, si Rose que es una mujer de 17 años que viaja en primera clase o Jack que es un hombre de 20 años viajando en tercera clase, evaluamos estos casos con nuestro modelo.

#armo un dataframe con las caracteristicas de las dos personas a evaluar. La primera es Rose y la segunda Jack
set.seed(1234)
rose_jack <- data.frame(
  Pclass = as.factor(c(1, 3)),
  Sex = c('female', 'male'),
  Age = c(17.0, 20.0))

# Ahora evaluamos la predicción para los dos casos, pidiendo probabilidad en la respuesta

prediccion <- predict(modelo1, newdata = rose_jack, type = 'response')
prediccion
##         1         2 
## 0.9704455 0.1107944

Las probabilidades de sobrevivir de Rose son sumamente altas (97%) mientras que las de Jack, son muy bajas (11%)

Nuevos modelos combinando variables

Para realizar combinaciones de variables en varios modelos será conveniente usar la función “formulas” del paquete “modelr”.

# Armamos varios modelos a ser evaluados en un solo script usando dicha función.

logit_formula <- formulas(.response = ~Survived, #Variable a predecir que siempre es la misma
                         modelo1_class_sex_age = ~Pclass+Sex+Age, #Modelo con las mismas variables que ya usamos anteriormente 
                         m2_pclass_sex= ~Pclass+Sex, #Sacando solo Age
                         m3_pclass_age = ~Pclass+Age, #Sacando solo Sex
                         m4_pclass = ~ Pclass ) #Solo con Pclass

Ahora armamos una tabla donde aplicamos los distintos modelos al dataset de train

modelos_train <- tibble(logit_formula) %>% # armamos un dataframe con todo
  mutate( modelName = names(logit_formula), #Agregamos una columna para saber el nombre de modelo utilizado
          expression = paste(logit_formula),# Columna con el detalle de variables usadas
          mod = map(logit_formula, ~glm(.,family = 'binomial', data = titanic_train)))# Se aplica cada formula al dataset de entrenamiento

Evaluación de modelos

Calculo ahora, para cada modelo, las medidas de evaluación y las ordeno

modelos_train <- modelos_train %>% 
  mutate(glance = map(mod,glance))

# Obtenemos las medidas de evaluación de interés
modelos_train %>%
  unnest(glance) %>%
  mutate(percExplainedDev = 1-deviance/null.deviance) %>% # Calculo el porcentaje de deviance explicada
  select(modelName, percExplainedDev, null.deviance, deviance, logLik, expression) %>% # Selecciono sólo algunas variables de interés
  arrange(deviance) # Ordeno de mayor a menor por deviance

El modelo con menor deviance es el que utiliza mas variables (modelo 1). A medida que se disminuyen variables vemos que el deviance aumenta.

AUC y curva ROC

Para graficar la curva ROC del modelo 1 con las 3 variables, primero es necesario calcular sus predicciones.

# Añadimos una columna con las predicciones a los datos
modelos_train <- modelos_train %>%
  mutate(pred = map(mod, augment, type.predict = "response"))

# Extraemos las predicciones de nuestro modelo 1
prediction_m1 <-
  modelos_train %>%
  filter(modelName == "modelo1_class_sex_age") %>%
  unnest(pred, .drop = TRUE)

Y con los datos obtenidos previamente, graficamos.

# Armo la curva ROC
roc_m1 <-
  roc(response = prediction_m1$Survived, predictor = prediction_m1$.fitted)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
# Grafico
ggroc(roc_m1, size = 1.2,color="blue") +
  geom_abline(slope = 1,
              intercept = 1,
              linetype = 'dashed',
              color="red") +
  
  theme_bw() +
  labs(
    title = 'Curva ROC',
    caption = paste('AUC:', round(roc_m1$auc, 3)),
    subtitle = "Modelo 1 usando 3 variables"
  )

En el eje Y, se observa la “sensitivity”, que representa la proporción de sobrevivientes sobre el total de casos. Es decir que la curva ofrece una vista a los casos que el algoritmo estimó como positivos y se condice con lo realmente ocurrido a medida que se amplia el punto de corte en la probabilidad estimada de supervivencia de cada sujeto. Por su parte en el eje X tenemos la “especificidad”, que es la proporción de casos que no sobrevivieron sobre el total cuando movemos dicho punto de corte (acá de derecha a izquierda).

En otras palabras, si ordenamos de mayor a menor a las personas con su probabilidad de supervivencia, lo que vemos en el gráfico (de izquierda a derecha en el eje x) es que si tomamos como punto de corte, por ejemplo, a la probabilidad en la persona que divide el 50% de los casos que sobrevivieron (eje Y), vemos que la curva azul apenas avanza hacia la derecha sino que mayormente crece sobre el eje Y. Esto significa que con ese punto de corte, los True Positive (verdaderos sobrevivientes) predichos son muchos en relación a los que no sobrevivieron. A partir de allí, para incrementar la cantidad de sobrevivientes reales nuestro modelo cometerá mas errores, es decir, empieza a incluir los casos que no sobrevivieron. La línea punteada roja representa una aproximación al comportamiento de predecir de forma aleatoria. Cuanto mayor sea el área que deja la curva ROC por sobre la linea roja punteada, mejor será la predicción. Dicho área se conoce como AUC (Area under de curve) y en este caso es de 0,874. Cabe mencionar, que no existe un mejor punto de corte, sino que este depende del problema al que nos enfrentemos y que tipo de error estemos dispuestos a cometer, ya que existe un “trade-off” a lo largo de la curva.

Gráfico violín

ggplot(prediction_m1, aes(x = Survived, y = .fitted, group = Survived, fill = factor(Survived))) + 
  geom_violin() +
  theme_bw() +
  guides(fill=FALSE) +
  labs(title='Gráfico violín', subtitle='Modelo 1 con tres variables', y='Probabilidad predicha')

La forma de leer este gráfico es la siguiente: Si ordenamos de mayor a menor la probabilidad predicha por nuestro modelo de sobrevivir de los sujetos, podríamos elegir cualquier punto de corte a lo largo del eje Y. La bondad del ajuste del modelo 1 se observa en que si eligieramos por ejemplo un punto con probabilidad 0,50 y trazamos una linea imaginaria horizontal, vemos que las personas que no sobrevivieron (en rojo en nuestro gráfico) quedarán (en mayoria) situadas por debajo de esta, mientras que a su vez la mayor proporción de los que sobrevivieron quedarán por encima. Cuanto mejor divide esta linea los positivos de los negativos, mejor será el modelo.

Elección del punto de corte

Hasta aquí evaluamos el modelo de forma general, pero si queremos avanzar en establecer un punto de corte concreto, se deberá utilizar las predicciones realizadas sobre validación. Para encontrar dicho punto se graficarán (en forma de curvas) los resultados de hacer distintos cortes a lo largo de todas las probabilidades predichas. Como resultados, agregaremos a la sensitividad (o recall) y especificidad, ya abordados, el “acurracy” (proporción de casos clasificados correctamente sobre el total) y “precision” (proporción de positivos clasificados correctamente encontrados entre el total de estos).

Entonces, primero calculamos las predicciones para la validación.

# Calculo predicciones al dataset de validación
modelo_validacion <- modelos_train %>%
  filter(modelName == 'modelo1_class_sex_age') %>%
  mutate(pred = map(mod, augment, newdata = titanic_validacion, type.predict = "response"))

# Obtengo las correspondientes al mejor modelo
prediction_validation_class_sex_age <-
  modelo_validacion %>%
  filter(modelName == "modelo1_class_sex_age") %>%
  unnest(pred, .drop = TRUE)
# Hago una función que clasifica las observaciones para diferentes puntos de corte y devuelve las métricas de performance vistas
performance_metrics <- function(cutoff, predictions = prediction_validation_class_sex_age) {
    table_res <- predictions %>%
      mutate(
        predicted_class = if_else(.fitted > cutoff, 1, 0) %>% as.factor(),
        Survived = factor(Survived)
      )
    
    confusionMatrix(table_res$predicted_class, table_res$Survived, positive = "1") %>%
      tidy() %>%
      select(term, estimate) %>%
      filter(term %in% c('accuracy', 'sensitivity', 'specificity',
                         'precision', 'f1')) %>%
      mutate(cutoff = cutoff)
  }

# Hago el vector con todos los cortes en intervalos de 0.005
cutoffs = seq(
  min(prediction_validation_class_sex_age$.fitted),
  max(prediction_validation_class_sex_age$.fitted),
  0.005
)

# Aplico la función al vector
logit_pred = map_dfr(cutoffs, performance_metrics) %>% mutate(term = as.factor(term))

# Grafico las métricas
logit_pred %>%
  ggplot(., aes(cutoff, estimate, group = term, color = term)) +
  geom_line(size = 1.5) +
  theme_minimal() +
  labs(title = 'Accuracy, Sensitivity, Specificity y Precision',
       subtitle = 'Modelo usando Pclass, Sex y Age',
       color = "") + 
  # Escalo ticks de los ejes
  scale_x_continuous(breaks = round(seq(0, 1, by = 0.1), 2)) +
  scale_y_continuous(breaks = round(seq(0, 1, by = 0.1), 2))+
  theme(legend.position = "bottom")
## Warning: Removed 3 rows containing missing values (geom_path).

En el gráfico anterior se ve claramente el trade-off entre “sensitividad” y “especificidad”, es decir, que a medida que subimos la probabilidad de corte la primera disminuye mientras que la segunda aumenta. A diferencia de por ejemplo predicciones de enfermedades graves donde uno preferiría pecar de excesivamente cauteloso para detectar la mayor cantidad de casos positivos incrementando los falsos positivos, en nuestro caso no tenemos algún tipo de preferencia particular. Por ello, trataremos de hacer un equilibrio para detectar los verdaderos positivos y negativos, con lo cual elegimos el punto donde “sensitivy” se une con “specificity”. Este punto da como corte un valor cercano a 0,46 y se puede observar que el “accuracy” esta cercano a su punto máximo.

Vemos ese punto en el siguiente gráfico de violín y comprobamos que efectivamente parece separar bien la mayoria de casos positivos de los negativos.

ggplot(prediction_validation_class_sex_age, aes(x = Survived, y = .fitted, group = Survived, fill = factor(Survived))) + 
  geom_violin() +
  geom_abline(slope = 0, intercept = 0.46, color = 'blue', linetype = 'dashed', size = 0.5) +
  theme_bw() +
  guides(fill=FALSE) +
  labs(title='Gráfico Violín', subtitle='Modelo seleccionado con datos de validación', y='Probabilidad predicha')

Matriz de confunsión para el punto de corte seleccionado

set.seed(1234)
#hago el mismo modelo (entrenado con datos de train)
modelo <-
  glm(Survived ~ Pclass + Sex + Age, family = 'binomial', data = titanic_train)
# Aplico modelo de train en los nuevos datas de validación y calculo las nuevas predicciones
table = augment(x = modelo,
                newdata = titanic_validacion,
                type.predict = 'response')
# Clasifico utilizando el punto de corte
table = table %>% mutate(
  predicted_class = if_else(.fitted >= 0.46, 1, 0) %>% as.factor(),
  Survived = factor(Survived)
)
# Creo la matriz de confusión
confusionMatrix(table$predicted_class, table$Survived, positive = "1")
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 131  22
##          1  35  80
##                                           
##                Accuracy : 0.7873          
##                  95% CI : (0.7334, 0.8347)
##     No Information Rate : 0.6194          
##     P-Value [Acc > NIR] : 2.775e-09       
##                                           
##                   Kappa : 0.5597          
##                                           
##  Mcnemar's Test P-Value : 0.112           
##                                           
##             Sensitivity : 0.7843          
##             Specificity : 0.7892          
##          Pos Pred Value : 0.6957          
##          Neg Pred Value : 0.8562          
##              Prevalence : 0.3806          
##          Detection Rate : 0.2985          
##    Detection Prevalence : 0.4291          
##       Balanced Accuracy : 0.7867          
##                                           
##        'Positive' Class : 1               
##

La matriz de confusión (por como fue armada) presenta en las filas los valores predichos y en las columnas los valores reales. Èsta predice correctamente 131 personas que no sobrevivieron y 80 que sí lo hicieron. El Accuracy máximo deberia estar cercano a las 211 casos correctos que fueron predichos correctamente con este corte. El accuracy es de 0,78 y poblacionalmente se encontraría, con un intervalo de confianza del 95%, entre 0,73 y 0,83.

Testing

Por primera vez leemos los datos de test y le realizamos las mismas transformación hechas al dataset de entrenamiento.

titanic_complete_test <- titanic_complete_test %>% 
  select(PassengerId, Survived, Pclass, Sex,
         Age, SibSp, Parch, Fare, Embarked) %>% #selecciono las variables
  mutate_at(vars(Survived, Pclass, Embarked), as_factor) #transformo en factor a estas variables
set.seed(1234)
# El modelo no hace falta volver a calcularlo porque será el entrenado con train (lo llamamos "modelo")
# Obtenemos las predicciones para el test con dicho modelo
pred_test = augment(x = modelo, newdata=titanic_complete_test,type.predict = 'response')

# Clasifico utilizando el punto de corte
pred_test = pred_test %>% 
  mutate(predicted_class = if_else(.fitted >= 0.46, 1, 0) %>% as.factor(),
  Survived = factor(Survived))

confusionMatrix(pred_test$predicted_class, pred_test$Survived, positive = "1")
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 203  46
##          1  58 111
##                                           
##                Accuracy : 0.7512          
##                  95% CI : (0.7069, 0.7919)
##     No Information Rate : 0.6244          
##     P-Value [Acc > NIR] : 2.419e-08       
##                                           
##                   Kappa : 0.4775          
##                                           
##  Mcnemar's Test P-Value : 0.2807          
##                                           
##             Sensitivity : 0.7070          
##             Specificity : 0.7778          
##          Pos Pred Value : 0.6568          
##          Neg Pred Value : 0.8153          
##              Prevalence : 0.3756          
##          Detection Rate : 0.2656          
##    Detection Prevalence : 0.4043          
##       Balanced Accuracy : 0.7424          
##                                           
##        'Positive' Class : 1               
##

Como puede verse de las últimas dos salidas, “Sensitivity”, se redujo de 0.78 a 0.71 desde el test a la validación mientras que “Specificity” se redujo apenas 1%. El accuracy del test se redujo de 0,78 a 0,75 respecto al de validación. Es decir que el modelo predijo correctamente el 75% de los casos presentes en test (sumando positivos y negativos).