Regresion Logistica

Cargamos librerías que vamos a usar

library(tidyverse)
library(broom)
library(ISLR)
library(GGally)
library(modelr)
library(pROC)
library(cowplot)
library(OneR)
library(rlang)
library(caret)
set.seed(1992)

Cargamos los datos de entrenamiento

titanic_complete_train <- read_csv("~/maestria/tps/sabado/titanic_complete_train.csv")

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

glimpse(titanic_complete_train)

Observations: 891
Variables: 12
$ PassengerId <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 2…
$ Survived    <int> 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0…
$ Pclass      <int> 3, 1, 3, 1, 3, 3, 1, 3, 3, 2, 3, 1, 3, 3, 3, 2, 3, 2, 3, 3, 2, 2, 3, 1, 3, 3, 3, 1, 3, 3, 1, 1, 3, 2…
$ Name        <chr> "Braund, Mr. Owen Harris", "Cumings, Mrs. John Bradley (Florence Briggs Thayer)", "Heikkinen, Miss. …
$ Sex         <chr> "male", "female", "female", "female", "male", "male", "male", "male", "female", "female", "female", …
$ Age         <dbl> 22.00000, 38.00000, 26.00000, 35.00000, 35.00000, 26.50759, 54.00000, 2.00000, 27.00000, 14.00000, 4…
$ SibSp       <int> 1, 1, 0, 1, 0, 0, 0, 3, 0, 1, 1, 0, 0, 1, 0, 0, 4, 0, 1, 0, 0, 0, 0, 0, 3, 1, 0, 3, 0, 0, 0, 1, 0, 0…
$ Parch       <int> 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 1, 0, 0, 5, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 5, 0, 2, 0, 0, 0, 0, 0, 0…
$ Ticket      <chr> "A/5 21171", "PC 17599", "STON/O2. 3101282", "113803", "373450", "330877", "17463", "349909", "34774…
$ Fare        <dbl> 7.2500, 71.2833, 7.9250, 53.1000, 8.0500, 8.4583, 51.8625, 21.0750, 11.1333, 30.0708, 16.7000, 26.55…
$ Cabin       <chr> NA, "C85", NA, "C123", NA, NA, "E46", NA, NA, NA, "G6", "C103", NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ Embarked    <chr> "S", "C", "S", "S", "S", "Q", "S", "S", "S", "C", "S", "S", "S", "S", "S", "S", "Q", "S", "S", "C", …

Modificamos el dataset para ajustarlo a los requerimientos

titanic_complete_train <- titanic_complete_train %>% select(PassengerId, Survived, Pclass, Sex, Age, SibSp, Parch, Fare, Embarked)
titanic_complete_train <- titanic_complete_train %>% mutate_at(vars(Survived, Pclass, Embarked), factor)
glimpse(titanic_complete_train)

Observations: 891
Variables: 9
$ PassengerId <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 2…
$ Survived    <fct> 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0…
$ Pclass      <fct> 3, 1, 3, 1, 3, 3, 1, 3, 3, 2, 3, 1, 3, 3, 3, 2, 3, 2, 3, 3, 2, 2, 3, 1, 3, 3, 3, 1, 3, 3, 1, 1, 3, 2…
$ Sex         <chr> "male", "female", "female", "female", "male", "male", "male", "male", "female", "female", "female", …
$ Age         <dbl> 22.00000, 38.00000, 26.00000, 35.00000, 35.00000, 26.50759, 54.00000, 2.00000, 27.00000, 14.00000, 4…
$ SibSp       <int> 1, 1, 0, 1, 0, 0, 0, 3, 0, 1, 1, 0, 0, 1, 0, 0, 4, 0, 1, 0, 0, 0, 0, 0, 3, 1, 0, 3, 0, 0, 0, 1, 0, 0…
$ Parch       <int> 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 1, 0, 0, 5, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 5, 0, 2, 0, 0, 0, 0, 0, 0…
$ Fare        <dbl> 7.2500, 71.2833, 7.9250, 53.1000, 8.0500, 8.4583, 51.8625, 21.0750, 11.1333, 30.0708, 16.7000, 26.55…
$ Embarked    <fct> S, C, S, S, S, Q, S, S, S, C, S, S, S, S, S, S, Q, S, S, C, S, S, Q, S, S, S, C, S, Q, S, C, C, Q, S…

