Enfoque Estadístico del Aprendizaje

Trabajo Práctico Nº3

federicomoreno613@gmail.com*
Trabajo sobre los datos de Kaggle de Titanic*

          used  (Mb) gc trigger  (Mb) max used  (Mb)
Ncells 3104495 165.8    5437717 290.5  5437717 290.5
Vcells 6834011  52.2   14786712 112.9 12255457  93.6

1 Preparación de los datos

# Carga de datos
train <- read_csv("C:/Users/fede_/Desktop/Maestria/6- Enfoque Estadistico del aprendizaje/Practico Sabados/titanic_complete_train.csv")

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()
)

a.Leer archivo y mostrar su estructura<br

El archivo se compone de 891 observaciones y 12 variables. La clase está en la columna Survived.

#Muestro estructura del dataset
glimpse(train)

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

b. Seleccion de varaibles.

Se seleccionan las variables PassengerId , Survived , Pclass , Sex , Age , SibSp , Parch, Fare y Embarked

train %<>% # Uso el operador de magrittr que reasigna a la variable tras transformar
  select(PassengerId, Survived, Pclass, Sex,
         Age, SibSp,Parch, Fare, Embarked)

c. Transformar las variables a factor.

train %<>% # Muto las variables elegidas a factor
  mutate_at(vars(Survived, Pclass, Embarked), as_factor)

#muestro nueva estructura
glimpse(train)

Observations: 891
Variables: 9
$ PassengerId <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,...
$ 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...
$ 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...
$ Sex         <chr> "male", "female", "female", "female", "male", "male", "male", "male", "female", "f...
$ Age         <dbl> 22.00000, 38.00000, 26.00000, 35.00000, 35.00000, 26.50759, 54.00000, 2.00000, 27....
$ SibSp       <dbl> 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...
$ Parch       <dbl> 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...
$ Fare        <dbl> 7.2500, 71.2833, 7.9250, 53.1000, 8.0500, 8.4583, 51.8625, 21.0750, 11.1333, 30.07...
$ 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...

# analisis de na,ceros y estructura.
funModeling::status(train)

Al analizar la cantidad de ceros y NA por columna vemos que la única columna con valores faltantes, en el orden del 77% de las filas, es la de Cabina. Las filas de SibSp y Parch presentan una gran cantidad de valores 0.

d. Realizar un gráfico de ggpairs para las variables Survived, Pclass, Sex, Age y Fare e interpretarlo

train %>%
  select(Survived, Pclass, Sex, Age, Fare) %>% # Elijo variables 
  ggpairs(., mapping = aes(color = Survived), legend = c(1,1)) + # Grafico ggpairs abriendo por clase en el color
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  theme_grey()+
  theme(legend.position = "bottom")

Figuara 1: Grafico de GGpairs de variables seleccionadas por clase.

En la figura 1 podemos analizar la distribución de clase según las variables seleccionadas. En rojo están los pasajeros no superiviventes y en azul los que lo hicieron.

Se observa que la clase mayoritaria es la de no sobreviviente.

Pclass determina una diferencia importante entre los pasajeros de la 3er clase y los demás. Por un lado son mayoriatria y por otro lado también son los que menos sobrevivieron. Ser mujer en primera o segunda clase casi aseguró la supervivencia, y si era pasajera de tercera la tasa estaba aproximadamente en el 50%. Los caballeros de primera y segunda clase tuvieron también una tasa de supervivencia alta

Sex la proporción de varones es mayor en el caso de los sobrevivientes. Sin embargo pueden estar sobrerpresentados entre los que no lograron salvar sus vidas.

Age y Price Analizando el coeficiente de correlacion de Pearson de 0.118 no es signficativo, aunque si se observa en la tercera clase la distribución por edad representa a personas más jovenes

e. Mostrar la distribución de clase (Sobrevivientes vs No Sobrevivientes)
Analicemos la distribución de la clase. La proporción de clase es de aprox 62% para los que no sobrevivieron y 38% de casos de que si sobrevivieron

# agrupo por variable target
train %>% group_by(Survived) %>% summarise(numero_casos=n(), perc = round(n() / nrow(train) *100,2))

NA

f. Datasets de entrenamiento y validación

# Separo training y validation con modelr
train_test <- train %>% resample_partition(c(train = 0.7, test = 0.3))

# Convierto en tibbles de tidyverse para construir el dataset
training <- train_test$train %>% as_tibble()
validation <- train_test$test %>% as_tibble()

# Muto a factores los campos correspondientes.
training %<>%
  mutate_at(vars(Survived, Pclass, Embarked), as_factor)

