1. Preparación de los datos

Se cargan los paquetes

library(pROC)
library(dplyr)
library(tidyr)
library(readr)
library(GGally)
library(ggthemes)
library(ggplot2)
library(broom)
library(purrr)
library(modelr)
library(caret)
theme_set(theme_bw()) #carga un theme para todos los gráficos
  1. Leer el archivo y mostrar su estructura
df <- read_csv("D:/OneDrive/Ciencia de Datos/06 Enfoque estadistico/TP01/TP03/titanic_complete_train.csv")
glimpse(df)
## Observations: 891
## Variables: 12
## $ PassengerId <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, ...
## $ Survived    <dbl> 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0...
## $ Pclass      <dbl> 3, 1, 3, 1, 3, 3, 1, 3, 3, 2, 3, 1, 3, 3, 3, 2, 3, 2, 3...
## $ Name        <chr> "Braund, Mr. Owen Harris", "Cumings, Mrs. John Bradley ...
## $ Sex         <chr> "male", "female", "female", "female", "male", "male", "...
## $ Age         <dbl> 22.00000, 38.00000, 26.00000, 35.00000, 35.00000, 26.50...
## $ SibSp       <dbl> 1, 1, 0, 1, 0, 0, 0, 3, 0, 1, 1, 0, 0, 1, 0, 0, 4, 0, 1...
## $ Parch       <dbl> 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 1, 0, 0, 5, 0, 0, 1, 0, 0...
## $ Ticket      <chr> "A/5 21171", "PC 17599", "STON/O2. 3101282", "113803", ...
## $ Fare        <dbl> 7.2500, 71.2833, 7.9250, 53.1000, 8.0500, 8.4583, 51.86...
## $ Cabin       <chr> NA, "C85", NA, "C123", NA, NA, "E46", NA, NA, NA, "G6",...
## $ Embarked    <chr> "S", "C", "S", "S", "S", "Q", "S", "S", "S", "C", "S", ...
  1. Seleccionar las variables PassengerId, Survived, Pclass, Sex, Age, SibSp,Parch, Fare y Embarked
df = select(df, PassengerId, Survived, Pclass, Sex, Age, SibSp,Parch, Fare, Embarked)
  1. Transformar las variables Survived, Pclass y Embarked a factor. (Se añade la variable Sex)
charcol = c('Survived', 'Pclass', "Embarked", 'Sex')
df[charcol] <- lapply(df[charcol], factor)
glimpse(df)
## Observations: 891
## Variables: 9
## $ PassengerId <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, ...
## $ Survived    <fct> 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0...
## $ Pclass      <fct> 3, 1, 3, 1, 3, 3, 1, 3, 3, 2, 3, 1, 3, 3, 3, 2, 3, 2, 3...
## $ Sex         <fct> male, female, female, female, male, male, male, male, f...
## $ Age         <dbl> 22.00000, 38.00000, 26.00000, 35.00000, 35.00000, 26.50...
## $ SibSp       <dbl> 1, 1, 0, 1, 0, 0, 0, 3, 0, 1, 1, 0, 0, 1, 0, 0, 4, 0, 1...
## $ Parch       <dbl> 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 1, 0, 0, 5, 0, 0, 1, 0, 0...
## $ Fare        <dbl> 7.2500, 71.2833, 7.9250, 53.1000, 8.0500, 8.4583, 51.86...
## $ Embarked    <fct> S, C, S, S, S, Q, S, S, S, C, S, S, S, S, S, S, Q, S, S...
  1. Realizar un gráfico de ggpairs para las variables Survived, Pclass, Sex, Age y Fare e interpretarlo
pairs = df %>% select(Survived, Pclass, Sex, Age, Fare)

#función para facilitar la visualización de la distribución de variables continuas
my_dens <- function(data, mapping, ..., low = "#132B43", high = "#56B1F7") {
  ggplot(data = data, mapping=mapping) +
    geom_density(..., alpha=0.7) }
  

ggpairs(pairs,mapping = aes(colour= Survived),diag=list(continuous=my_dens)) + 
    theme(axis.text.x = element_text(angle = 90, hjust = 1))

De la imagen anterior se puede interpretar lo siguiente:

-Existe una proporción mayor de personas que no sobrevivieron que las que sobrevivieron.