Exploratorias

Realizamos un gráfico exploratorio completo para ver el comportamiento y las relaciones entre las variables. El color rojo designa a quienes no sobreviven y el azul a los que sí.

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

Se observa que las variables Pclass y Sex son buenos discriminantes de Survived, mientras que Age y Fare no.

Analicemos la distribución de la clase de supervivencia

titanic_complete_train %>% group_by(Survived) %>% summarise(numero_casos=n())

El 61% de los pasajeros mueren mientras que el 39% sobreviven. Es un problema con un pequeño desbalance en las clases.

División entre training y validación

Dividimos un 70% del dataset para train y un 30% para test.

train_test  <- titanic_complete_train %>% resample_partition(c(train=0.7, validacion=0.3))
train       <- train_test$train %>% as_tibble()
validacion  <- train_test$validacion %>% as_tibble()

Veamos la distribución de casos en cada uno:

train      %>% group_by(Survived) %>% summarise(numero_casos=n())

validacion %>% group_by(Survived) %>% summarise(numero_casos=n())

La distribución de casos tanto para train como para validacion es muy similar a la original.

Problema

Queremos estimar \(P(Survived=Yes|X)=P(X)\) para cada individuo y partir de ello poder definir un punto de corte para predecir quienes son los que van a sobrevivir.

Regresión logística

Usamos la funcion logistica

Creación de fórmulas.

Generamos varias fórmulas para usar a lo largo del trabajo.

logit_formulas <- formulas(.response = ~Survived, # único lado derecho de las formulas.
                         pclass_sex_age = ~ Pclass + Sex + Age, 
                         pclass_sex = ~Pclass + Sex,
                         pclass = ~Pclass,  
                         todas = ~ Survived + Pclass + Sex + Age + SibSp + Parch + Fare + Embarked
                         )

Creación de modelos

Procedemos a crear los modelos a partir de estas fórmulas

models <- data_frame(logit_formulas) %>% # dataframe a partir del objeto formulas
  mutate(models = names(logit_formulas), # columna con los nombres de las formulas
         expression = paste(logit_formulas), # columna con las expresiones de las formulas
         mod = map(logit_formulas, ~glm(.,family = 'binomial', data = train))) # En mod va a estar el modelo fiteado 

Usamos family = ‘binomial’ porque estamos trabajando con un fenómeno que tiene una distribución binomial. Veamos el primero que es el que necesitamos para este punto.

models %>% 
  filter(models %in% c('pclass_sex_age')) %>%
  mutate(tidy = map(mod,tidy)) %>%  # Qué realizamos en este paso? Que va a tener esta columna?
  unnest(tidy, .drop = TRUE) %>% 
  mutate(estimate=round(estimate,5),
         p.value=round(p.value,4))

Vemos que todos los coeficientes son significativos. Recordando que un coeficiente positivo indica que frente a aumentos de dicha variable la probabilidad aumenta mientras que un coeficiente negativo indica lo contrario:

• El coeficiente del intercept implica que es muy elevada la probabilidad de supervivencia para un miembro de Pclass1 femenino (alrededor de 0.97).
• El coeficiente de Pclass2 negativo indica que disminuye la probabilidad de supervivencia respecto de Pclass1.
• El coeficiente de Pclass3 negativo indica que disminuye su probabilidad de supervivencia respecto de Pclass1 más aún que Pclass2.
• El coeficiente Sexmale negativo indica que ser hombre disminuye la probabilidad de sobrevivir.
• El coeficiente Age es muy cercano a 0, pero indica que a mayor edad, menor la probabilidad de sobrevivir (Mujeres y niños (y ricos!) primero).

Predecimos la probabilidad de sobrevivir de Rose y Jack

rose_jack <- data.frame(Pclass = c("1", "3"), Sex = c("female", "male"), Age = c(17, 20))
prediccion_rose_jack <- augment(x=models$mod$pclass_sex_age, newdata=rose_jack, type.predict='response') 
prediccion_rose_jack

Como era de esperar, la probabilidad de sobrevivir de Rose era mucho mayor que la de Jack (la de Rose cercana a 1 mientras que la de Jack a 0.1).