validation %<>%
  mutate_at(vars(Survived, Pclass, Embarked), as_factor)

# Uno las cuentas y proporciones por clase del dataset original con las particiones hechas--
rbind(
  train %>%
    group_by(Survived) %>%
    summarise(
      dataset = "dataset_completo",
      count = n(),
      prop = n() / nrow(train)
    ) %>%
    select(dataset, everything())
  ,
  training %>%
    group_by(Survived) %>%
    summarise(
      dataset = "training",
      count = n(),
      prop = n() / nrow(training)
    ) %>%
    select(dataset, everything())
  ,
  validation %>%
    group_by(Survived) %>%
    summarise(
      dataset = "validation",
      count = n(),
      prop = n() / nrow(validation)
    ) %>%
    select(dataset, everything())
)

NA

2 Predicciones

a. Generación del primer modelo

Genero el listado de fórmulas con los modelos que voy a utilizar cómo optativos

logit_formula <- 
  # creamos un objeto que contiene todas las fórmulas que vamos a utilizar
  formulas(.response = ~Survived,
                         modelobase = ~Pclass+Sex+Age,
                         modelo1= ~Embarked+Fare+Pclass,
                         modelo2 = ~Sex+SibSp,
                         modelo3 = ~Pclass
                          )

Creación del modelo

# Genero los modelos de regresión logística mapeando a glm las distintas fórmulas recién creadas
models <- data_frame(logit_formula) %>%
  mutate(model = names(logit_formula),
         expression = paste(logit_formula),
         mod = map(logit_formula, ~glm(.,family = 'binomial', data = training))) 

b. Interpretación de coeficientes

models %>%
  # Filtro solo el modelo base
  filter(model == "modelobase") %>%
  # Mapeo tidy de la librería broom para obtener coeficientes
  mutate(tidy = map(mod,tidy)) %>% 
  # Hago el unnest de la información obtenida
  unnest(tidy, .drop = TRUE) %>% 
  # Elijo columnas de interés
  select(term, estimate, std.error, statistic, p.value)

Todos los coeficientes resultan significativos respecto a la predicción “survived”

Los coeficientes asociados a las variables tienen signo negativo, a lo cual resulta que un incremento, en dichas variables, disminuye la probabilidad de sobrevivir

PClass Es mayor la probabilidad de no sobrevivir siendo de tercera clase que de segunda clase.
Sex La categoría basal respecto al sexo es mujer, lo cual significa que el hecho de ser hombre reduce la probabilidad de sobrevivir y si bien el coeficiente de Age es pequeño es significativo y también es negativo, por lo tanto a mayor edad disminuye la probabilidad de sobrevivir. Esto daría verosimil a la frase “Women and children first”

c. Probabilidad de supervivencia de Rose y Jack

A priori observando los cieficientes, dado que Jack es de tercera clase y hombre debería ser más probable que no sobreviva que rose. Hago el testeo:

# Construyo nuevo dataset con tribble
new_data <- tribble(~Pclass, ~Sex, ~Age, ~Name,
                   "1", "female", 17, "Rose",
                   "3", "male", 20, "Jack")


models %>%
  filter(model == "modelobase") %>%
  # Mapeo augment al nuevo modelo, usando la nueva data para obtener predicciones
  mutate(pred = map(mod, augment, newdata = new_data, type.predict = "response")) %>%
  # Hago el unnest de los resultados obtenidos del augment
  unnest(pred) %>%
  # Elijo variables
  select(Name, Sex, Pclass, Age, .fitted)

El resultado avala nuestra interpretación

3 Generación de modelos

a. Generación de 3 modelos adicionales

Se trabaja con los 3 modelos adicionales generados antes en la variable ‘logit_formula’:

Me interesa conocer todos los atributos económicos de la persona por eso genero el modelo1

modelo1= ~Embarked+Fare+Pclass

En el modelo 2 me interesa saber la relación entre genero y relaciones familiares SibSp

modelo2 = ~Sex+SibSp

En el tercer modelo me interesa sólo saber la relación con la clase económica o de pasaje y su probabilidad de supervivencia<br

modelo3 = ~Pclass

models %>%
  # Filtro modelo de interés
  filter(model != "modelo1") %>%
  # Obtengo coeficientes
  mutate(tidy = map(mod, tidy)) %>%
  # Unnest de la info de tidy
  unnest(tidy, .drop = TRUE) %>%
  # Elijo variables
  select(model, term, estimate, std.error, statistic, p.value) %>%
  # Chequeo la significatividad de los coeficientes
  mutate(is_sign = p.value <= 0.05) %>%
  filter(term != "(Intercept)")

