TP 3 - Carlos Arana
Preparación de los datos
- Leer el archivo titanic_complete_train.csv y mostrar su estructura
library(readr)
library(tidyverse)## -- Attaching packages ---------------------------------------------------- tidyverse 1.2.1 --## v ggplot2 3.2.1 v purrr 0.3.3
## v tibble 2.1.3 v dplyr 0.8.3
## v tidyr 1.0.0 v stringr 1.4.0
## v ggplot2 3.2.1 v forcats 0.4.0## -- Conflicts ------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()ds <- read.csv("C:/Users/charly/Desktop/FCEyN- MAESTRIA/Materia - ENFOQUE ESTADISTICO DEL APRENDIZAJE/trabajos_practicos/TP-3/titanic_complete_train.csv")
glimpse(ds)## Observations: 891
## Variables: 12
## $ PassengerId <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,...
## $ Survived <int> 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0,...
## $ Pclass <int> 3, 1, 3, 1, 3, 3, 1, 3, 3, 2, 3, 1, 3, 3, 3, 2, 3,...
## $ Name <fct> "Braund, Mr. Owen Harris", "Cumings, Mrs. John Bra...
## $ Sex <fct> male, female, female, female, male, male, male, ma...
## $ Age <dbl> 22.00000, 38.00000, 26.00000, 35.00000, 35.00000, ...
## $ SibSp <int> 1, 1, 0, 1, 0, 0, 0, 3, 0, 1, 1, 0, 0, 1, 0, 0, 4,...
## $ Parch <int> 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 1, 0, 0, 5, 0, 0, 1,...
## $ Ticket <fct> A/5 21171, PC 17599, STON/O2. 3101282, 113803, 373...
## $ Fare <dbl> 7.2500, 71.2833, 7.9250, 53.1000, 8.0500, 8.4583, ...
## $ Cabin <fct> NA, C85, NA, C123, NA, NA, E46, NA, NA, NA, G6, C1...
## $ Embarked <fct> S, C, S, S, S, Q, S, S, S, C, S, S, S, S, S, S, Q,...- Seleccionar las variables PassengerId, Survived, Pclass, Sex, Age, SibSp,Parch, Fare y Embarked
ds = ds %>% select(PassengerId, Survived, Pclass, Sex, Age, SibSp,Parch, Fare, Embarked)- Transformar las variables Survived, Pclass y Embarked a factor
ds[,c("Survived", "Pclass", "Embarked")] = lapply(ds[,c("Survived", "Pclass", "Embarked")], factor)
str(ds)## 'data.frame': 891 obs. of 9 variables:
## $ PassengerId: int 1 2 3 4 5 6 7 8 9 10 ...
## $ Survived : Factor w/ 2 levels "0","1": 1 2 2 2 1 1 1 1 2 2 ...
## $ Pclass : Factor w/ 3 levels "1","2","3": 3 1 3 1 3 3 1 3 3 2 ...
## $ Sex : Factor w/ 2 levels "female","male": 2 1 1 1 2 2 2 2 1 1 ...
## $ Age : num 22 38 26 35 35 ...
## $ SibSp : int 1 1 0 1 0 0 0 3 0 1 ...
## $ Parch : int 0 0 0 0 0 0 0 1 2 0 ...
## $ Fare : num 7.25 71.28 7.92 53.1 8.05 ...
## $ Embarked : Factor w/ 3 levels "C","Q","S": 3 1 3 3 3 2 3 3 3 1 ...- Realizar un gráfico de ggpairs para las variables Survived, Pclass, Sex, Age y Fare e interpretarlo
library(GGally)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2##
## Attaching package: 'GGally'## The following object is masked from 'package:dplyr':
##
## nasads %>%
select(Survived, Pclass, Sex, Age, Fare) %>%
ggpairs(.,mapping = aes(colour= Survived)) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + theme_bw()## `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 la luz del correlograma procederemos a analizar “grosso modo” las relaciones entre la variable Survive y las restantes variables explicativas, a saber:
-Pclass: la variable Pclass y Survived no son independientes, pues es notable la diferencia entre la probabilidad condicional de Sobreviviente dado clase que de sobreviviente (marginal). Esta desigualdad matemática implica que no se mantiene la porción de sobrevivientes-No sobrevivientes en cada una de las clases de Class, notándose a simple vista que hay una mayor tendencia a no sobrevivir si el viajero se encuentra en la clase más baja (Pclass 3). La distribución marginal de Survived se invierte en los caso de Clases 1 y 2, siendo mayoritaria la clase “Sobreviviente”
-Género: misma situación matemática que la anterior, al no haber igualdad entra la condicional y la marginal, implicando falta total de independencia. Esta situación se nota claramente al observar la notable diferencia entre las distribuciones condicionales en las que se aprecias que la tasa de supervivencia en la clase “Femenino” de Sex invierte la distribución de la marginal “Survived”, o sea, se salvan en muchísima mayor proporción las mujeres que los hombres
-Age: aquí esperaríamos observar lo mismo que en el caso anterior, haciendo honor a la tradición de salvataje en siniestros y tragedias que reza “Mujeres y niños primero”. Pero no fue así: no se nota diferencia sustancial entre las distribuciones condicionales de Survived condicionada a edad vs la distrinucion marginal de Surviced, por lo que no se evidencia que salvar a los niños haya sido una preferencia.
-Fare: aun habiendo alguna diferencia entre las condicionaels y la marginal no pareceiera ser a simple vista sustancial y/o siginificativa. Deberñia someter esta relación a un análisis más echaustivo.
- Mostrar la distribución de clase (Sobrevivientes vs No Sobrevivientes)
table(ds$Survived)/nrow(ds)##
## 0 1
## 0.6161616 0.3838384ds %>%
select(Survived) %>%
ggpairs(., mapping = aes(color = Survived))
- 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)
set.seed(111)
muestra=sample(nrow(ds), 0.7*nrow(ds))
train=ds[muestra,]
validation=ds[-muestra,]
# distribución de clase en train_test
table(train$Survived)/nrow(train)##
## 0 1
## 0.6211878 0.3788122# distribución de clase en test_test
table(validation$Survived)/nrow(validation)##
## 0 1
## 0.6044776 0.3955224Podemos observar que la distribución se mantiene, pues en ambos casos (entrenamiento y validación) rondan las proporciones .60 - .40 para las clase 0 y 1 (No Sobreviviente-Sobreviviente) respectivamentede la Variable objetivo “Survived”
———————————- FIN PUNTO 1 —————————————-
- Predicciones (Trabajar con dataset de ENTRENAMIENTO)
- 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
library(broom)
library(stats)
rl <- glm(formula = Survived ~ Pclass+Sex+Age, family= "binomial", data = train)
rl_tdy <- rl %>% tidy()
rl_tdy- Dar una breve interpretación de los coeficientes y su significatividad
A la luz de evidencia aportada por los datos de este dataset todas las variables resultaron ser significativas a la hora de explicar la probabilidad de que un pasajero sobreviva o no mediante un modelo de regresión logística
El signo y los valores de los coeficientes que arroja el ajuste efectuado por el algoritmo glm iindican que:
-EL cambio de clase 1 a 2 y a 2 implica en ambos casos (evidenciado por los coeficientes asociados a las variables dummies Pclass 2 y Pclass3) una disminución en la probabilidad de sobrevivir.
-El ser hombre implica una disfunción de la probabilidad de supervivencia respecto al nivel basal de esta variable, que es “Female”
-El incremento de La edad, también implica, aunque en un una cuantía bastante menor, una disminución en la probabilidad de supervivencia
- ¿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
str(ds$Sex)## Factor w/ 2 levels "female","male": 2 1 1 1 2 2 2 2 1 1 ...pr_Rose = predict(rl,type = 'response', newdata =data.frame(Pclass="1", Sex="female",Age=17))
pr_Jack = predict(rl,type = 'response', newdata =data.frame(Pclass="3", Sex="male",Age=20))
data.frame(Pasajero=c("Rose","Jack"), Prob_Sobrevive=c(pr_Rose , pr_Jack))Aquí vemos claramente todo lo indicado previamente respecto a la distrubuión de la variable Survived condicinada al sexo. Aún estabdo condicionado también por la edad, la probabilidad de que Rose sobreviva es mucho mayor a la de que Jack sobreviva pues P(Survived/Female) >> P(P(Survived/male)
———————————- FIN PUNTO 2 —————————————-
Generación de modelos (Trabajar con dataset de ENTRENAMIENTO)
- 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
library(modelr)##
## Attaching package: 'modelr'## The following object is masked from 'package:broom':
##
## bootstraplogit_formulas <- formulas(.response = ~Survived,
Pclassbal= ~Pclass,
Sex_class= ~ Pclass + Sex,
Sex_class_Age_Fare= ~ Sex + Pclass + Age + Fare
)Generé tres modelos cuyas formulas son : Survived ~ Pclass, Survived ~ Pclass + Sex,
Survived ~ Sex + Pclass + Age + Fare
- 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
logit_formulas_b = formulas(.response = ~Survived,
Pclass= ~Pclass,
Sex_class= ~ Pclass + Sex,
Sex_class_Age = ~ Sex + Pclass + Age,
Sex_class_Age_Fare= ~ Sex + Pclass + Age + Fare
)
models <- data_frame(logit_formulas_b) %>% # dataframe a partir del objeto formulas
mutate(models = names(logit_formulas_b), # columna con los nombres de las formulas
expression = paste(logit_formulas_b), # columna con las expresiones de las formulas
mod = map(logit_formulas_b, ~glm(.,family = 'binomial', data = train))) # corremos cada modelo de 'fornulas' y lo cargamos su resultado como una lista en la columna mod## Warning: `data_frame()` is deprecated, use `tibble()`.
## This warning is displayed once per session.models <- models %>%
mutate(glance = map(mod,glance))
models %>%
unnest(glance, .drop = TRUE) %>%
mutate(perc_explained_dev = 1-deviance/null.deviance) %>%
select(-c(models, df.null, AIC, BIC)) %>%
arrange(deviance)## Warning: 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.El mejor mopdelo es Sex_class_Age_Fare, cuya fórmula es: Survived ~ Sex + Pclass + Age + Fare
———————————- FIN PUNTO 3 —————————————-
- Evaluación del modelo (Trabajar con dataset de ENTRENAMIENTO)
- Realizar el gráfico de curva ROC y obtener el AUC para el modelo elegido. Interpretar el gráfico
library(pROC)## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, varmodels <- models %>%
mutate(pred= map(mod,augment, type.predict = "response"))
models$pred[1]## $Pclass
## # A tibble: 623 x 10
## .rownames Survived Pclass .fitted .se.fit .resid .hat .sigma .cooksd
## <chr> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 507 1 2 0.496 0.0455 1.18 0.00826 1.11 2.85e-3
## 2 724 0 2 0.496 0.0455 -1.17 0.00826 1.11 2.76e-3
## 3 175 0 1 0.605 0.0403 -1.36 0.00680 1.11 3.53e-3
## 4 793 0 3 0.245 0.0228 -0.750 0.00282 1.11 3.07e-4
## 5 699 0 1 0.605 0.0403 -1.36 0.00680 1.11 3.53e-3
## 6 69 1 3 0.245 0.0228 1.68 0.00282 1.11 2.91e-3
## 7 675 0 2 0.496 0.0455 -1.17 0.00826 1.11 2.76e-3
## 8 584 0 1 0.605 0.0403 -1.36 0.00680 1.11 3.53e-3
## 9 794 0 1 0.605 0.0403 -1.36 0.00680 1.11 3.53e-3
## 10 305 0 3 0.245 0.0228 -0.750 0.00282 1.11 3.07e-4
## # ... with 613 more rows, and 1 more variable: .std.resid <dbl>prediction_Sex_class_Age_Fare <- models %>%
filter(models=="Sex_class_Age_Fare") %>%
unnest(pred, .drop=TRUE)
roc_full <- roc(response=prediction_Sex_class_Age_Fare$Survived, predictor=prediction_Sex_class_Age_Fare$.fitted)## Setting levels: control = 0, case = 1## Setting direction: controls < casesprediction_bad <- models %>%
filter(models=="Pclass") %>%
unnest(pred, .drop=TRUE)
roc_bad <- roc(response=prediction_bad$Survived, predictor=prediction_bad$.fitted)## Setting levels: control = 0, case = 1
## Setting direction: controls < casesggroc(list(Sex_class_Age_Fare=roc_full, Pclass=roc_bad), size=1) + geom_abline(slope = 1, intercept = 1, linetype='dashed') + theme_bw() + labs(title='Curvas ROC', color='Modelo')
Esencialmente l curva ROC mide, para modelos de clasificación con variable objetivo dicotómica, la performance conjunta del modelo en términos de capacidad de clasificar a los casos positivos y a los negativos correctamente (sensitividad y especificidad respectivamente). La forma de medir esta performance es mediante la AUC (are andar ROC Curve) que mide el área debajo de la curva ROC. Como la curva que con mejor performance es la que más se acerca a los ejes izquierdo y superior del gráfico, la AUC asociada a ella tenderá a tomar un valor de 1. SI el modelo no tuviera capacidad predictiva y eligiera la clase la variable objetivo al azar, la curva ROC tenderá a ser la bisectriz del gráfico, y por ende su AUC igual a 0.5
Como podemos notar en nuestro gráfico ROC-AUC el modelo “Sex + Pclass + Age + Fare” es el que mayor valor de AUC posee.
- Realizar un violin plot e interpretar
library(cowplot)##
## ********************************************************## Note: As of version 1.0.0, cowplot does not change the## default ggplot2 theme anymore. To recover the previous## behavior, execute:
## theme_set(theme_cowplot())## ********************************************************violin_full=ggplot(prediction_Sex_class_Age_Fare, 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')
violin_bad=ggplot(prediction_bad, aes(x=Survived, y=.fitted, group=Survived, fill=factor(Survived))) +
geom_violin() +
theme_bw() +
guides(fill=FALSE) +
labs(title='Violin plot', subtitle='Modelo malo', y='Predicted probability')
plot_grid(violin_bad, violin_full)
Aquí podemos ver que el modelo ajusta bien a los datos y generaliza bien en nuevos datos, pues para cada clase el cuerpo del violín se encuentro invertido, indicando que el modelo que está realizando predicciones correctas.
———————————- FIN PUNTO 4 —————————————-
- Elección del punto corte (Trabajar con dataset de VALIDACION)
- Sobre el dataset de validación realizar un gráfico de Accuracy, Specificity, Recall y Precision en función del punto de corte
# predicciones del mejor modelo en el set de validación
models_val <- models %>%
filter(models == 'Sex_class_Age_Fare') %>%
mutate(pred = map(mod, augment, newdata = validation, type.predict = "response"))# Obtenemos las correspondientes a nuestro mejor modelo
prediction_validation_Sex_class_Age_Fare <-
models_val %>%
filter(models == "Sex_class_Age_Fare") %>%
unnest(pred, .drop = TRUE)library(caret)## Loading required package: lattice##
## Attaching package: 'caret'## The following object is masked from 'package:purrr':
##
## liftprediction_metrics <- function(cutoff, predictions=prediction_validation_Sex_class_Age_Fare ){
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)
}
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 completo', color="")
Al no tener específicamente indicada qué métrica se quiere maximizar vamos a optar por un umbral de clasificación que iguala la sensitividad con la especificidad, obteniendo así un balance entre los casos clasificados correctamente como Survived y no Survived Este umbral podemos ver que se encuentra aproximadamente en el punto 0.45
- Obtener la matriz de confusión con el modelo y punto de corte elegidos. Interpretarla
rl_optimo <- glm(formula = Survived ~ Pclass+Sex+Age+Fare, family= "binomial", data = train)
# Agrego las predicciones al dataset de validación
table = augment(x = rl_optimo,
newdata = validation,
type.predict = 'response')
# Clasifico utilizando el umbral elegido
table = table %>% mutate(
predicted_class = if_else(.fitted >= 0.45, 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 129 21
## 1 33 85
##
## Accuracy : 0.7985
## 95% CI : (0.7454, 0.8449)
## No Information Rate : 0.6045
## P-Value [Acc > NIR] : 8.598e-12
##
## Kappa : 0.5867
##
## Mcnemar's Test P-Value : 0.1344
##
## Sensitivity : 0.8019
## Specificity : 0.7963
## Pos Pred Value : 0.7203
## Neg Pred Value : 0.8600
## Prevalence : 0.3955
## Detection Rate : 0.3172
## Detection Prevalence : 0.4403
## Balanced Accuracy : 0.7991
##
## 'Positive' Class : 1
## Mediante la aplicación de un modelos de regresión logística a los datos observados de entrenamiento hemos podido conseguir predecir en forma correcta si un pasajero del Titanic sobreviviría o no conociendo las variables explicativos “Sexo” “Clase” “Edad” y “Tarifa” en el 79.85 % de los casos.
Nuestra capacidad de predecir que un pasajero se va a salvar y que esto realmente ocurra es de un 80.19%, levemente mejor que la capacidad de predecir en forma correcta la suerte del pasajero. Sin embargo la capacidad de predecir que no sobrevivirá y que esto realmente ocurra es levemente menor, siendo su probabilidad de 0.7963. Estas dos proporciones (0.8019 y 0.7963) se denominan sensibilidad y especificidad, y sus formulas son: Sensibilidad: VP/(VP+FN) Especificidad: VN /(VN + FP)
———————————- FIN PUNTO 5 —————————————-
Dataset de testeo (Trabajar con dataset de TESTEO)
- Leer el archivo titanic_complete_test.csv y transformar las variables Survived, Pclass y Embarked a factor
test <- read.csv("C:/Users/charly/Desktop/FCEyN- MAESTRIA/Materia - ENFOQUE ESTADISTICO DEL APRENDIZAJE/trabajos_practicos/TP-3/titanic_complete_test.csv")
test[,c("Survived", "Pclass", "Embarked")] = lapply(test[,c("Survived", "Pclass", "Embarked")], factor)- Con el modelo y punto de corte elegidos clasificar a las personas del dataset de testing.
# Agrego la predicciones al dataset de testeo
table_test = augment(x = rl_optimo,
newdata = test,
type.predict = 'response')
# Clasifico utilizando el punto de corte
table_test = table_test %>% mutate(
predicted_class = if_else(.fitted >= 0.45, 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)- Obtener la matriz de confusión y comparar con la obtenida en el punto 5)c).
confusionMatrix(table_test$predicted_class, table_test$Survived, positive = "1")## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 206 47
## 1 55 110
##
## Accuracy : 0.756
## 95% CI : (0.7119, 0.7964)
## No Information Rate : 0.6244
## P-Value [Acc > NIR] : 6.988e-09
##
## Kappa : 0.485
##
## Mcnemar's Test P-Value : 0.4882
##
## Sensitivity : 0.7006
## Specificity : 0.7893
## Pos Pred Value : 0.6667
## Neg Pred Value : 0.8142
## Prevalence : 0.3756
## Detection Rate : 0.2632
## Detection Prevalence : 0.3947
## Balanced Accuracy : 0.7450
##
## 'Positive' Class : 1
## Aquí podemos observar leve variación de la accuracy, lo cuál es lógico pues esta muestra no fue utilizada ni para entrenar el modelo (dataset “train”) ni para delimitar el punto de corte óptimo (dataset “validation”).
Las interpretaciones de las métricas más importantes, sensibilidad y especificidad, son las mismas, y los valores que toman en este test de prueba son:
Sensitivity : 0.7070
Specificity : 0.7893
Como podemos observar la performance, evaluada en términos de estas dos métricas, se invierte (ahora es mayor la capacidad de predecir “No sobrevivientes” que “Sobrevivientes”), efecto debido a la misma causa que en el caso de la métrica “accuracy”, la mera aleatoriedad a la hora de seleccionar el dataset “test”.
———————————- FIN TP —————————————-