TP 3: Regresión Logística

1. Preparación de los datos

a. Leer el archivo titanic_complete_train.csv y mostrar su estructura

## Parsed with column specification:
## cols(
##   PassengerId = col_double(),
##   Survived = col_double(),
##   Pclass = col_double(),
##   Name = col_character(),
##   Sex = col_character(),
##   Age = col_double(),
##   SibSp = col_double(),
##   Parch = col_double(),
##   Ticket = col_character(),
##   Fare = col_double(),
##   Cabin = col_character(),
##   Embarked = col_character()
## )
## 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", ...

b. Seleccionar las variables PassengerId, Survived, Pclass, Sex, Age, SibSp,Parch, Fare y Embarked

d y e. Realizar un gráfico de ggpairs para las variables Survived, Pclass, Sex, Age y Fare e interpretarlo. Mostrar la distribución de clase (Sobrevivientes vs No Sobrevivientes)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

A partir del gráfico realizado se pueden extraer, por ejemplo, las siguientes conclusiones:

  • Es mayor la cantidad de personas que subrevivió en comparación a la que falleció en el Titanic.
  • Una proporción elevada de los pasajeron que murieron viajaban en tercera clase.
  • Las personas que sobrevivieron pagaron en general tarifas más altas que las personas que murieron (relacionado al punto anterior)
  • Murieron proporcionalmente más hombres que mujeres.

e. Dividir al dataset en conjunto de entrenamiento (70% de los datos) y validación (30% de los datos). Volver a analizar la distribución de clase para chequear que sea aproximadamente igual entre ambos conjuntos y respecto a la distribución de clase que obtuvieron para todo el dataset en el punto 1)e)

# Divido el dataset
Particion <- createDataPartition(df_training_titanic$Survived, p = .7, 
                                  list = FALSE, 
                                  times = 1)

entrenamiento <- df_training_titanic[ Particion,]  
validacion  <- df_training_titanic[-Particion,] 

# Revisar distribución en los distintos conjuntos
distribucion_training <- df_training_titanic %>%
  group_by(Survived) %>%
  summarise(numero_casos=n())
distribucion_entrenamiento <- entrenamiento %>%
  group_by(Survived) %>%
  summarise(numero_casos=n())
distribucion_validacion <- validacion %>%
  group_by(Survived) %>%
  summarise(numero_casos=n())

a <- ggplot(distribucion_training, aes(x=Survived, y=numero_casos, fill = Survived)) + geom_col() + labs(subtitle = "Dataset completo de training", x = "Clasificación supervivencia", y = "cantidad de pasajeros") + guides(fill = "none") + theme_bw()
b <- ggplot(distribucion_entrenamiento, aes(x=Survived, y=numero_casos, fill = Survived)) + geom_col() + labs(subtitle = "Dataset Entrenamiento", x = "Clasificación  supervivencia", y = "cantidad de pasajeros") + guides(fill = "none") + theme_bw()
c <- ggplot(distribucion_validacion, aes(x=Survived, y=numero_casos, fill = Survived)) + geom_col() + labs(subtitle = "Dataset Validación", x = "Clasificación supervivencia", y = "cantidad de pasajeros") + guides(fill = "none") + theme_bw()

grid.arrange(a,b,c,
  nrow = 1,
  top = "Distribución de los datasets"
  #bottom = textGrob(
  #  "En rojo se presentaron las personas que fallecieron",
  #  gp = gpar(fontface = 4, fontsize = 10),
  #  hjust = 1,
  #  x = 1
  #)
)

2. Predicciones

a. Realizar un modelo de regresión logística para predecir la supervivencia en función de Pclass, Sex y Age. Usar solo el dataset de entrenamiento

