TP3
Edgar Sánchez
1 de diciembre de 2019
Regresión lógistica
En este trabajo se usan datos del Titanic, los cuales tienen información sobre los pasajeros, con una características de cada uno y un target por si sobrevivieron o no.
Librerías
library(tidyverse)
library(dplyr)
library(tidyr)
library(ggplot2)
library(skimr)
library(broom)
library(OneR)
library(rlang)
library(caret)
library(gridExtra)
library(GGally)
library(modelr)
library(pROC)
library(cowplot)
library(knitr)
library(kableExtra)
library(rmarkdown)
set.seed(10291)Lectura y preparación de los datos
El archivo es un texto plano separado por comas.
## Leer el archivo contenido en formato csv
titanic <- read_csv("titanic_complete_train.csv")
paged_table(titanic[1:4,]) ##Imprimir primeras filasInicialmente es necesario aplicar algunas trasnformaciones sobre los datos, ya que existen variables categoricas y estas deben ser tratads como factores; por otro lado se eliminan variables que no harán parte del proceso, como el nombre del pasajero, el número de boleto y el numero de cabina.
## Modificar variables como factor
titanic <- titanic %>%
mutate( Survived = factor(Survived),
Pclass = factor(Pclass ),
Embarked = factor(Embarked )) %>%
select( -c( Name, Ticket, Cabin))A contiuación se presentan gráficos, con el fin de observar el comportamiento de los datos, algunos puntos que llaman la atención:
- No existe correlación lineal entre las dos variables continuas.
- se observa que la mayor parte de las mujeres sobrevivieron; caso contrario en los hombres, donde la mayor parte murió.
- En cuanto a la clase en que viajaban se observa que a medida que disminuye la clase, aumenta la porción de personas que mueren.
- La tarifa de las personas que no sobrevivieron es baja, algo acorde a la clase mencionada anteriormente.
- La purta de emabrque también parece tener una asociación con la cantidad de personas que sobrevivieron.
Hay que destacar que la clase, la tarifa y la puerta de emarque están asociadas de gran manera.
En cuanto a la porción de personas que sobrevivieron, se observa que solo cerca del 40% de las personas sobrevivieron.
ggpairs(data=titanic, ## Grafico
mapping = aes(color = Survived), ## Color dado a la supervivencia
columns=c('Survived', 'Pclass', 'Sex', 'Age','Fare' )) ##Incluir las columnas
ggplot(titanic, aes(Survived)) + ## Grafico de la variable survived
geom_bar(aes(fill = Survived)) ## Color dado a la supervivencia
Para evaluar el modelo se opta por dividir el dataset de entrenamiento en un 70% para ajusatar el modelo, y un 30% para validación del mismo. En ambas particiones se observa que la proporción de la clase permanece al mismo nivel en el que se encontraba inicialmente.
## Dividir data en 70-30
splited_data <- titanic %>% resample_partition(c(train=0.7,test=0.3))
## Generar datasets de cada set
training_set <- splited_data$train %>% as_tibble()
validation_set <- splited_data$test %>% as_tibble()
## Plot por clase en cada set
ggplot(training_set, aes(Survived)) +
geom_bar(aes(fill = Survived))
ggplot(validation_set, aes(Survived)) +
geom_bar(aes(fill = Survived))
Modelos de regresión lógistica
Inicialmente se construye un modelo, para estimar la probabilida de sobrevivir, en función de la clase, el sexo y la edad.
## Modelo glm sobre titanic dada la formuka
model <- glm(Survived ~ Pclass+ Sex + Age ,family=binomial(link='logit'),data=training_set)
## Summary del modelo
paged_table(tidy(model))Todas las variables usadas en este modelo resultan significativas. La interpretación de los coeficientes está ligada a la manera en que las variables influyen en la probabilidad estimada. Por ejemplo, los hombres tienen menos probabilidad de sobrevivir, algo acorde a lo que se observó anteriormente. En tercera clase es donde mas se dismunuye la probabilidad de sobrevivir; la segunda clase también disminuye la probabilidad de vivir. La edad es la variable que menos afecta la probabilidad.
A continuación se simula la probabilidad de supervivencia de dos sujetos fictisios y con ninguna relación a la realidad. * Individuo 1: Rose que es una mujer de 17 años que viaja en primera clase * Individuo 2: Jack que es un hombre de 20 años viajando en tercera clase
## Generar df para testear
sujetos <- data.frame(
Pclass = as.factor(c(1, 3)),
Sex = c('female', 'male'),
Age = c(17.0, 20.0)
)
## Aplicar el modelo al dataset creado
pred <- predict(model, newdata = sujetos, type = 'response')
pred## 1 2
## 0.94061295 0.09702124En este caso Rose tiene cerca del 95.1% de probabilidad de sobrevivir, mientras que Jack tiene un 9.4%, quizá esto condicionó el destino.
Generación de modelos
Ahora se incluirán tres modelos, con el fin de establecer una comparación para el poder predictivo de cada uno y de esta manera poder encontrar aquel que estime de mejor manera el problema.
Para ello se definen las formulas de los diferentes modelos, estas formulas serán aplicadas mas adelante para ajustar los modelos.
logit_formulas <- formulas(.response = ~Survived, # Variable respuesta de los modelos
base = ~ Pclass + Sex + Age, ## Formula de los modelos creados
modelo1 = ~ Pclass + Sex + Age + Parch,
modelo2 = ~ Pclass + Sex + Age + Embarked,
modelo3 = ~ Pclass + Sex + Age + Embarked + Fare
)Una vez definidas las formulas de los diferentes modelos, se procede a aplicar cada una al conjunto de entrenamiento. Cada formula es mapeada a un modelo lineal generalizado, estos modelos serán de la familia binomial.
models <- tibble(logit_formulas) %>% ## Df de las formulas de los modelos
mutate(
models = names(logit_formulas), ## Nombres de los modelos creados
expression = paste(logit_formulas), ## Formulas para cada modelo
mod = map(logit_formulas, ~glm(.,family = 'binomial', data = training_set)) ##Aplicar cada modelo al training_set
) Una vez mapeado esto, el siguiente paso es comparar el ajuste de cada modelo
coeficientes <- models %>%
filter(models %in% c('base','modelo1','modelo2','modelo3')) %>% ## Modelos a evaluar
mutate(tidy = map(mod,tidy)) %>% ## Coeficientes para cada modelo
unnest(tidy, .drop = TRUE) %>% ## separar la información
mutate(estimate=round(estimate,5), ##Calcular significatividad
p.value=round(p.value,4)) %>%
select(c(models,expression,term, estimate,p.value, std.error, statistic)) ## Variables a mostrar
paged_table(coeficientes)En el segundo modelo la variable de la puerta de embarque no resultó significativa; en el modelo tres la variable de la puerta de emabarque y la tarifa del pasajero tampoco resultaron significativas. Ambas variables están directamente relacionadas con la clase en la que viajaban los pasajeros, así que es muy probable que esté influenciado por esto.
Para obtener mas información sobre la bondad del ajuste se usa la función Glance, esta permite obtener resúmenes de los modelos, como: medidas de bondad de ajuste, p-values para pruebas de hipótesis sobre residuos o información de convergencia del modelo.
eval_models <- models %>%
mutate(glance = map(mod,glance)) %>% ## Aplico glance sobre los modelos
unnest(glance, .drop = TRUE) %>% ## Separar lo obtenido
mutate(perc_explained_dev = 1-deviance/null.deviance) %>% ## Calcular de la deviance explicada
select(c(expression,
logLik,deviance,df.residual,perc_explained_dev))%>% ## Seleccionar campos
arrange(deviance)
kable(eval_models)%>%
kable_styling(bootstrap_options = c("striped", "hover"))| expression | logLik | deviance | df.residual | perc_explained_dev |
|---|---|---|---|---|
| Survived ~ Pclass + Sex + Age + Embarked + Fare | -280.0811 | 560.1623 | 615 | 0.3164049 |
| Survived ~ Pclass + Sex + Age + Embarked | -280.1754 | 560.3509 | 616 | 0.3161747 |
| Survived ~ Pclass + Sex + Age + Parch | -281.0025 | 562.0049 | 617 | 0.3141562 |
| Survived ~ Pclass + Sex + Age | -281.9690 | 563.9380 | 618 | 0.3117971 |
Los modelos minimizan la tasa de deviance de manera similar, sin embargo, se opta por elegir el modelo $Survived <- Pclass + Sex + Age $ ya que presenta un deviance similar a los otros modelos, pero en este modelo todas las variables resultan significativas.
Evaluación del modelo
Sobre este modelo seleccionado se realizan las predicciones para el conjunto de datos, y de esta manera se obtiene el area bajo la curva y un diagrama de violín, en el conjunto de entrenamiento.
# Añadir las predicciones sobre los datos de entrenamiento
models <- models %>%
mutate(pred= map(mod,augment, type.predict = "response"))
# Modelo completo de entrenamiento
prediction_full <- models %>%
filter(models=="base") %>% ## Modelo seleccionado
unnest(pred, .drop=TRUE) ##Separar los datos
## Calcular area bajo la curva en entrenamiento
roc_full <- roc(response=prediction_full$Survived, predictor=prediction_full$.fitted)
## Plot de curva auc
ggroc(list(full=roc_full), size=1) + geom_abline(slope = 1, intercept = 1, linetype='dashed') + theme_bw() + labs(title='Curvas ROC', color='Modelo')
## Imprimir area bajo la curva en entrenamiento
print(paste('AUC: training_set', roc_full$auc))## [1] "AUC: training_set 0.839331234899031"Lo primero que se puede observar es que el resultado está por encima de la recta del azar. En el eje Y se encuentra la tasa de verdaderos positivos, mientras que en X la tasa de falsos positivos, el modelo ideal se situaría en 1 en TPR (eje Y) y 0 en FPR (eje X). Se busca un punto de equilibrio donde haya una buena tasa de TPR pero sin aumentar mucho la tasa de FPR. En este caso cerca de 0.75 en X y de 0.77 en Y parece ser un buen punto.
## Grafico de violin para la clase vs la distribuciónde probabilidades
ggplot(prediction_full,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')
En este caso no hay claro punto para establecer como umbral, en caso de querer tener mayor ‘precisión’ al momento de estimar las personas que no van a sobrevivir, un buen punto de corte sería 0.75, con ese umbral habrían muchos casos en que se estimaría que una presona no sobrevive pero realmente si lo hace. Si por el contrario se quiere mayor acierto al estimar los que van a sobrevivir un buen punto parece estar cerca de 0.25.
Elección punto de corte
Con el fin de evitar un sobreajuste en el modelo se procede a obtener el punto de corte de manera óptima a través de le evaluación en el dataset de validación.
# Añadimos predicciones pero esta vez del dataset de validación
validation_model <- models %>%
filter(models == 'base') %>% ## modelo a aplicar
mutate(pred = map(mod, augment, newdata = validation_set, type.predict = "response"))
# Obtenemos las correspondientes a nuestro mejor modelo
validation_prediction <- validation_model %>%
filter(models == "base") %>% ## mejor modelo
unnest(pred, .drop = TRUE)
prediction_metrics <- function(cutoff, predictions=validation_prediction){
table <- predictions %>% ## clasificar cada punto si es mayor al punto de corte como positivo
mutate(predicted_class=if_else(.fitted>cutoff, 1, 0) %>% as.factor(),
Survived= factor(Survived))
confusionMatrix(table(table$predicted_class, table$Survived), positive = "1") %>%
tidy() %>% ## matriz de confusión que se obtendría con ese punto de corte
select(term, estimate) %>%
filter(term %in% c('accuracy', 'sensitivity', 'specificity', 'precision','recall')) %>%
mutate(cutoff=cutoff)
}
## Rango del punto de corte
cutoffs = seq( min(validation_prediction$.fitted),
max(validation_prediction$.fitted),
0.01)
## calcular las métricas sobre los rangos
logit_pred= map_dfr(cutoffs, prediction_metrics)%>% mutate(term=as.factor(term))
## gtafico de cada una de las metricas en cada punto de corte añadiendo linea del cruce
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="") +
geom_vline(aes(xintercept = 0.345), color = 'black', linetype = 'dashed')
Como se había mencionado anteriormente se busca un punto que maximice la tasa de verdaderos positivos y se disminuya la tasa de falsos positivos. El punto de corte entre estos está en 0.35, así que este es el punto seleccionado.
## clasificar laos puntos por encima del punto de corte calculado
table= validation_prediction %>%
mutate(predicted_class=if_else(.fitted>0.345, 1, 0) %>% as.factor(),
Survived= factor(Survived))
## matriz de confusion sobre esas estimaciones
confusionMatrix(table(table$Survived, table$predicted_class), positive = "1")## Confusion Matrix and Statistics
##
##
## 0 1
## 0 123 32
## 1 24 89
##
## Accuracy : 0.791
## 95% CI : (0.7374, 0.8381)
## No Information Rate : 0.5485
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.5756
##
## Mcnemar's Test P-Value : 0.3496
##
## Sensitivity : 0.7355
## Specificity : 0.8367
## Pos Pred Value : 0.7876
## Neg Pred Value : 0.7935
## Prevalence : 0.4515
## Detection Rate : 0.3321
## Detection Prevalence : 0.4216
## Balanced Accuracy : 0.7861
##
## 'Positive' Class : 1
## En este caso se obtiene un accuracy de 0.791 con un interalo de confianza del 95% que nos dice que el accuracy está entre 0.7374 y 0.8381. Se puede ver que el modelo distingue bien ambos casos, sin inclinar el modelo hacia una sola clase, es decir tiene un buen rendimiento para detecetar los que sobreviven y los que no, ambos cerca del 0.79.
Dataset testeo
Se lee el nuevo dataset, este no ha sido usado para entrenar el modelo, con el fin de evaluar de manera correcta el ajuste del modelo, se realizan las transformaciones necesarias sobre el mismo.
##leer nuevo dataset
titanic_test <- read_csv("titanic_complete_test.csv")
paged_table(titanic_test[1:4,]) ##Imprimir primeras filas## trasnformaciones necesarias
titanic_test <- titanic_test %>%
mutate( Survived = factor(Survived),
Pclass = factor(Pclass ),
Embarked = factor(Embarked )) %>%
select( -c( Name, Ticket, Cabin))Se procede a entrenear el modelo final, usando todo el conjunto de entranemiento, para clasificar los nuevos datos, usando como punto de corte, el encontrado anteriormente.
## punto de corte encontrado
sel_cutoff = 0.345
# Creamos el modelo con todo el conjunto de entrenamiento
full_model <- glm(logit_formulas$base, family = 'binomial', data = titanic)
# Agregamos la predicciones al dataset de testeo (nunca usado en el modelo)
table= augment(x=full_model, newdata=titanic_test, type.predict='response')
# Clasificamos utilizamos el punto de corte calculado
table=table %>% mutate(predicted_class=if_else(.fitted>0.345, 1, 0) %>% as.factor(),
Survived= factor(Survived))
# Creamos la matriz de confusión sobre la clasificación
confusionMatrix(table(table$Survived, table$predicted_class), positive = "1")## Confusion Matrix and Statistics
##
##
## 0 1
## 0 185 76
## 1 36 121
##
## Accuracy : 0.7321
## 95% CI : (0.6869, 0.7739)
## No Information Rate : 0.5287
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.4563
##
## Mcnemar's Test P-Value : 0.0002286
##
## Sensitivity : 0.6142
## Specificity : 0.8371
## Pos Pred Value : 0.7707
## Neg Pred Value : 0.7088
## Prevalence : 0.4713
## Detection Rate : 0.2895
## Detection Prevalence : 0.3756
## Balanced Accuracy : 0.7257
##
## 'Positive' Class : 1
## En este caso el accuracy es de 0.7321 con un intervalo de confianza del 95% está entre 0.6869 y 0.7739, como se esperaba es menor al encontrado en el conjunto de validación, se observa que la precisión es mayor en el caso de los positivos que en los negativos, es decir el modelo está asignando algunos casos que no son positivos como positivos en mayor medida.
Sin embargo, es un buen rendimiento, con un modelo que presenta todas sus variables significativas, y que resulta fácil de interpretar debido a las pocas variables que tiene. En este modelo la probabilidad de sobrvivir aumenta a medida que crece la clase en que viaja el pasajero, si es hombre disminuye la probabilidad de vivir, la edad disminuye la probabilidad pero en menor medida.