Veamos qué modelo minimiza la deviance.

models <- models %>% 
  mutate(glance = map(mod,glance))
# Obtener las medidas de evaluacion de interes
models %>% 
  unnest(glance, .drop = TRUE) %>%
  # Calculo de la deviance explicada
  mutate(perc_explained_dev = 1-deviance/null.deviance) %>% 
  select(-c(models, df.null, AIC, BIC)) %>% 
  arrange(deviance)

El modelo que utiliza todas las variables es el que menor deviance obtiene y mayor porcentaje de la deviance explica.

models %>% 
  filter(models %in% c('todas')) %>%
  mutate(tidy = map(mod,tidy)) %>%  # Qué realizamos en este paso? Que va a tener esta columna?
  unnest(tidy, .drop = TRUE) %>% 
  mutate(estimate=round(estimate,5),
         p.value=round(p.value,4))

En este modelo el único coeficiente significativo nuevo que se agrega es EmbarkedS, y nuevamente ser de primera clase y sexo femenino da una altísima probabilidad de sobrevivir, mientras que cualquier otra situación la disminuye.

Predecimos con entrenamiento y curvas ROC

# Añadimos las predicciones
models <- models %>% 
  mutate(pred = map(mod,augment, type.predict = "response"))
#Observaciones con probabilidad más baja
models$pred$todas %>% arrange(.fitted) %>% head(10)

#Observaciones con probabilidades más altas
models$pred$todas %>% arrange(desc(.fitted)) %>% head(10)

Curvas ROC

# Calculamos curvas ROC
roc_todas <- roc(response=models$pred$todas$Survived, predictor=models$pred$todas$.fitted)

Setting levels: control = 0, case = 1
Setting direction: controls < cases

Graficamos

ggroc(list(todas=roc_todas), size=1) + geom_abline(slope = 1, intercept = 1, linetype='dashed') + theme_bw() + labs(title='Curvas ROC', color='Modelo')

print(paste('AUC: Modelo todas', roc_todas$auc))

[1] "AUC: Modelo todas 0.852773131150746"

El AUC es de 0.85, bastante alto. Quizás convendría utilizar precision-recall por ser un dataset con desbalanceo. La sensibilidad aumenta hasta alrededor de 0.7 sin un incremento en los falsos positivos y a partir de ahí encontrar nuevos positivos requiere de muchas falsas detecciones.

Veamos el gráfico de violín:

ggplot(models$pred$todas, aes(x=Survived, y=.fitted, group=Survived,fill=factor(Survived))) + 
  geom_violin() +
  theme_bw() +
  guides(fill=FALSE) +
  labs(title='Violin plot', subtitle='Modelo todas', y='Predicted probability')

Se observa que las clases se separan muy bien, estándo la mayoría de los que no sobreviven por debajo del 0.5 y los que si por arriba. Podemos usar 0.5 como punto de corte.

Elección del punto de corte

Para elegir el punto de corte vamos a usar el dataset de validación. Predecimos la probabilidad de sobrevivir en validación:

prediccion_validacion <- augment(x=models$mod$todas, newdata=validacion, type.predict='response') 
prediccion_validacion

#Observaciones con probabilidades más altas
prediccion_validacion %>% arrange(desc(.fitted)) %>% head(10)

Como sospechábamos, las pasajeras de primera clase más jóvenes tienen las más altas probabilidad de sobrevivir.

#Observaciones con probabilidades más altas
prediccion_validacion %>% arrange(.fitted) %>% head(10)

Mientras que los hombres de tercera más viejos la menor.

Veamos como son las métricas en función del cutoff

prediction_metrics <- function(cutoff, predictions=prediccion_validacion){
    tabla <- predictions %>% 
    mutate(predicted_class=if_else(.fitted>cutoff, 1, 0) %>% as.factor(),
           Survived = factor(Survived))
  
    confusionMatrix(table(tabla$predicted_class, tabla$Survived), positive = "1") %>%
    tidy() %>%
    select(term, estimate) %>%
    filter(term %in% c('accuracy', 'sensitivity', 'specificity', 'precision','recall')) %>%
    mutate(cutoff=cutoff)
  
}
cutoffs = seq(0.05,0.95,0.01)
logit_pred= map_dfr(cutoffs, prediction_metrics) %>% mutate(term=as.factor(term))
ggplot(logit_pred, aes(cutoff,estimate, group=term, color=term)) + geom_line(size=1) +
  theme_bw() +
  labs(title= 'Accuracy, Sensitivity, Specificity, Recall y Precision', subtitle= 'Modelo todas', color="")