-La clase mayoritaria es la #3, la proporción de no sobrevivientes aumenta conforme baja la clase. La tercera clase es la que tiene proporcionalmente menos sobrevivientes.

-Una gran proporción de las mujeres de la primera y segunda clase sobrevivieron. Los hombres de la primera clase sobrevivieron proporcionalmente más que los de segunda y tercera.

-La mayoría de las mujeres sobrevivieron y la mayoría de los hombres perecieron. En el viaje habían más viajeros hombres que mujeres.

-La variable edad agregada no parece ser un buen indicador de la sobrevivencia.

-Si la edad se divide según la clase, se puede observar que las personas más jóvenes en cada clase sobrevivieron más que las más viejas.

-En la primera y segunda clase existe una mayor dispersión de las edades.

-Los hombres sobrevivientes tienen una mayor edad que los no sobrevivientes. No parece haber una diferencia marcada en las edades de las mujeres sobrevivientes y no sobrevivientes.

-La edad disminuye de la primera a la tercera clase. Siendo los más jóvenes los de tercera clase.

-Al comparar la distribución de la edad de los viajeros y su condición de sobrevivencia se puede ver que los niños sobrevivieron en mayor proporción y las personas cerca de los 20 años murieron en mayor proporción.

-La variable tarifa permite hacer un análisis similar al de las clases en donde se observa que los que pagaron una baja tarifa murieron en mayor proporción.

  1. Mostrar la distribución de clase (Sobrevivientes vs No Sobrevivientes)
prop.total = prop.table(table(df$Survived))*100
print(prop.total)
## 
##        0        1 
## 61.61616 38.38384

Se observa que el 62% de los tripulantes murieron en el evento.

  1. 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 este balanceado
set.seed(0)
train_val <- df %>% resample_partition(c(train=0.7,val=0.3))

df_train <- train_val$train %>% as_tibble()
df_val <- train_val$val %>% as_tibble()

prop.train = prop.table(table(df_train$Survived))*100
prop.val = prop.table(table(df_val$Survived))*100

proporciones = rbind(prop.total, prop.train, prop.val)
colnames(proporciones) = c('No sobreviviente', 'Sobreviviente')
proporciones
##            No sobreviviente Sobreviviente
## prop.total         61.61616      38.38384
## prop.train         61.95827      38.04173
## prop.val           60.82090      39.17910

Se observa que la distribución de las clases es muy similar entre los grupos.

2. Predicciones

  1. 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.
mod01 <- glm(Survived ~ Pclass + Sex + Age, data = df_train, family = "binomial")
  1. Dar una breve interpretación de los coeficientes y su significatividad.
tidy(mod01)
## # A tibble: 5 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   3.74     0.467        8.01 1.17e-15
## 2 Pclass2      -1.16     0.317       -3.65 2.58e- 4
## 3 Pclass3      -2.67     0.316       -8.44 3.08e-17
## 4 Sexmale      -2.55     0.225      -11.3  1.39e-29
## 5 Age          -0.0376   0.00929     -4.05 5.14e- 5

En primer lugar se observa que todos los coeficientes estimados son menores a 0,05 y por ende se puede concluir que son significativamente diferentes de 0.

Todos los coeficientes de las variables elgidas son negativas. Esto implica cuando la variable en cuestión aumenta, la probabilidad estimada de sobrevir disminuye. En estos casos nótese que las clases basales son Pclase1 y Sexo Femenino.

Si se miran las clases se puede establecer esta relación en cuanto a la probabilidad estimada de sobrevivir Clase 1 > Clase 2 > Clase 3 (manteniendo todo constante).

Ser mujer aumenta las probablidades estimadas de sobrevir respecto a los hombres manteniendo todo constante. A sy vez, aumento en la edad es perjudicial para las probabilidades de sobrevivencia manteniendo las demás variables constantes.

  1. ¿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.
#se genera una matriz con las propiedades de los individuos
rose_vs_jack = data.frame(id = c(1,2),
                        Pclass = as.factor(c(1,3)),
                        Sex = as.factor(c('female', 'male')),
                        Age = c(17, 20))

#se predice probabilidad de sobrevivencia
predict(mod01, rose_vs_jack, type = 'response')
##          1          2 
## 0.95674501 0.09698184

Las predicciones indican que Rose tiene 95.7% probabilidades de sobrevivir y Jack solo el 9.7%.

3. Generación de modelos

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. 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

