TP03 - Estefania De Marzio

library(dplyr)
library(tidymodels)
library(GGally)
library(purrr)
library(kableExtra)
library(modelr)
library(pROC)
library(caret)

Preparación de datos

Se debe armar un modelo de regresión logística para clasificar si una persona que viajaba en el Titanic sobrevivió o no.

Leer archivo

Se leen los datos y se muestra su estructura.

#leer archivo y mostrar estructura
titanic_train <- read.csv(file = "titanic_complete_train.csv")
glimpse(titanic_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...
$ 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,...
$ 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,...
$ Name        <fct> "Braund, Mr. Owen Harris", "Cumings, Mrs. John Bradley (Florence Briggs Thayer)", "Heikki...
$ Sex         <fct> male, female, female, female, male, male, male, male, female, female, female, female, mal...
$ Age         <dbl> 22.00000, 38.00000, 26.00000, 35.00000, 35.00000, 26.50759, 54.00000, 2.00000, 27.00000, ...
$ 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,...
$ 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,...
$ Ticket      <fct> A/5 21171, PC 17599, STON/O2. 3101282, 113803, 373450, 330877, 17463, 349909, 347742, 237...
$ Fare        <dbl> 7.2500, 71.2833, 7.9250, 53.1000, 8.0500, 8.4583, 51.8625, 21.0750, 11.1333, 30.0708, 16....
$ Cabin       <fct> NA, C85, NA, C123, NA, NA, E46, NA, NA, NA, G6, C103, NA, NA, NA, NA, NA, NA, NA, NA, NA,...
$ 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,...

Seleccionar Variables

Se seleccionan las variables solicitadas y se transforman a factor las features de Survived, Pclass y Embarked.

titanic_train <-
  titanic_train %>% select(
    'PassengerId',
    'Survived',
    'Pclass',
    'Sex',
    'Age',
    'SibSp',
    'Parch',
    'Fare' ,
    'Embarked'
  )
#transformar las variables a factor
titanic_train[, c('Survived', 'Pclass', 'Embarked')] <-
  map(titanic_train[, c('Survived', 'Pclass', 'Embarked')], as.factor)

Gráfico de GGPairs

Se realiza gráfico de ggpairs y se observa que los que viajaban en 3rd clase sobrevivieron en menor proporción que los de 2da y 1era clase. Sobrevivieron más las mujeres que los hombres. No se observan diferencias significativas por Edad y precio del boleto.

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

Distribución de la clase

Se puede ver que se tiene una distribucion desbalancea con 62% de negativos y 38% de positvos.

titanic_train %>% group_by(Survived) %>% summarise(numero_casos=n()/nrow(titanic_train))

ggplot(data=titanic_train, aes(titanic_train$Survived)) + geom_bar() + theme_bw() +
  labs(title="Titanic Train", x="Survived")

#Distribucion de la clase de titanic train.
titanic_train %>% group_by(Survived) %>% summarise(numero_casos=n()/nrow(titanic_train)*100)

Dividir el conjunto de Train

Se observa que luego deñ resample_partition, los datos tienen la misma distribución de clase que la orginal.

#Separo en train y validacion
train_test <- titanic_train %>% resample_partition(c(train=0.7,test=0.3))

titanic_train2 <- train_test$train %>% as_tibble()
titanic_validation <- train_test$test %>% as_tibble()
ggplot(data=titanic_train2, aes(titanic_train2$Survived)) + geom_bar() + theme_bw() +
  labs(title="Titanic Train 70%", x="Survived")

ggplot(data=titanic_validation, aes(titanic_validation$Survived)) + geom_bar() + theme_bw() +
  labs(title="Titanic Validation 30%", x="Survived")

#Distribucion de la clase de titanic train2.
titanic_train2 %>% group_by(Survived) %>% summarise(numero_casos=n()/nrow(titanic_train2)*100)

#Distribucion de la clase de titanic validation
titanic_validation %>% group_by(Survived) %>% summarise(numero_casos=n()/nrow(titanic_validation)*100)

NA

Predicciones

Modelo Pclass+Sex+Age

Se generan varios modelos incluyendo el que se solicita en este punto para mayor practicidad.

logit_formulas <- formulas(.response = ~Survived, 
                         pc_sex_age=~Pclass+Sex+Age, 
                         pc_fare_emb=~Pclass+Fare+ Embarked , 
                         sex_fare= ~ Sex + Fare,
                         pclass= ~Pclass
                         )
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 = titanic_train2)))