## 
## Call:
## glm(formula = Survived ~ Pclass + Sex + Age, family = "binomial", 
##     data = entrenamiento)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.8437  -0.6337  -0.4191   0.6072   2.5251  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.112277   0.490574   8.383  < 2e-16 ***
## Pclass2     -1.270566   0.319276  -3.980 6.91e-05 ***
## Pclass3     -2.566463   0.320177  -8.016 1.09e-15 ***
## Sexmale     -2.785316   0.230051 -12.107  < 2e-16 ***
## Age         -0.043327   0.009746  -4.445 8.77e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 832.49  on 624  degrees of freedom
## Residual deviance: 539.89  on 620  degrees of freedom
## AIC: 549.89
## 
## Number of Fisher Scoring iterations: 5

b. Dar una breve interpretación de los coeficientes y su significatividad

En primera instancia es importante dar cuenta que todos los coeficientes del modelo realizaro resultaron significativos para la predicción de la supervivencia en funcion de la clase del pasajero, sy sexo y su edad. Tanto Pclass2 como Pclass3 presentaron coeficientes negativos. El estimador del coeficiente para Pclass2 fue de -1.04, lo cual puede ser interpretado como que es 104% menos probable que una persona que viajaba en segunda clase sobreviva en comparación con alguien que viajaba en primera clase (categoria de referencia). Del mismo modo, el estimador del coeficiente para Pclass 3 fue de -2.23, lo cual puede ser interpretado como que es 223% más probable que alguien que viajaba en tercera clase muriera en comparación con alguien que viajaba en prmera clase. También uno puede analizar esta predicción en relación al sexo de los pasajeros. El estimador del coeficiente para sexo masculino ha sido de -2.48. Lo cual significaria que es 248% más probable que muera alguien de sexo masculino comparado con alguien de sexo femenino. Finalmente, y en relación a la edad, el estimador del coeficiente para dicha variable es de -0.03. Esto puede ser interpretado como que es 3% más probable que una persona muera por cada año adicional que tiene dicha persona. Es decir que a mayor edad, mayor probabilidad de morir en el Titanic.

c. ¿Quién tiene una mayor probabilidad de supervivencia? 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

Nombre del pasajero Clase Sexo Edad Predicción de supervivencia
Rose 1 female 17 0.9669
Jack 3 male 20 0.1085

Como muestra esta tabla, existe muchas más posibilidades de que Rose sobreviva (94,10% de probabilidad de supervencia) en comparación con Jack (11.55% de probabilidad de supervivencia).

La segunda vez que una tabla le juega a Jack una mala pasada.

#indignado

#indignado

3. Generación de modelos

a. Generar 3 modelos de regresión logística sobre el dataset de entrenamiento utilizando diferentes combinaciones de variables. Al menos dos modelos deben ser multivariados

Modelo 1

El primer modelo intentará predecir la supervivencia en función unicamente de la tarifa pagada por los pasajeros.

## 
## Call:
## glm(formula = Survived ~ Fare, family = "binomial", data = entrenamiento)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7149  -0.8737  -0.8371   1.3604   1.5661  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.005873   0.116364  -8.644  < 2e-16 ***
## Fare         0.017740   0.002864   6.195 5.84e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 832.49  on 624  degrees of freedom
## Residual deviance: 775.24  on 623  degrees of freedom
## AIC: 779.24
## 
## Number of Fisher Scoring iterations: 5

Como se puede ver, la tarifa del pasaje es un predictor de la supervivencia (dicho coeficiente resulto significativo). Por cada unidad monetaria extra pagada por la tarifa, aumenta 1.4% la probabilidad de supervivencia en el Titanic.

Modelo 2

El segundo modelo intentará predecir la supervivencia en función de la tarifa pagada y el sexo del pasajero.

## 
## Call:
## glm(formula = Survived ~ Fare + Sex, family = "binomial", data = entrenamiento)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3582  -0.5886  -0.5597   0.7651   1.9687  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.739424   0.187376   3.946 7.94e-05 ***
## Fare         0.013046   0.003132   4.165 3.11e-05 ***
## Sexmale     -2.614942   0.210947 -12.396  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 832.49  on 624  degrees of freedom
## Residual deviance: 592.97  on 622  degrees of freedom
## AIC: 598.97
## 
## Number of Fisher Scoring iterations: 5