Vemos que a medida que aumenta el cutoff aumenta la precisión y especificidad y disminuyen la sensibilidad y la accuracy (el recall es lo mismo que la sensibilidad y por eso no se visualiza). Vemos que entre 0.25 y 0.75 de cutoff la accuracy se plancha. El punto de corte depende de qué estemos buscando en la predicción. Podríamos elegir un cutoff en el cruce de precision - recall, alrededor de 0.45, con una buena capacidad de detección y confianza. Si necesitamos más confianza en la predicción, podríamos usar un cutoff de aproximadamente 0.75 pero perdemos bastante capacidad de recall.

Matriz de confusión

tabla <- prediccion_validacion %>% 
mutate(predicted_class=if_else(.fitted>0.45, 1, 0) %>% as.factor(),
       Survived = factor(Survived))
table(tabla$predicted_class, tabla$Survived)

   
      0   1
  0 139  29
  1  28  72

Vemos que detectamos correctamente 72 de los 101 sobrevivientes (un 72%). Nos perdemos de detectar 29 sobrevivientes y detectamos erroneamente 28 muertos como sobrevivientes.

Veamos como clasificamos los pasajeros de test

titanic_complete_test <- read_csv("~/maestria/tps/sabado/titanic_complete_test.csv")

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

titanic_complete_test <- titanic_complete_test %>% select(PassengerId, Survived, Pclass, Sex, Age, SibSp, Parch, Fare, Embarked)
titanic_complete_test <- titanic_complete_test %>% mutate_at(vars(Survived, Pclass, Embarked), factor)
glimpse(titanic_complete_test)

Observations: 418
Variables: 9
$ PassengerId <int> 892, 893, 894, 895, 896, 897, 898, 899, 900, 901, 902, 903, 904, 905, 906, 907, 908, 909, 910, 911, …
$ Survived    <fct> 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1…
$ Pclass      <fct> 3, 3, 2, 3, 3, 3, 3, 2, 3, 3, 3, 1, 1, 2, 1, 2, 2, 3, 3, 3, 1, 3, 1, 1, 1, 3, 1, 3, 1, 3, 2, 2, 3, 3…
$ Sex         <chr> "male", "female", "male", "male", "female", "male", "female", "male", "female", "male", "male", "mal…
$ Age         <dbl> 34.50000, 47.00000, 62.00000, 27.00000, 22.00000, 14.00000, 30.00000, 26.00000, 18.00000, 21.00000, …
$ SibSp       <int> 0, 1, 0, 0, 1, 0, 0, 1, 0, 2, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 2, 1, 2, 1, 1…
$ Parch       <int> 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 3, 0, 1, 0, 0, 0, 0, 0, 2, 2…
$ Fare        <dbl> 7.8292, 7.0000, 9.6875, 8.6625, 12.2875, 9.2250, 7.6292, 29.0000, 7.2292, 24.1500, 7.8958, 26.0000, …
$ Embarked    <fct> Q, S, Q, S, S, S, Q, S, C, S, S, S, S, S, S, C, Q, C, S, C, C, S, S, C, C, S, C, C, S, C, S, S, S, S…

Predecimos los test con el modelo todas y cutoff 0.45.

prediccion_test <- augment(x=models$mod$todas, newdata=titanic_complete_test, type.predict='response') 
prediccion_test

Matriz de confusión

tabla <- prediccion_test %>% 
mutate(predicted_class=if_else(.fitted>0.45, 1, 0) %>% as.factor(),
       Survived = factor(Survived))
table(tabla$predicted_class, tabla$Survived)

   
      0   1
  0 204  43
  1  57 114

Estamos detectando correctamente un 53%, bastante menos que con el dataset de validación. Esto es razonable porque usamos validación para encontrar el mejor punto de corte, que no tiene por qué ser el de testing. Nos perdemos 43 sobrevivientes y detectamos falsamente 57 muertos como sobrevivientes.