#indico las combinaciones de los modelos. el mod 01 corresponde al modelo indicado en el TP.
logit_formulas <- formulas(.response = ~Survived, 
                           mod01= ~Pclass + Sex + Age, 
                           mod02= ~Pclass + Sex + Embarked,
                           mod03= ~Pclass + Sex + Parch + SibSp + Age,
                           mod04= ~Pclass + Sex + Age + SibSp + Parch + Fare + Embarked)

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 = df_train))) 
model_tidy = models %>% 
  mutate(tidy = map(mod,tidy)) %>% 
  unnest(tidy) %>% 
  mutate(estimate=round(estimate,4),
         p.value=round(p.value,4)) %>%
  mutate(p.adjusted = p.adjust(p.value))%>%
  filter(term != '(Intercept)')%>%
  select(-c(logit_formulas, mod))

#creo un dataset con los terminos no significantes
no_significante = model_tidy %>%
  mutate(significante = p.adjusted < 0.05)%>%#se aplica el p valor ajustado por la cantidad de modelos que se realizan.
  filter(significante == FALSE)%>%
  select(models, expression, term, significante)

no_significante
## # A tibble: 11 x 4
##    models expression                                        term    significante
##    <chr>  <chr>                                             <chr>   <lgl>       
##  1 mod02  Survived ~ Pclass + Sex + Embarked                Pclass2 FALSE       
##  2 mod02  Survived ~ Pclass + Sex + Embarked                Embark~ FALSE       
##  3 mod02  Survived ~ Pclass + Sex + Embarked                Embark~ FALSE       
##  4 mod03  Survived ~ Pclass + Sex + Parch + SibSp + Age     Parch   FALSE       
##  5 mod03  Survived ~ Pclass + Sex + Parch + SibSp + Age     SibSp   FALSE       
##  6 mod04  Survived ~ Pclass + Sex + Age + SibSp + Parch + ~ Pclass2 FALSE       
##  7 mod04  Survived ~ Pclass + Sex + Age + SibSp + Parch + ~ SibSp   FALSE       
##  8 mod04  Survived ~ Pclass + Sex + Age + SibSp + Parch + ~ Parch   FALSE       
##  9 mod04  Survived ~ Pclass + Sex + Age + SibSp + Parch + ~ Fare    FALSE       
## 10 mod04  Survived ~ Pclass + Sex + Age + SibSp + Parch + ~ Embark~ FALSE       
## 11 mod04  Survived ~ Pclass + Sex + Age + SibSp + Parch + ~ Embark~ FALSE

Cuando se agregan más variables al modelo, empiezan a aparecer niveles y variables no significantes. En el modelo 02, los niveles según el sitio de embarcación no son significantes respecto al sitio de embarca basal. Tampoco logra encontrar una diferencia significativa en la Clase 02 respecto a la Clase 01.

El modelo 03 es igual al modelo 01 pero añade las variables Parch y SibSp que resultan no significativas.

El modelo 04 contiene todas las variables pero no encuentra significancia con las variables SibSP, Parch, Fare y con los niveles EmbarkedQ y Pclass02.

#ploteo los coeficientes para ver su variación según los modelos.
ggplot(model_tidy, aes(term, estimate, colour = as.factor(models)))+
  geom_jitter(height=0, width=0.2)+ 
  labs(title ="Variación de coeficientes", x = "Variable", y = "Coeficiente")+
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 10))+
  theme(legend.key = element_blank())

#ploteo la significancia de los coeficientes
ggplot(model_tidy, aes(term, p.adjusted, colour = as.factor(expression)))+
  geom_jitter(height=0, width=0.3)+ 
  labs(title ="Significancia de los coeficientes", x = "Variable", y = "P valor ajustado")+
  geom_hline(yintercept=0.05, linetype="dashed", color = "red")+
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 10))+
  theme(legend.key = element_blank())

El gráfico de significancia de los coeficientes muestra como en los modelos explorados, existen coeficientes que siempre son significativos y otros que no. El único coeficiente que pierde significatividad en presencia de otras variables es el PClass 2.

Respecto a la variación de los coeficientes se puede observar una gran variación en el coeficiente PClass2 y PClass3, los otros coeficientes se mantienen en rangos similares.

  1. Ordenar por la deviance los modelos generados