En este segundo modelo realizado, también todos sus coeficientes fueron significativos. Por cada unidad monetaria aumenta 1.03% la probabilidad de supervivencia. Mientras que siendo hombre las probabilidades de sobrevivir son %239.94 inferiores que siendo mujer.

Modelo 3

El tercer modelo intentará predecir la supervivencia en función de la edad, el sexo y la clase del pasaje obtenido.

## 
## Call:
## glm(formula = Survived ~ Age + Sex + Pclass + SibSp, family = "binomial", 
##     data = entrenamiento)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.8681  -0.5651  -0.3964   0.5833   2.5020  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.49335    0.52295   8.592  < 2e-16 ***
## Age         -0.04842    0.01016  -4.765 1.89e-06 ***
## Sexmale     -2.86728    0.23582 -12.159  < 2e-16 ***
## Pclass2     -1.34622    0.32423  -4.152 3.29e-05 ***
## Pclass3     -2.58102    0.32036  -8.057 7.85e-16 ***
## SibSp       -0.29997    0.12357  -2.428   0.0152 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 832.49  on 624  degrees of freedom
## Residual deviance: 533.04  on 619  degrees of freedom
## AIC: 545.04
## 
## Number of Fisher Scoring iterations: 5

Para este último modelo todos los coeficientes también resultaron significativos. Por cada año de edad que tiene extra un individuo, la probabilidad de supervivencia baja 0.95% según este modelo. Del mismo modo, siendo hombre la probabilidad de supervivencia es 262.1% menor que siendo mujer. Además, si uno viaja en segunda o tercera clase tiene mayores probabilidades de morir en el Titanic en comparación con viajar en primera clase. Finalmente, por cada hermano que uno tiene en el Titanic, las probabilidades de morir aumentan %36.05 (según el presente modelo).

Presento a continuación los estimadores para los tres modelos contruidos previamente.

## Warning: `data_frame()` is deprecated, use `tibble()`.
## This warning is displayed once per session.
## # A tibble: 11 x 7
##    models  expression          term   estimate std.error statistic  p.value
##    <chr>   <chr>               <chr>     <dbl>     <dbl>     <dbl>    <dbl>
##  1 modelo1 Survived ~ Fare     (Inte~  -1.01     0.116       -8.64 5.42e-18
##  2 modelo1 Survived ~ Fare     Fare     0.0177   0.00286      6.19 5.84e-10
##  3 modelo2 Survived ~ Fare + ~ (Inte~   0.739    0.187        3.95 7.94e- 5
##  4 modelo2 Survived ~ Fare + ~ Fare     0.0130   0.00313      4.17 3.11e- 5
##  5 modelo2 Survived ~ Fare + ~ Sexma~  -2.61     0.211      -12.4  2.74e-35
##  6 modelo3 Survived ~ Pclass ~ (Inte~   4.49     0.523        8.59 8.52e-18
##  7 modelo3 Survived ~ Pclass ~ Pclas~  -1.35     0.324       -4.15 3.29e- 5
##  8 modelo3 Survived ~ Pclass ~ Pclas~  -2.58     0.320       -8.06 7.85e-16
##  9 modelo3 Survived ~ Pclass ~ Age     -0.0484   0.0102      -4.76 1.89e- 6
## 10 modelo3 Survived ~ Pclass ~ Sexma~  -2.87     0.236      -12.2  5.17e-34
## 11 modelo3 Survived ~ Pclass ~ SibSp   -0.300    0.124       -2.43 1.52e- 2

b. Ordenar por la deviance los 3 modelos creados en el punto 3)a) y el creado en el punto 2)a) y seleccionar el mejor modelo en términos de la deviance explicada