Todos los coeficientes resultan significativos excepto Pclass2 para el modelo3 que sólo tenia la clase a la que pertenecia<br

b. Selección del mejor modelo<br

Hago elección de modelo considerando cómo mejor aquel que minimiza su valor o que es lo mismo, maximiza la deviance explicada. En este caso sería el modelo base<br

models %<>%
  # Mapeo glance para obtener medidas de performance de los modelos
  mutate(glance = map(mod, broom::glance))

models %>%
  # Unnest de estas medidas
  unnest(glance, .drop = TRUE) %>%
  # Calculo el porcentaje de deviance explicada por cada modelo
  mutate(perc_explained_dev = 1 - deviance / null.deviance) %>%
  # Elijo variables
  select(expression, null.deviance, deviance, perc_explained_dev) %>%
  # Ordeno por deviance
  arrange(deviance)

4 Evaluación del modelo

a.Curva ROC y AUC<br

Calculo preddiciones<br

# Añadimos predicciones a los datos de training
models <- models %>%
  mutate(pred = map(mod, augment, type.predict = "response"))

# Unnest de las de nuestro modelo de interés
prediction_mbase <-
  models %>%
  filter(model == "modelobase") %>%
  unnest(pred, .drop = TRUE)

Se dibuja la curva roc del mejor modelo<br

# Construyo la curva ROC
roc_mbase <-
  roc(response = prediction_mbase$Survived, predictor = prediction_mbase$.fitted)

# La grafico
ggroc(roc_mbase, size = 1.2) +
  geom_abline(slope = 1,
              intercept = 1,
              linetype = 'dashed') +
  theme_grey() +
  labs(
    title = 'Curva ROC',
    caption = paste('AUC:', round(roc_pc_sex_age$auc, 3)),
    subtitle = "Modelo usando Pclass, Sex y Age"
  )

Error in paste("AUC:", round(roc_pc_sex_age$auc, 3)) : 
  object 'roc_pc_sex_age' not found

La curva Roc esta dibujada sobre los ejes que representan la especificidad (True negative rate) y Sensitivity (True positive rates), es decir, existe un trade off entre la proporción de positivos (sobrevivientess) clasificados correctamente versus la proporcion de negativos (no sobrevivientes)

Se puede observar que nuestro modelo está por encima de la linea punteada que representa la probabilidad del azar.

b. Violin Plot

# Grafico predicciones
ggplot(prediction_pc_sex_age,
       aes(
         x = Survived,
         y = .fitted,
         group = Survived,
         fill = factor(Survived)
       )) +
  geom_violin() +
  theme_grey() +
  guides(fill = FALSE) +
  labs(title = 'Violin plot',
       subtitle = 'Modelo usando Pclass, Sex y Age',
       y = 'Predicted probability') +
  # Agrego puntos con jittering
  geom_jitter(
    shape = 16,
    size = 1,
    position = position_jitter(0.2),
    color = "black",
    alpha = .7
  )

NA

Se observa que el modelo es bastante bueno para realizar predicciones dado que para cada clase los “violines” se encuentran invertidos. A priori pareciera que 0.5 es un buen punto de corte para efectuar predicciones.

5 Elección del punto corte

a.Gráficos de métricas de performance<br

Analizamos distintas métricas a continuación, para analizar el punto de corte.

Recall o Sensitivity (Proporción de verdaderos positivosclasificados cómo tal) Accuracy (proporción de casos clasificados correctamente sobre el total de casos) Precision (proporción del total de casos clasificados como positivos clasificados correctamente).

Primero se generan las predicciones para el dataset de validación

# Obtenemos las correspondientes a nuestro mejor modelo
prediction_validation_modelobase <-
  models_val %>%
  filter(model == "modelobase") %>%
  unnest(pred, .drop = TRUE)

Error in eval(lhs, parent, parent) : object 'models_val' not found

Luego se grafica:

b. Elección del punto de corte óptimo<br

El punto de corte a elegir es en donde se cruzan precision y recall (sensitivity). La idea es predecir la mayor cantidad de supervivientes correctamente. Apróximadamente el punto de corte óptimo es de 0.48

# Me quedo con registros cuya diferencia entre sensitividad y especificidad es mínima, dado que ninguno es estrictamente 0
logit_pred %>% 
  # Función nueva de tidyr que reemplaza a spread()
  pivot_wider(names_from = term, values_from = estimate) %>% 
  arrange(-accuracy) %>% 
  filter(abs(sensitivity - specificity)==min(abs(sensitivity - specificity)) ) %>% 
  select(-c("precision","f1"))

# Idem con registros en que la accuracy es máxima
logit_pred %>% 
  pivot_wider(names_from = term, values_from = estimate) %>% 
  filter(accuracy == max(accuracy))%>% 
  select(-c("precision","f1"))