# Calcular las medidas de evaluación para cada modelo
models <- models %>% 
  mutate(glance = map(mod,glance))

# Obtener las medidas de evaluacion de interes
models %>% 
  unnest(glance) %>%
  # Calculo de la deviance explicada
  mutate(perc_explained_dev = 1-deviance/null.deviance) %>% 
  select(-c(models, df.null, AIC, BIC, mod, logit_formulas)) %>% 
  arrange(deviance)
## # A tibble: 4 x 6
##   expression          null.deviance logLik deviance df.residual perc_explained_~
##   <chr>                       <dbl>  <dbl>    <dbl>       <int>            <dbl>
## 1 Survived ~ Pclass ~          828.  -270.     540.         613            0.348
## 2 Survived ~ Pclass ~          828.  -271.     543.         616            0.344
## 3 Survived ~ Pclass ~          828.  -277.     554.         618            0.331
## 4 Survived ~ Pclass ~          828.  -283.     567.         617            0.315

Según el TP se elige el modelo con la mejor deviance. En este caso corresponde al modelo 04.

4. Evaluación del modelo

  1. Realizar el gráfico de curva ROC y obtener el AUC para el modelo elegido. Interpretar el gráfico
# Añadir las predicciones
models <- models %>% 
  mutate(pred= map(mod,augment, type.predict = "response"))

# Modelo completo
pred_mod <- models %>% 
  filter(models=="mod04") %>% 
  unnest(pred)


# Calculamos curvas ROC
roc_mod <- roc(response=pred_mod$Survived, predictor=pred_mod$.fitted)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
print(paste('AUC: Modelo', roc_mod$auc))
## [1] "AUC: Modelo 0.865082748518834"
#grafico curva roc
ggroc(list(Modelo_4=roc_mod), size=1) + geom_abline(slope = 1, intercept = 1, linetype='dashed') + labs(title='Curva ROC', color='Modelo')

El gráfico de la curva ROC muestra como el modelo se aleja de la diagonal punteada del azar. El área bajo la curva es de 0.865 cuando el AUC del azar es 0.5 y el ajuste perfecto un AUC = 1.

  1. Realizar un violin plot e interpretar
#violin plot
violin_full= ggplot(pred_mod, aes(x=Survived, y=.fitted, group=Survived,fill=factor(Survived))) + 
  geom_violin() +
  guides(fill=FALSE) +
  labs(title='Violin plot', subtitle='Modelo 04', y='Predicted probability')

print(violin_full)

En el violin plot se puede ver como el modelo predice de mejor manera aquellos pasajeros que no sobrevivieron. Si se elige un punto de corte cerca de 0.25, se predicirá relativamente bien la mayoría de los no sobrevivientes. En cuanto a los sobrevivientes, estos presentan una distribución más homogénea en en las diferentes probabilidades y como es de esperar se ensancha en los valores altos de probabilidada.

5. Elección del punto corte

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

#modelo entrenado sobre los datos de train sin validation
mod04_train<- glm(logit_formulas$mod04, family = 'binomial', data = df_train)

pred_mod_val = augment(x=mod04_train, newdata=df_val, type.predict='response') 



cutoff = seq(0.01,0.95,0.01)
#pred_mod_val %>% 
    # mutate(predicted_class=if_else(.fitted>cutoff, 1, 0) %>% as.factor(),
    #        survived= factor(Survived))
#se crea una función de que genere las matrices de confusión para múltiples puntos de corte.
prediction_metrics <- function(cutoff, predictions=pred_mod_val){
  table <- predictions %>% 
    mutate(predicted_class=if_else(.fitted>cutoff, 1, 0) %>% as.factor(),
           survived= factor(Survived))
  
  confusionMatrix(table(table$predicted_class, table$survived), positive = "1") %>%
    tidy() %>%
    select(term, estimate) %>%
    filter(term %in% c('accuracy', 'sensitivity', 'specificity', 'precision','recall')) %>%
    mutate(cutoff=cutoff)
  
}
#se crea el dataset con la función anterior para graficar las métricas