## # A tibble: 4 x 6
##   expression     null.deviance logLik deviance df.residual perc_explained_~
##   <chr>                  <dbl>  <dbl>    <dbl>       <int>            <dbl>
## 1 Survived ~ Pc~          832.  -267.     533.         619           0.360 
## 2 Survived ~ Pc~          832.  -270.     540.         620           0.351 
## 3 Survived ~ Fa~          832.  -296.     593.         622           0.288 
## 4 Survived ~ Fa~          832.  -388.     775.         623           0.0688

De esta forma, podemos ver que el modelo que maximiza la deviance explicada es el modelo 3, que considera las variable de clase, edad, sexo y SibSp". La deviancia explicada por este modelo es de 0.3143.

3. Evaluación del modelo

a. Realizar el gráfico de curva ROC y obtener el AUC para el modelo elegido. Interpretar el gráfico

## Setting levels: control = murió, case = sobrevivió
## Setting direction: controls < cases

## [1] "AUC: Mejor modelo 0.863847402597403"

Para el modelo presentado el AUC es de 0.8445. A partir de esto podemos ver que dicho modelo se aleja del 0.5 de AUC que representaria una predicción como la presente en el azar. Ademas, a partir del análisis de la curva ROC podemos ver que el modelo performa mejor en relación a la especificidad (cantidad de predicciones negativas correctas sobre la cantidad total de eventos negativos reales) en comparación con la sensibilidad (cantidad de predicciones positivas correctas sobre la cantidad total de eventos positivos reales). Es decir que detecta relativamente mejor los verdaderos negativos que los verdaderos positivos.

5. Elección del punto corte

a. Sobre el dataset de validación realizar un gráfico de Accuracy, Specificity, Recall y Precision en función del punto de corte

## # A tibble: 266 x 14
##    models expression .rownames PassengerId Survived Pclass Sex     Age
##    <chr>  <chr>      <chr>           <dbl> <fct>    <fct>  <fct> <dbl>
##  1 model~ Survived ~ 2                   2 sobrevi~ 1      fema~  38  
##  2 model~ Survived ~ 5                   5 murió    3      male   35  
##  3 model~ Survived ~ 9                   9 sobrevi~ 3      fema~  27  
##  4 model~ Survived ~ 22                 22 sobrevi~ 2      male   34  
##  5 model~ Survived ~ 25                 25 murió    3      fema~   8  
##  6 model~ Survived ~ 27                 27 murió    3      male   26.5
##  7 model~ Survived ~ 28                 28 murió    1      male   19  
##  8 model~ Survived ~ 30                 30 murió    3      male   26.5
##  9 model~ Survived ~ 33                 33 sobrevi~ 3      fema~  21.8
## 10 model~ Survived ~ 35                 35 murió    1      male   28  
## # ... with 256 more rows, and 6 more variables: SibSp <dbl>, Parch <dbl>,
## #   Fare <dbl>, Embarked <fct>, .fitted <dbl>, .se.fit <dbl>
## Warning: package 'plyr' was built under R version 3.6.1
## -------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## -------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following object is masked from 'package:purrr':
## 
##     compact
prediction_validation_modelo3$Survived <- revalue(prediction_validation_modelo3$Survived, c("sobrevivió"="1"))
prediction_validation_modelo3$Survived <-
revalue(prediction_validation_modelo3$Survived, c("murió"="0"))

# Función que clasifica las observaciones para diferentes puntos de corte y devuelve métricas de performance asociadas
prediction_metrics <-
  function(cutoff, predictions = prediction_validation_modelo3) {
    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)
  }

# Construyo vector de puntos de corte desde el mínimo al máximo de nuestras probabilidades con intervalo de 0.005
cutoffs = seq(
  min(prediction_validation_modelo3$.fitted),
  max(prediction_validation_modelo3$.fitted),
  0.001
)