# Comparando ambos encuentro aprox. el punto de corte óptimo

Elegir el punto de corte depende de cual es el objetivo de nuestra investigación. Si nuestro foco está en los positivos queremos aumentar la sensitividad y la precisión deberia elegir el punto de corte de 0.48. En cambio si está en los negativos queremos maximizar el recall y deberia elegir el cutoff de 0.35.

Una manera de establecer el equilibrio es elegir a partir de la métrica de F Score y por lo tanto me quedaría con el punto de corte en 0.48.

# F1 Score
logit_pred %>%
  pivot_wider(names_from = term, values_from = estimate) %>%
  arrange(-f1) %>%
  head(1) %>%
  select(cutoff, f1, sensitivity, precision)

c. Matriz de confusión de validación<br

La matriz de confusión permite comparar la clase predicha contra la clase que realmente era. Predije correctamente 147 no sobrevivientes de 170 (82%) y predije correctamente 75 sobrevivientes de 98 (76%) PRedigo mejor a los no sobrevivientes que eran la clase mayoritaria

sel_cutoff = 0.480963   

# Creo el modelo nuevamente, para simplificar
modelo <-
  glm(logit_formula$modelobase, family = 'binomial', data = training)
# Agrego las predicciones al dataset de validación
table = augment(x = modelo,
                newdata = validation,
                type.predict = 'response')
# Clasifico utilizando el punto de corte
table = table %>% mutate(
  predicted_class = if_else(.fitted >= sel_cutoff, 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 147  23
         1  23  75
                                          
               Accuracy : 0.8284          
                 95% CI : (0.7778, 0.8715)
    No Information Rate : 0.6343          
    P-Value [Acc > NIR] : 2.385e-12       
                                          
                  Kappa : 0.63            
                                          
 Mcnemar's Test P-Value : 1               
                                          
            Sensitivity : 0.7653          
            Specificity : 0.8647          
         Pos Pred Value : 0.7653          
         Neg Pred Value : 0.8647          
             Prevalence : 0.3657          
         Detection Rate : 0.2799          
   Detection Prevalence : 0.3657          
      Balanced Accuracy : 0.8150          
                                          
       'Positive' Class : 1               
                                          

6 Dataset de testeo

a. Lectura y transformación del dataset de testing<br

# Leo el dataset de test, transformo a factores y me quedo con las variables relevantes

test <- read_csv("titanic_complete_test.csv")

test %<>%
  select(PassengerId, Survived, Pclass, Sex,
         Age, SibSp, Parch, Fare, Embarked) %>%
  mutate_at(vars(Survived, Pclass, Embarked), as_factor)

head(test)

b. Clasificación de individuos<br

Clasfico a los pasajeros del dataset de test utilizando el modelo base con el punto de corte elegido

# Agrego la predicciones al dataset de testeo
table_test = augment(x = modelo,
                     newdata = test,
                     type.predict = 'response')
# Clasifico utilizando el punto de corte
table_test = table_test %>% mutate(
  predicted_class = if_else(.fitted >= sel_cutoff, 1, 0) %>% as.factor(),
  Survived = factor(Survived)
)

# Muestro la clasificación de cada persona

table_test %>%
  select(PassengerId, Pclass, Sex, Age, prob = .fitted, predicted_class) %>%
  arrange(-prob)

NA
NA

c. Matriz de confusión de test<br

Generamos la matriz de confusión del test

# Creo la matriz de confusión
confusionMatrix(table_test$predicted_class, table_test$Survived, positive = "1")$table

          Reference
Prediction   0   1
         0 212  50
         1  49 107

# Genero métricas para validación nuevamente
logit_val = map_dfr(sel_cutoff, prediction_metrics) %>%
  mutate(
    estimate = round(estimate, 3),
    dataset = "validation",
    term = as.factor(term)
  ) %>%
  pivot_wider(names_from = term, values_from = estimate) %>%
  select(-cutoff)


# Genero métricas para test
logit_test = map_dfr(sel_cutoff, prediction_metrics, table_test) %>% mutate(
  estimate = round(estimate, 3),
  dataset = "test",
  term = as.factor(term)
) %>%
  pivot_wider(names_from = term, values_from = estimate) %>%
  select(-cutoff)

# Las coloco una abajo de la otra
rbind(logit_val, logit_test)

NA

El modelo obtenidos con el dataset de validación tiene mejor performace que el que solo utiliza el dataset de training. Esto es porque optimizamos el punto de corte con el de validación. El acuracy es 6% mayor y además todas las metricas suben