logit_pred= map_dfr(cutoff, 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 04 en validation set', color="")

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

En este caso las clases tienen una proporción en el dataset aproximada del 60% (no sobrevive) y 40% (sobrevive). Por estos valores se puede decir que las clases no tienen un problema de desbalanceo. Dependiendo de la aplicación puede convenir utilizar la sensibilidad o la especificidad u otra métrica. Puesto que la clase es balanceada y no hay preferencia por la predicción de algún tipo específico de clase, se utiliza un valor de la probabilidad de corte de 0,5 que coincide con el máximo valor del accuracy.

  1. Obtener la matriz de confusión con el modelo y punto de corte elegidos. Interpretarla
sel_cutoff = 0.5

# Agregamos la predicciones al dataset de validacion
table= augment(x=mod04_train, newdata=df_val, type.predict='response') 
# Clasificamos utilizamos el punto de corte
table=table %>% mutate(predicted_class=if_else(.fitted>sel_cutoff, 1, 0) %>% as.factor(),
           survived= factor(Survived))

# Creamos la matriz de confusión
confusionMatrix(table(table$survived, table$predicted_class), positive = "1")
## Confusion Matrix and Statistics
## 
##    
##       0   1
##   0 148  15
##   1  34  71
##                                           
##                Accuracy : 0.8172          
##                  95% CI : (0.7656, 0.8616)
##     No Information Rate : 0.6791          
##     P-Value [Acc > NIR] : 2.736e-07       
##                                           
##                   Kappa : 0.6036          
##                                           
##  Mcnemar's Test P-Value : 0.01013         
##                                           
##             Sensitivity : 0.8256          
##             Specificity : 0.8132          
##          Pos Pred Value : 0.6762          
##          Neg Pred Value : 0.9080          
##              Prevalence : 0.3209          
##          Detection Rate : 0.2649          
##    Detection Prevalence : 0.3918          
##       Balanced Accuracy : 0.8194          
##                                           
##        'Positive' Class : 1               
##

La matriz de confusión indica que de los 182 pasajeros que murieron, se logró clasificar correctamente 148 y de los 86 sobrevivientes 71 fueron correctamente clasificados. Se observa como las predicciones eróneas por clase son aproximadamente 1/5 de las correctas en ambos casos. Esto resulta de esta manera por elegir como métrica de punto de corte el accuracy.

6. Dataset de testeo

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

df_test <- read_csv("D:/OneDrive/Ciencia de Datos/06 Enfoque estadistico/TP01/TP03/titanic_complete_test.csv")
df_test = select(df_test, PassengerId, Survived, Pclass, Sex, Age, SibSp,Parch, Fare, Embarked)
df_test[charcol] <- lapply(df_test[charcol], factor)

Con el modelo y punto de corte elegidos clasificar a las personas del dataset de testing. Obtener la matriz de confusión y comparar con la obtenida en el punto 5)c)

#entreno con todos los datos del train (validation y train)
mod04_train_val<- glm(logit_formulas$mod04, family = 'binomial', data = df)

# Agregamos la predicciones al dataset de test
table= augment(x=mod04_train_val, newdata=df_test, type.predict='response') 
# Clasificamos utilizamos el punto de corte
table=table %>% mutate(predicted_class=if_else(.fitted>sel_cutoff, 1, 0) %>% as.factor(),
           survived= factor(Survived))
# Creamos la matriz de confusión
confusionMatrix(table(table$survived, table$predicted_class), positive = "1")
## Confusion Matrix and Statistics
## 
##    
##       0   1
##   0 215  46
##   1  48 109
##                                          
##                Accuracy : 0.7751         
##                  95% CI : (0.732, 0.8143)
##     No Information Rate : 0.6292         
##     P-Value [Acc > NIR] : 1.055e-10      
##                                          
##                   Kappa : 0.5193         
##                                          
##  Mcnemar's Test P-Value : 0.9179         
##                                          
##             Sensitivity : 0.7032         
##             Specificity : 0.8175         
##          Pos Pred Value : 0.6943         
##          Neg Pred Value : 0.8238         
##              Prevalence : 0.3708         
##          Detection Rate : 0.2608         
##    Detection Prevalence : 0.3756         
##       Balanced Accuracy : 0.7604         
##                                          
##        'Positive' Class : 1              
##

En este caso, se puede observar como el accuracy que es la métrica elegida, pasó del 82% al 78% al utilizar los datos del test. Resulta interesante notar que el Sensitivity disminuyó de 83% a 70% pero el Specificity se mantuvo práctimente igual. Esto quiere decir que en el test, el modelo predice mejor los no sobrevivientes que los sobrevivientes en comparación con los datos de validación.