library(tidyverse)
library(GGally)
library(tidymodels)
library(glmnet)
library(modelr)
library(pROC)
set.seed(2048)
'data.frame': 418 obs. of 12 variables:
$ PassengerId: int 892 893 894 895 896 897 898 899 900 901 ...
$ Pclass : int 3 3 2 3 3 3 3 2 3 3 ...
$ Name : Factor w/ 418 levels "Abbott, Master. Eugene Joseph",..: 210 409 273 414 182 370 85 58 5 104 ...
$ Sex : Factor w/ 2 levels "female","male": 2 1 2 2 1 2 1 2 1 2 ...
$ Age : num 34.5 47 62 27 22 14 30 26 18 21 ...
$ SibSp : int 0 1 0 0 1 0 0 1 0 2 ...
Error in gregexpr(calltext, singleline, fixed = TRUE) :
regular expression is invalid UTF-8
$ Parch : int 0 0 0 0 1 0 0 1 0 0 ...
$ Ticket : Factor w/ 363 levels "110469","110489",..: 153 222 74 148 139 262 159 85 101 270 ...
$ Fare : num 7.83 7 9.69 8.66 12.29 ...
$ Cabin : Factor w/ 76 levels "A11","A18","A21",..: NA NA NA NA NA NA NA NA NA NA ...
$ Embarked : Factor w/ 3 levels "C","Q","S": 2 3 2 3 3 3 2 3 1 3 ...
$ Survived : int 0 1 0 0 1 1 0 1 1 0 ...
El dataset Titanic
tiene 418 observaciones y 12 variables. Resulta conveniente tratas como factor a la variables Survived y Pclass.
Se seleccionan las variables PassengerId, Survived, Pclass, Sex, Age, SibSp, Parch, Fare y Embarked, y se transforman a factor las variables Survived, Pclass y Embarked.
En el siguiente gráfico vemos la correlación y las distribuciones de las variables Survived, Pclass, Sex, Age y Fare. En color rojo se destacan los que no sobrevivieron al hundimiento y en color celeste a aquellos que sí. Algunos puntos a destacar:
Calculamos la distribución de clase para el dataset completo.
train %>%
group_by(Survived) %>%
summarise(
numero_casos = n(),
porcentaje = numero_casos / nrow(train)
)
ratio_supervivencia <- function(data) {
tmp <- data %>% filter(Survived==1) %>% dplyr::select(1) %>% nrow()
tmp2 <- data %>% dplyr::select(1) %>% nrow()
paste("sobrevive el", round(tmp/tmp2,3)*100,"%")
}
train <- train %>% mutate(tag = case_when(Survived == 1 ~ "Sobrevive",
Survived == 0 ~ "No sobrevive"))
ggplot(train)+ geom_bar( aes(x = Survived, fill=tag))+
labs(title = "Titanic. Sobrevivientes y no sobrevivientes",
subtitle = "Fuenet: Kaggle") +
labs(x = "Sobreviviente", y = "# personas")+
theme(legend.position = "bottom", legend.direction = "horizontal", legend.title = element_blank()) + guides(fill = guide_legend(reverse = T))
[1] "sobrevive el 37.6 %"
Dividimos el dataset en entrenamiento (70%) y validación (30%)
#Separo en train y test
set.seed(2048)
titanic_split <- initial_split(train, prop = 0.7)
titanic_split
<293/125/418>
Verificamos la distribución de clase nuevamente.
#TRAINING
titanic_split %>%
training() %>%
group_by(Survived) %>%
summarise(
numero_casos = n(),
porcentaje = numero_casos / nrow(training(titanic_split))
)
[1] "sobrevive el 38.9 %" "en train"
#VALIDACION
#TRAINING
titanic_split %>%
testing() %>%
group_by(Survived) %>%
summarise(
numero_casos = n(),
porcentaje = numero_casos / nrow(training(titanic_split))
)
[1] "sobrevive el 34.4 %" "en validación"
El dataset de validacion tiene 39.1% de supervivientes y validacion 34.4%. Para el tamaño del dataset (418 obs.) es un balanceo adecuado.
\[Survived_i = \beta_0 + \beta1 Pclass_i + \beta_2 Sex_i + \beta_3 Age_i \]
mod_0 <- glm(Survived ~ Pclass + Sex + Age, data = training(titanic_split), family = "binomial")
summary(mod_0)
Call:
glm(formula = Survived ~ Pclass + Sex + Age, family = "binomial",
data = training(titanic_split))
Deviance Residuals:
Min 1Q Median 3Q Max
-2.0901 -0.6147 -0.5156 0.7591 2.2017
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.30916 0.65893 5.022 5.11e-07 ***
Pclass2 -1.67309 0.46407 -3.605 0.000312 ***
Pclass3 -2.00682 0.41455 -4.841 1.29e-06 ***
Sexmale -2.28525 0.29939 -7.633 2.29e-14 ***
Age -0.03457 0.01370 -2.523 0.011625 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 391.64 on 292 degrees of freedom
Residual deviance: 288.00 on 288 degrees of freedom
AIC: 298
Number of Fisher Scoring iterations: 4
Pclass2, Pclass3: coeficientes negativos y significativos al 1%. Es decir que pertenecer a las clases 2 y 3 reduce la probabilidad de sobrevivir.
Sexomale es negativo y significativo a 1%: ser hombre reduce la probabilidad de sobrevivir.
Age es negativo y significativo al 5%: mayor edad reduce la probabilidad de sobrevivir.
newdata = data.frame("Pclass"=c("1","3"),"Sex"=c("female","male"),"Age"=c(17,20))
predict.glm(object = mod_0, newdata,type = "response")
1 2
0.9382844 0.1578577
Rose tiene más probabilidades de sobrevivir (93.8%), mientras que Jack tiene 15.7% de probabilidades de sobrevivir.
Se generaron tres modelos:
\[Survived_i = \beta_0 + \beta1 Pclass_i + \beta_2 Sex_i + \beta_3 Age_i + \beta_4 Embarked_i\]
\[Survived_i = \beta_0 + \beta1 Pclass_i + \beta_2 Sex_i + \beta_3 Age_i + \beta_4 Embarked_i+ \beta_5 Parche_i\]
\[Survived_i = \beta_0 + \beta_1 Sex_i + \beta_2 Fare_i + \beta_3 Embarked_i\] ## 3.b. Selección del mejor modelo por deviance explicada
#Genero formulas de modelo 2_a y tres nuevos modelos
logit_formulas <- formulas(.response = ~Survived,
mod_2a= ~ Pclass + Sex + Age,
mod_1= ~ Pclass + Sex + Age + Embarked,
mod_2= ~Pclass + Sex + Age + Parch + Embarked,
mod_3= ~ Sex + Fare + Embarked
)
#corro lo modelos logit
models <- tibble(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 = training(titanic_split))))
#agrego estadisticas resumen de los modelos
models <- models %>%
mutate(glance = map(mod,glance))
#tabla salida ordenada por deviance
models %>%
unnest(glance) %>% arrange(deviance)
Como era esperable el modelo con mayor cantidad de regresores es que presenta la menor deviance. mientras que el modelo del punto 2a al ser el más parsimonioso logra captar una menor proporción de la deviance explicada
La curva ROC mide la relación entre la tasa de verdaderos positivos y la tasa de falsos positivos en distintos umbrales de probabilidad. Asimismo, el área bajo la curva ROC (AUC) nos permite comparar la capacidad de clasificación de distintos modelos. En el caso del dataset Titanic, la curva AUB registra un valor e 0.82, por lo que podemos decir que el modelo 2 es un buen modelo para predecir la supervivencia de los pasajeros.
# Prediccion Modelo elegido en TRAIN
prediction_mod_2 <- models %>%
mutate(pred= map(mod,augment,newdata=training(titanic_split), type.predict = "response")) %>%
filter(models=="mod_2") %>%
unnest(pred)
# Calculamos curvas ROC
roc_mod_2 <- roc(response=prediction_mod_2$Survived, predictor=prediction_mod_2$.fitted)
Setting levels: control = 0, case = 1
Setting direction: controls < cases
[1] "AUC: Modelo 2 0.822919729491326"
El violin plot presentado abajo, es un tipo de gráfico que permite observar la distribución de probabilidad predicha para la dos clases del dataset: supervivientes y no supervivientes. Como se puede observar la mayor parte de los que no sobreviven registran probabilidades predichas por debajo de 0.25, mientras que en el grupo de lo supervivientes la mayor parte presenta una probabilidad predicha superior al 50%.
Este gráfico esta estrechamente vinculado al gráfico de la curva ROC presentado arriba, ya que de definirse un umbral de probabilidad predicha más alto se obtendrá un menor recall y un mayor sensitivity (o TPR) a costa de mayor specificity (o FPR).
Sobre el dataset de validación se calcularon las curvas de acuracy, Specificity, Recall y Precision en función del punto de corte. En el punto de corte igual a 0.35 encontramos la intersección entre las curvas acuracy, recall, sensitivity y specificity. Es decir que en 0.35 tendremos un balance entre las distintas métricas y por ello elegimos ese punto como el punto de corte.
library(caret)
library(e1071)
# Prediccion Modelo elegido en VALIDACION
prediction_mod_2_validacion <- models %>%
mutate(pred= map(mod,augment,newdata=testing(titanic_split), type.predict = "response")) %>%
filter(models=="mod_2") %>%
unnest(pred)
prediction_metrics <- function(cutoff, predictions=prediction_mod_2_validacion){
table <- predictions %>%
mutate(predicted_class=if_else(.fitted>cutoff,"1", "0") %>% as.factor(),
Survived= factor(Survived)) %>% select(predicted_class,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)
}
seq(0.10,0.95,0.05)
[1] 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95
cutoffs = seq(0.10,0.95,0.1) #Da error si hago con menor separacion de cutoff
logit_pred= map_dfr(cutoffs, prediction_metrics) %>% mutate(term=as.factor(term))
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, alpha=0.5) +
theme_bw() +
labs(title= 'Accuracy, Sensitivity, Specificity, Recall y Precision', subtitle= 'Modelo completo', color="")
A continuación se presenta la matríz de confusión para el dataset de validación. De los 43 casos postitivos se logra clasificar de forma correcta a 33; mientras que del total de casos no positivos (82) se clasifica correctamente a 61. De esta forma el accuracy es de 75.2% con un intervalo de confianza del 95% entre 0.66 y 0.82, y que el mismo es significativamente (p-valor pequeño) superior que el NIR (No Information Rate). Podemos decir también que predecimos levemente mejor los casos de supervivencia que los de no supervivencia (0.76 vs 0.74).
tabla <- prediction_mod_2_validacion %>%
mutate(predicted_class=if_else(.fitted>0.35,"1", "0") %>% as.factor(),
Survived= factor(Survived)) %>% select(predicted_class,Survived)
confusionMatrix(table(tabla$predicted_class, tabla$Survived), positive = "1")
Confusion Matrix and Statistics
0 1
0 61 10
1 21 33
Accuracy : 0.752
95% CI : (0.6668, 0.8249)
No Information Rate : 0.656
P-Value [Acc > NIR] : 0.01355
Kappa : 0.482
Mcnemar's Test P-Value : 0.07249
Sensitivity : 0.7674
Specificity : 0.7439
Pos Pred Value : 0.6111
Neg Pred Value : 0.8592
Prevalence : 0.3440
Detection Rate : 0.2640
Detection Prevalence : 0.4320
Balanced Accuracy : 0.7557
'Positive' Class : 1
A continuación se presenta la matríz de confusión obtenida con las predicciones del modelo elegido sobre el dataset de testing.
test<-read.csv("titanic_complete_test.csv")
#Transformo las variables Survived, Pclass y Embarked a factor
test <- test %>%
dplyr::select(Survived, Pclass, Sex, Age, SibSp,Parch, Fare,Embarked) %>%
mutate( Survived = as.factor(Survived),
Pclass = as.factor(Pclass),
Embarked = as.factor(Embarked) )
str(test)
'data.frame': 418 obs. of 8 variables:
$ Survived: Factor w/ 2 levels "0","1": 1 2 1 1 2 2 1 2 2 1 ...
$ Pclass : Factor w/ 3 levels "1","2","3": 3 3 2 3 3 3 3 2 3 3 ...
$ Sex : Factor w/ 2 levels "female","male": 2 1 2 2 1 2 1 2 1 2 ...
$ Age : num 34.5 47 62 27 22 14 30 26 18 21 ...
$ SibSp : int 0 1 0 0 1 0 0 1 0 2 ...
$ Parch : int 0 0 0 0 1 0 0 1 0 0 ...
$ Fare : num 7.83 7 9.69 8.66 12.29 ...
$ Embarked: Factor w/ 3 levels "C","Q","S": 2 3 2 3 3 3 2 3 1 3 ...
Para el caso de los supervivientes, se logra clasificar correctamentea 119 de 157 asos; mientras que de los no spervivientes se clasifica correctamente a 195 de 261 casos. En suma, se obtuvo un accuracy de 75.12% levemente por debajo del 75.2% obtenido en validación.
Confusion Matrix and Statistics
0 1
0 195 38
1 66 119
Accuracy : 0.7512
95% CI : (0.7069, 0.7919)
No Information Rate : 0.6244
P-Value [Acc > NIR] : 2.419e-08
Kappa : 0.4878
Mcnemar's Test P-Value : 0.008107
Sensitivity : 0.7580
Specificity : 0.7471
Pos Pred Value : 0.6432
Neg Pred Value : 0.8369
Prevalence : 0.3756
Detection Rate : 0.2847
Detection Prevalence : 0.4426
Balanced Accuracy : 0.7525
'Positive' Class : 1