`data_frame()` is deprecated, use `tibble()`.
This warning is displayed once per session.

Interpretar coeficientes

Todos los coeficientes resultaron significativos en este modelo. Las variables tienen coeficientes negativos lo cual indica que disminuye la probabilidad esperada de sobrevivir a medida que éstas aumentan.

models %>% 
  filter(models == 'pc_sex_age') %>%
  mutate(tidy = map(mod,tidy)) %>%  
  unnest(tidy, .drop = TRUE) %>% 
  select(term, estimate, std.error, statistic, p.value)

The `.drop` argument of `unnest()` is deprecated as of tidyr 1.0.0.
All list-columns are now preserved.
This warning is displayed once per session.
Call `lifecycle::last_warnings()` to see where this warning was generated.

Quién tiene mayor probabilidad de Supervivencia

Se observa que Jack tenia muy bajas chances de sobrevivir aunque el trozo de madera que utlizo Rose podria haber salvado a ambos.

newdata <- tibble(Age = c(17,20), Pclass = c('1','3'), Sex = c('female','male'))

models %>% 
  filter(models == 'pc_sex_age') %>%
  mutate(pred= map(mod,augment, newdata= newdata,type.predict = "response")) %>%  
  unnest(pred, .drop = TRUE) %>%
  select(Age, Pclass, Sex, .fitted)

NA

Generación de Modelos

Generar 3 modelos

Se generaron tres modelos:

• Pclass+Fare+ Embarked
• Sex + Fare
• Pclass

Se observa que el modelo ~Pclass+Fare+ Embarked, EmbarkedQ y Pclass2 no son significativos.

models %>% 
  filter(models != 'pc_sex_age') %>%
  mutate(tidy = map(mod,tidy)) %>%  
  unnest(tidy, .drop = TRUE) %>% 
  select(term, estimate, std.error, statistic, p.value)

Ordernar por deviance

Elijo el modelo que maximice el porcentaje de la desviance explicada. En este caso, el mejor modelo es ~ Pclass + Sex + Age.

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

models %>% 
  unnest(glance, .drop = TRUE) %>%
  mutate(perc_explained_dev = 1-deviance/null.deviance) %>% 
  select(expression,null.deviance, deviance, perc_explained_dev) %>% 
  arrange(deviance)

Evaluacion del Modelo

Curva ROC y AUC

Se obtuvo un auc de 0.8488 que es mayor al azar 0.5; aunque este valor se obtiene sobre el conjunto de entrenamiento. La curva ROC visualiza la sensibilidad de la predicción versus la especificidad de la misma. La sensiblidad son las predicciones positivas sobre la cantidad total de positivos y la especificidad son las predicciones negativas sobre la cantidad total de negativos.

models <- models %>% 
  mutate(pred= map(mod,augment, type.predict = "response"))


# Calculamos curvas ROC
prediction_pc_sex_age <- models %>% 
  filter(models=="pc_sex_age") %>% 
  unnest(pred, .drop=TRUE)

roc_pc_sex_age <- roc(response=prediction_pc_sex_age$Survived, predictor=prediction_pc_sex_age$.fitted)

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

roc_pc_sex_age 

Call:
roc.default(response = prediction_pc_sex_age$Survived, predictor = prediction_pc_sex_age$.fitted)

Data: prediction_pc_sex_age$.fitted in 389 controls (prediction_pc_sex_age$Survived 0) < 234 cases (prediction_pc_sex_age$Survived 1).
Area under the curve: 0.8485

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

Violin Plot

Se puede observar a través del gráfico que las probabilidades de las observaciones de la clase no sobreviviente tiene un concentracion entre 0.25 y 0. Mientras que las probabilidades de las observaciones de los sobrevivientes tienen mayor concentracion en 1. Definir un punto de corte que prediga bien ambas clase de forma tal de cometer el mejor error posible, pareciera no dislumbrase clarametne en el grafico.

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

Eleccion Punto de Corte

Gráfico en funcion del punto de corte

pred_val<-models %>% 
  filter(models == 'pc_sex_age') %>%
  mutate(pred= map(mod,augment, newdata= titanic_validation,type.predict = "response")) %>% 
  unnest(pred, .drop=TRUE)

prediction_metrics <- function(cutoff=0.95, predictions=pred_val){
  table <- predictions %>% 
    mutate(predicted_class=if_else(.fitted>cutoff, 1, 0) %>% as.factor(),
           Survived= factor(Survived))
  
  confusionMatrix(table$predicted_class, table$Survived, positive = "1") %>%
    tidy() %>%
    select(term, estimate) %>%
    filter(term %in% c('accuracy', 'sensitivity', 'specificity', 'precision','recall')) %>%
    mutate(cutoff=cutoff)
  
}
cutoffs = seq(0.01,0.95,0.01)