# Uso la función con este vector de puntos de corte
logit_pred = map_dfr(cutoffs, prediction_metrics) #%>% mutate(term = as.factor(term))

# Grafico las métricas de performance en función del punto de corte
logit_pred %>%
  ggplot(., aes(cutoff, estimate, group = term, color = term)) +
  geom_line(size = 1.5) +
  theme_minimal() +
  labs(title = 'Accuracy, Sensitivity, Specificity y Precision',
       color = "") + 
  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") +
  geom_vline(xintercept=0.4, linetype="dashed", color = "black")

b. Elegir un punto de corte y explicar su decisión

Se eligió como punto de corte 0.4, teniendo en cuenta que es el punto en donde (aproximadamente) se cruzan sensitivity, specificity y accuracy

c. Obtener la matriz de confusión con el modelo y punto de corte elegidos. Interpretarla

## Confusion Matrix and Statistics
## 
##    
##       0   1
##   0 124  40
##   1  26  76
##                                           
##                Accuracy : 0.7519          
##                  95% CI : (0.6955, 0.8026)
##     No Information Rate : 0.5639          
##     P-Value [Acc > NIR] : 1.461e-10       
##                                           
##                   Kappa : 0.4885          
##                                           
##  Mcnemar's Test P-Value : 0.1096          
##                                           
##             Sensitivity : 0.6552          
##             Specificity : 0.8267          
##          Pos Pred Value : 0.7451          
##          Neg Pred Value : 0.7561          
##              Prevalence : 0.4361          
##          Detection Rate : 0.2857          
##    Detection Prevalence : 0.3835          
##       Balanced Accuracy : 0.7409          
##                                           
##        'Positive' Class : 1               
##

Eligiendo el punto de corte 0.4, es posible observarlos valores para las métricas en el conjunto de validación son (tomando como clase positiva la supervivencia en el Titanic):

  • Accuracy: 0.7519
  • Specificity: 0.8267
  • Recall: 0.7519
  • Precision: 0.7451

Como puede observarse, aquello que mejor se predice son los True Negative (n=124)

6. Dataset de testeo

a. Leer el archivo titanic_complete_test.csv y transformar las variables Survived, Pclass y Embarked a factor

## Parsed with column specification:
## cols(
##   PassengerId = col_double(),
##   Pclass = col_double(),
##   Name = col_character(),
##   Sex = col_character(),
##   Age = col_double(),
##   SibSp = col_double(),
##   Parch = col_double(),
##   Ticket = col_character(),
##   Fare = col_double(),
##   Cabin = col_character(),
##   Embarked = col_character(),
##   Survived = col_double()
## )

b. Con el modelo y punto de corte elegidos clasificar a las personas del dataset de testing.

c. Obtener la matriz de confusión y comparar con la obtenida en el punto 5)c).

## Confusion Matrix and Statistics
## 
##    
##       0   1
##   0 199  62
##   1  43 114
##                                           
##                Accuracy : 0.7488          
##                  95% CI : (0.7044, 0.7897)
##     No Information Rate : 0.5789          
##     P-Value [Acc > NIR] : 3.156e-13       
##                                           
##                   Kappa : 0.4771          
##                                           
##  Mcnemar's Test P-Value : 0.07898         
##                                           
##             Sensitivity : 0.6477          
##             Specificity : 0.8223          
##          Pos Pred Value : 0.7261          
##          Neg Pred Value : 0.7625          
##              Prevalence : 0.4211          
##          Detection Rate : 0.2727          
##    Detection Prevalence : 0.3756          
##       Balanced Accuracy : 0.7350          
##                                           
##        'Positive' Class : 1               
##

Para el conjunto de test las predicciones para todas las metricas son inferiores que para el conjunto de entrenamiento. Sin embargo, los descensos para el conjunto de test no fueron de gran tamaño. Podemos ver que el modelo performa bien en cuanto a la especificidad (82.23%)