logit_pred= map_dfr(cutoffs, prediction_metrics)%>% mutate(term=as.factor(term))

Levels are not in the same order for reference and data. Refactoring data to match.Levels are not in the same order for reference and data. Refactoring data to match.Levels are not in the same order for reference and data. Refactoring data to match.

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 Pclass+Sex+Age', color="") +
  scale_x_continuous(breaks =seq(0.05,0.95,0.1))

Elegir el Punto de Corte

La elección se basa en la métrica que queremos maximizar. En este caso buscamos la intersección de sensivivity y specifity de forma tal de obtener un balance entre los positivos y negativos. Elegimos como punto de corte 0.4.

Matriz de Confusión

La tabla nos muestra que se clasifico correctamente 117 personas como no sobrevivientes y 84 personas como sobrevivientes. 24 personas las clasifico como no sobrevientes cuando enrealidad si sobrevivieron y 43 personas como sobrevivientes cuando no sobrevivieron.

table <- pred_val %>%
  mutate(
    predicted_class = if_else(.fitted > 0.4, 1, 0) %>% as.factor(),
    Survived = factor(Survived)
  )

confusionMatrix(table$predicted_class, table$Survived, positive = "1")

Confusion Matrix and Statistics

          Reference
Prediction   0   1
         0 117  24
         1  43  84
                                          
               Accuracy : 0.75            
                 95% CI : (0.6937, 0.8007)
    No Information Rate : 0.597           
    P-Value [Acc > NIR] : 1.033e-07       
                                          
                  Kappa : 0.4949          
                                          
 Mcnemar's Test P-Value : 0.02787         
                                          
            Sensitivity : 0.7778          
            Specificity : 0.7312          
         Pos Pred Value : 0.6614          
         Neg Pred Value : 0.8298          
             Prevalence : 0.4030          
         Detection Rate : 0.3134          
   Detection Prevalence : 0.4739          
      Balanced Accuracy : 0.7545          
                                          
       'Positive' Class : 1               
                                          

Testeo

Leeer archivo de Testeo

titanic_test <- read.csv(file = "titanic_complete_test.csv")
titanic_test[, c('Survived', 'Pclass', 'Embarked')] <-
  map(titanic_test[, c('Survived', 'Pclass', 'Embarked')], as.factor)

Clasificar a las personas de Testing

pred_test<-models %>% 
  filter(models == 'pc_sex_age') %>%
  mutate(pred= map(mod,augment, newdata= titanic_test,type.predict = "response")) %>% 
  unnest(pred, .drop=TRUE)

table <- pred_test %>%
  mutate(
    predicted_class = if_else(.fitted > 0.4, 1, 0) %>% as.factor(),
    Survived = factor(Survived)
  )

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

Matriz de Confusión

En el conjunto de test el modelo clasificó correctamente 194 personas como no sobrevivientes y 115 personas como sobrevivientes. 42 personas las clasifico como no sobrevientes cuando enrealidad si sobrevivieron y 67 personas como sobrevivientes cuando no sobrevivieron. En comparacion con las métricas obtenidas sobre el conjunto de validacion, la sensibilidad, especificidad y el acurrancy son similares.

• Test: Acurrancy: 0.7392 , Sensitivity: 0.7325 y Specificity: 0.7433
• Validacion: Acurrancy: 0.75, Sensitivity: 0.7778 y Specificity: 0.7312

confusionMatrix(table$predicted_class, table$Survived, positive = "1")

Confusion Matrix and Statistics

          Reference
Prediction   0   1
         0 194  42
         1  67 115
                                          
               Accuracy : 0.7392          
                 95% CI : (0.6943, 0.7807)
    No Information Rate : 0.6244          
    P-Value [Acc > NIR] : 4.341e-07       
                                          
                  Kappa : 0.4611          
                                          
 Mcnemar's Test P-Value : 0.02152         
                                          
            Sensitivity : 0.7325          
            Specificity : 0.7433          
         Pos Pred Value : 0.6319          
         Neg Pred Value : 0.8220          
             Prevalence : 0.3756          
         Detection Rate : 0.2751          
   Detection Prevalence : 0.4354          
      Balanced Accuracy : 0.7379          
                                          
       'Positive' Class : 1