En este trabajo ajustaremos una serie de modelos de regresión logística utilizando un datset sobre pasajeros del Titanic, barco que se hundió en 1912.
El objetivo, en todos los casos, consistirá en predecir la cantidad sobrevivientes al naufragio. Los modelos que ajustemos nos proporcionarán la probabilidad de supervivencia. Elegiremos uno y, luego de establecer un punto de corte, podremos predecir la clase a la que pertenece cada pasajero.
Asimismo, se analizarán las características de los datasets utilizados, como las de los modelos generados.
# se cargan las librerías a utilizar
library("tidyverse")
library("data.table")
library("ggplot2")
library("GGally")
library("modelr")
library("gridExtra")
library("caret")
library("broom")
library("pROC")
# se elige una semilla
set.seed(1912)
2. Preparación de los datos
En primer lugar, se carga el dataset. Este cuenta con 891 registros y 12 variables:
PassengerId: Id único por pasajero
Survived: sobreviviente (1) o fallecido (0)
Pclass: clase del viaje (1ra, 2da o 3ra)
Name: nombre del pasajero
Sex: sexo del pasajero
Age: edad del pasajero
SibSp: cantidad de hermanos y/o cónyuge a bordo
Parch: cantidad de padres/hijos a bordo
Ticket: número de ticket
Fare: tarifa pagada por el ticket
Cabin: número de cabina
Embarked: lugar de abordo (C: Cherbourg, S: Southampton, Q: Queenstown)
# se carga el dataset
dataset <- fread("titanic_complete_train.csv")
glimpse(dataset)
De éstas, nos seleccionamos con PassengerId, Sex, Age, SibSp, Parch, Fare, Pclass, Survived y Embarked, que son las que utilizaremos para el ajuste de los distintos modelos, y convertimos las últimas tres a factores.
El nuevo dataset consta de 891 registros y 9 variables.
# se seleccinan las variables deseadas y se convierten algunas a factores
dataset <- dataset %>%
select(PassengerId, Survived, Pclass, Sex, Age, SibSp, Parch, Fare, Embarked) %>%
mutate(Survived = as_factor(Survived),
Pclass = as_factor(Pclass),
Embarked = as_factor(Embarked))
dataset
Con el objetivo de visualizar la relación entre distintas variables, se realiza un gráfico de ggpairs de los atributos Survived, Pclass, Sex, Age y Fare.
En primer lugar, se observa que la variable Sobreviviente contiene clases desbalanceadas: la cantidad de fallecidos supera a la de sobrevivientes.
Además, al relacionar esta variable con la clase en la cual viajaban los pasajeros, vemos tres agrupaciones. En la primera, posiblemente correspondiente a la primera clase, hay mayor cantidad de sobrevivientes que de fallecidos; en la segunda clase, ambas cantidades parecen ser similares, y en la tercera la relación se presenta a la inversa: hay más fallecidos que sobrevivientes.
Asimismo, al vincular la supervivencia con el sexo de los pasajeros se evidencian dos comportamientos diferenciados: en las mujeres hay más sobrevivientes que fallecidos, y en los hombres ocurre lo contrario.
Si también relacionamos esta información con la que aporta la variable de la edad, vemos que la media de las personas a bordo del Titanic se encontraba entre los 20 y 30 años y que también aquí se encuentra la media de los fallecidos, en su mayoría, varones de la tercera clase.
Apoyamos estas primeras aproximaciones a los datos con tablas que nos permiten confirmar nuestras hipotesis.
# porcentaje de sobrevivientes
dist.dataset <- round((342 / nrow(dataset)) * 100,0)
print(paste("Porcentaje de sobrevivientes en el dataset:", dist.dataset,"%"))
[1] "Porcentaje de sobrevivientes en el dataset: 38 %"
Con el objetivo de poder entrenar un modelo y validarlo previamente a realizar nuestras predicciones, dividimos el dataset en un conjunto de entrenamiento (70% de los datos) y otro de validación (30%). Para ello usamos la función resample_partition de modelr.
Constatamos también que ambos conjuntos presentan aproximadamente la misma distribución en la variable Sobreviviente que la que mostró el dataset completo, dado que esta será nuestra clase target y queremos asegurarnos de entrenar y validar con conjuntos de datos semejantes.
# se particiona el dataset en entrenamiento y validación
train_valid <- dataset %>% resample_partition(c(train=0.7,valid=0.3))
trainset <- train_valid$train %>% as_tibble()
validationset <- train_valid$valid %>% as_tibble()
# se grafican las distribuciones de clases de ambos conjuntos
trainplot <- ggplot(trainset,
aes(x=Survived, fill=Survived)) +
geom_bar() +
geom_text(stat='count', aes(label=..count..), position = position_fill(vjust = 50)) +
labs(title = "Conjunto de entrenamiento",
y = "Frecuencia absoluta") +
scale_x_discrete("Sobrevivientes", breaks = c("0","1"),
labels=c("No sobrevivientes","Sobrevivientes")) +
theme(legend.position = "none")
validplot <- ggplot(validationset,
aes(x=Survived, fill=Survived)) +
geom_bar() +
geom_text(stat='count', aes(label=..count..), position = position_fill(vjust = 22)) +
labs(title = "Conjunto de validación",
y = "Frecuencia absoluta") +
scale_x_discrete("Sobrevivientes", breaks = c("0","1"),
labels=c("No sobrevivientes","Sobrevivientes")) +
theme(legend.position = "none")
grid.arrange(
trainplot,
validplot,
nrow = 1,
top = "Sobrevivientes - Distribución de clases"
)
# porcentaje de sobrevivientes en set de entrenamiento
dist.train <- round((240 / nrow(trainset)) * 100,0)
print(paste("Porcentaje de sobrevivientes en el conjunto de entrenamiento:", dist.train,"%"))
[1] "Porcentaje de sobrevivientes en el conjunto de entrenamiento: 39 %"
# porcentaje de sobrevivientes en set de validación
dist.val <- round((102 / nrow(validationset)) * 100,0)
print(paste("Proporción de sobrevivientes en el conjunto de validación:", dist.val,"%"))
[1] "Proporción de sobrevivientes en el conjunto de validación: 38 %"
3. Generación de modelos
Procedemos a ajustar 4 modelos de regresión logísitca utilizando el conjunto de entrenamiento:
modelo 0: predice la probabilidad de sobrevivir en función de la clase, el sexo y la edad
modelo 1: predice la probabilidad de sobrevivir en función de la tarifa pagada por el pasaje y el sexo
modelo 2: predice la probabilidad de sobrevivir en función de la clase, el sexo, la edad y la tarifa pagada por el pasaje
modelo 3: predice la probabilidad de sobrevivir en función del sexo y la cantidad de padres y/o hijos abordo del barco
# se definen las fórmulas para los 4 modelos
logit_formulae <- formulas(.response = ~Survived,
model0 = ~ Pclass + Sex + Age,
model1 = ~ Fare + Sex,
model2 = ~ Pclass + Sex + Age + Fare,
model3 = ~ Sex + Parch
)
# se genera un data frame con la información y los modelos ajustados
models <- data_frame(logit_formulae) %>%
mutate(models = names(logit_formulae),
expression = paste(logit_formulae),
mod = map(logit_formulae, ~glm(.,family = 'binomial', data = trainset)))
models
Nos detenemos a analizar los coeficientes del primer modelo ajustado, el modelo 0.
En una regresión logística, los coeficientes positivos indican que, al aumentar su variable, la probabilidad predicha aumenta también. Por el contrario, un coeficiente negativo indica que la probabilidad disminuye cuando la variable incrementa su valor.
En el caso particular del modelo 0, la probabilidad de sobrevivir disminuye cuando aumentan la variable Pclass toma los valores 2 o 3 y, en este último caso, la disminución es más pronunciada. Del mismo modo, si el pasajero resulta ser un varón, su probabilidad de sobrevivir se reduce. Y, conforme la persona sea de edad más avanzada, la probabilidad de sobrevivir también decaerá.
summary(models$mod$model0)
Call:
glm(formula = ., family = "binomial", data = trainset)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.6713 -0.7017 -0.4466 0.6619 2.4407
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.611951 0.461715 7.823 5.16e-15 ***
Pclass2 -1.188844 0.310697 -3.826 0.000130 ***
Pclass3 -2.512577 0.308365 -8.148 3.70e-16 ***
Sexmale -2.392464 0.218841 -10.932 < 2e-16 ***
Age -0.036296 0.009367 -3.875 0.000107 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 830.54 on 622 degrees of freedom
Residual deviance: 577.86 on 618 degrees of freedom
AIC: 587.86
Number of Fisher Scoring iterations: 4
Realizamos dos predicciones que nos permiten ver mejor la relación entre coeficientes y probabilidad predicha. Por un lado, predecimos la probabilidad de sobrevivir de Rose, una mujer de 17 años que viaja en primera clase. Por otro, evaluamos el caso de Jack, varón de 20 años que viaja en tercera clase.
rose <- data.frame(Pclass = as.factor(1),
Sex = "female",
Age = 17)
rose.prob <- predict(object = models$mod$model0, newdata = rose, type = "response")
print(paste("Probabilidad de que Rose sobreviva:", round(rose.prob,3), sep= " "))
[1] "Probabilidad de que Rose sobreviva: 0.952"
jack <- data.frame(Pclass = as.factor(3),
Sex = "male",
Age = 20)
jack.prob <- predict(object = models$mod$model0, newdata = jack, type = "response")
print(paste("Probabilidad de que Jack sobreviva:", round(jack.prob,3), sep = " "))
[1] "Probabilidad de que Jack sobreviva: 0.117"
Como era de esperar, la probabilidad de sobrevivir predicha para Jack es menor. En el caso de Rose, la probabilidad no se ve disminuida por los coeficientes de las vairbales Pclass ni Sex, sino solo por el de Age. Jack, en cambio, no solo se ve afectado por este último (y en mayor medida dado que su edad es mayor), sino también por los coeficientes asociados a Pclass3 (Pclass al adoptar el valor 3) y Sexmale (Sex cuando toma el valor male), que tienen una influencia negativa en la probabilidad predicha para él.
Modelo 1: tanto la variable Fare como Sex resultan significativas. Mientras que la variable Fare (la tarifa pagada por el pasaje), al aumentar, incrementa la proabilidad de sobrevivir, la variable Sex produce ele fecto contrario al tomar el valor male.
Modelo 2: aquí todas las variables, a excepción de Fare, resultan significativas. El coeficiente asociados a Fare presenta un valor mayor a 0,05, lo que indica que esta variable no aporta información al modelo. En el caso de Pclass, Sex y Age se observa una influencia en la probabilidad similar a la encontrada en el modelo 0.
Modelo 3: aquí la variable Sex aporta infromación pero la variable Parch no. En el caso de Sex, tiene la misma influencia sobre la probabilidad predicha que la que tenía en el modelo 0: al adoptar el valor male, reduce la probabilidad. Parch, por su parte, influye negativamente al aumentar su valor.
A continuación, se ordenan los modelos por deviance explicada. Podemos ver que el modelo que minimiza la deviance es el modelo 2, lo que hace que su porcentaje de deviance explicada sea mayor al resto: aproximadamente 30%.
En este apartado procedemos a evaluar el mejor modelo ajustado en el apartado anterior, esto es, el modelo 2.
En primer lugar, realizamos una curva ROC. En ella graficamos el Ratio de Verdaderos Negativos o specificity en el eje de abscisas y el Ratio de Verdaderos Positivos o sensitivity en el eje de ordenadas, cuyas fórmulas son:
\(sensitivity = \frac{TP}{TP+FN}\)
\(specificity = \frac{TN}{TN+FP}\)
Donde:
TP = true positive o verdaderos positivos (casos que el modelo predijo como positivos y realmente lo eran)
FP = false positive o falsos positivos (casos que el modelo predijo como positivos y en realidad eran negativos)
TN = true negative o verdadernos negativos (casos que el modelo predijo como negativos y realmente lo eran)
FN = false negative o falsos negativos (casos que el modelo predijo como negativos y en realidad eran positivos)
La sensitivity, a la cual también se la suele denominar recall, indica la cuántos casos realmente positivos se lograron predecir correctamente, es decir, cuánta cobertura tiene el modelo. Para ello, se considera el total de casos que el modelo predijo como positivos y que realmente lo son (TP) sobre el total de positivos relaes (TP+FN), hayan sido predichos como tales por el modelo (TP) o hayan sido catalogados errónemanete como negativos (FN). En este trabajo, la clase positiva son los sobrevivientes.
La specificity, por otro lado, indica lo mismo pero con respecto a los negativos: cuántos casos que en realidad eran negativos fueron predichos de esta forma por nuestro modelo. Aquí la cuenta que se realiza toma los negativos reales que fueron predichos de tal forma (TN) y los divide por el total real de negativos (TN+FP). En este trabajo, los casos negtivos corresponden a los “no sobrevivientes”.
Además, podemos calcular el área bajo la curva graficada. El modelo 2 cuenta con un área del 83%, por lo que sus predicciones superan las del azar, representado por la línea punteada y cuya área es del 50%.
Graficamos también un violin plot que nos permita ver la distribución de las clases “sobreviviente” y “no sobreviviente” en relación a la probabilidad predicha por el modelo.
En él se puede observar que, a medida que la probabilidad asignada se incrementa, la clase “no sobreviviente” disminuye su densidad y la clase “sobreviviente” la aumenta. Esto nos ayuda a pensar un posible punto de corte a partir del cual podemos establecer que una persona sobrevive, intentando minimizar además el error de predicción (i.e. catalogar a un “no sobrevivientes” como “sobreviviente” o viceversa).
Dado que el quiebre entre disminución de la densidad de “no sobrevivientes” y aumento de la densidad de “sobrevivientes” parece estar entre 0.30 y 0.45, podría pensarse un punto de corte que esté en ese intervalo.
En este apartado procederemos a probar nuestro modelo con los datos que hemos reservado para la validación y, a partir de esos resultados, elegiremos un punto de corte que consideremos apropiado para predecir si alguien sobrevide.
# se testean el modelo con los datos de VALIDACIÓN reservados
valid <- models %>%
filter(grepl('2',models)) %>%
mutate(val= map(mod,augment, newdata=validationset, type.predict = "response"))
valid.results <- valid %>%
unnest(val, .drop=TRUE)
valid.results
Para ayudarnos, realizamos un gráfico donde se visualice accuracy, precision, recall y specificity.
Ya hemos hablado de los últimos dos, detengámonos un momento en los primeros:
\(accuracy = \frac{TP+TN}{TP+FP+TN+FN}\)
\(precision = \frac{TP}{TP+FP}\)
El accuracy permite medir la performance general del modelo. Para ello, se debe calcular cuántos casos fueron predichos correctamente, sean estos positivos (aquí, “sobrevivientes”) o negativos (“no sobrevivientes”), y relacionarlos con el total de casos (TP+FP+TN+FN), de modo que la cuenta resultante refleja proporción de los registros fue predicha de manera correcta respecto del total.
La precision, a su vez, indica cuán exacto es el modelo, dado que podría ocurrir que éste capturase correctamente los casos positivos pero a costa de predecir también como positivos los registros negativos. Para medir esto es que se usa la precision, que calcula la cantidad de verdaderos positivos que predijo el modelo sobre la cantidad total de positivos predichos, sean estos realmente positivos o no.
# función que calcula las métricas de las predicciones de acuerdo a los puntos de corte establecidos
prediction_metrics <- function(cutoff, predictions=valid.results){
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 definen puntos de corte
cutoffs = seq(0.05,0.95,0.01)
logit_pred= map_dfr(cutoffs, prediction_metrics) %>% mutate(term=as.factor(term))
# se seleccionan los puntos de corte con mayor accuracy
cutoff.max_ac <- logit_pred %>%
filter(term=='accuracy') %>%
filter(estimate == max(estimate)) %>%
select(cutoff)
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= 'Model2 - Sobre conjunto de validación', color="") +
geom_vline(xintercept = min(cutoff.max_ac), linetype="dotted", color = "black", size=1) +
geom_vline(xintercept = max(cutoff.max_ac), linetype="dotted", color = "black", size=1)
En el gráfico se observa que el mejor accuracy se obtiene entre los puntos de corte 0.45 y 0.46. A partir de allí, la línea que marca el accuracy comienza a decaer.
También es cerca de ese punto donde se intersectan las líneas de precision y recall (o sensitivity). Esta intersección da cuenta de la compensación o tensión en la que se encuentran la cobertura y la exactitud de un modelo: la capacidad de predecir correctamente mayor cantidad de casos positivos (i.e. de tener un mayor recall) está acompañada de una menor precisión, por la cual el modelo predice como positivos registros que en realidad no lo son. Y a medida que esta precisión aumenta, su precisión disminuye, es decir, predice correctamente menor cantidad de casos positivos.
De manera similar, la sensitivity se encuentra en tensión con la specificity: en tanto la primera disminuye, la segunda aumenta como consecuencia de que el modelo, al incrementar los falsos positivos y, en consecuencia, predecir menor cantidad de casos positivos correrctamente, reduce también las predicciones negativas que verdaderamente son de esta clase.
Un punto de corte adecuado para predecir si una persona es o no sobreviviente debería lograr cierta armonía entre estas tensiones. En este trabajo elegimos como punto de corte el valor 0.46 por ser aquel que logra el mayor accuracy con el conjunto de validación y encontrarse más cercano a los puntos donde se intersectan recall, precision y specificity.
A continuación, realizamos la matriz de confusión obtenida con tal punto de corte.
# se establece el punto de corte
sel_cutoff = max(cutoff.max_ac)
table <- valid.results %>%
mutate(predicted_class=if_else(.fitted>sel_cutoff, 1, 0) %>% as.factor(),
Survived = factor(Survived))
# se crea la matriz de confusión
confusionMatrix(table(table$Survived, table$predicted_class), positive = "1")
Confusion Matrix and Statistics
0 1
0 141 25
1 19 83
Accuracy : 0.8358
95% CI : (0.7859, 0.8781)
No Information Rate : 0.597
P-Value [Acc > NIR] : <2e-16
Kappa : 0.6557
Mcnemar's Test P-Value : 0.451
Sensitivity : 0.7685
Specificity : 0.8812
Pos Pred Value : 0.8137
Neg Pred Value : 0.8494
Prevalence : 0.4030
Detection Rate : 0.3097
Detection Prevalence : 0.3806
Balanced Accuracy : 0.8249
'Positive' Class : 1
Esta matriz nos permite visualizar la cantidad esperada de valores positivos y negativos y las predicciones realizadas por el modelo para dichas clases del siguiente modo:
Negativos esperados
Positivos esperados
Negativos predichos
Verdaderos negativos
Falsos negativos
Positivos predichos
Falsos positivos
Verdaderos positivos
Así, podemos ver que nuestro modelo predice adecuadamente 141 casos negativos y 83 positivos, y hay 25 casos que cataloga como negativos cuando en realidad son positivos y 19 con los que sucede lo contrario.
6. Testeo
En este apartado procedemos a testear el modelo seleccionado con un nuevo dataset no utilizado para entrenamiento ni testeo.
Cargamos el dataset y realizamos el mismo preprocesamiento que hicimos con el conjunto de entrenamiento: convertimos en factor las variables Survived, Pclass y Embarked.
El conjunto de testeo utilizado cuenta con 418 registros y presenta una distribución de clases similar a aquella con la cual entrenamos el modelo.
# se carga el conjunto de testeo
testset <- fread("titanic_complete_test.csv")
# se convierten las variables a factores
testset <- testset %>%
select(PassengerId, Survived, Pclass, Sex, Age, SibSp, Parch, Fare, Embarked) %>%
mutate(Survived = as_factor(Survived),
Pclass = as_factor(Pclass),
Embarked = as_factor(Embarked))
glimpse(testset)
# se calcula el porcentaje de sobrevivientes en el conjunto de testeo
surv.testset <- nrow(testset[testset$Survived==1,])
dist.testset <- round((surv.testset/nrow(testset))*100,0)
print(paste("Porcentaje de sobrevivientes en el conjunto de testeo:", dist.testset,"%"))
[1] "Porcentaje de sobrevivientes en el conjunto de testeo: 38 %"
Dado que ya hemos elegido el punto de corte y realizamos la validación, volvemos a entrenar nuestro modelo con todo el conjunto de entrenamiento, de modo que este cuente con una mayor cantidad de registros para realizar las predicciones.
# se entrena el modelo 2 con TODOS los datos para ENTRENAMIENTO
model2test <- glm(logit_formulae$model2, family = 'binomial', data = dataset)
# se calculan las predicciones para el conjunto de TESTEO
test.table = augment(x=model2test, newdata=testset, type.predict='response')
# se realiza la clasificación en clases según el punto de corte elegido
test.table = test.table %>%
mutate(predicted_class=if_else(.fitted>sel_cutoff, 1, 0) %>% as.factor(),
Survived= factor(Survived))
test.table
Finalmente, realizamos la matriz de confusión y obtenemos las métricas para las nuevas predicciones.
Es posible ver que el accuracy disminuyó, lo que era esperable dado que las predicciones se realizaron sobre un conjunto de datos no visto por el modelo previamente. De todos modos, un 75% de los casos de predicen correctamente; de estos, 206 corresponden a registros negativos y representan el 81% del total de casos negativos y 110 corresponden a registros positivos y, de acuerdo a la sentivity, constituyen el 66% del total de positivos.
Se ve de este modo que nuestor modelo predice mejor los casos negativos antes que los poritivos.
# se crea la matriz de confusión
confusionMatrix(table(test.table$Survived, test.table$predicted_class), positive = "1")
Confusion Matrix and Statistics
0 1
0 206 55
1 47 110
Accuracy : 0.756
95% CI : (0.7119, 0.7964)
No Information Rate : 0.6053
P-Value [Acc > NIR] : 5.475e-11
Kappa : 0.485
Mcnemar's Test P-Value : 0.4882
Sensitivity : 0.6667
Specificity : 0.8142
Pos Pred Value : 0.7006
Neg Pred Value : 0.7893
Prevalence : 0.3947
Detection Rate : 0.2632
Detection Prevalence : 0.3756
Balanced Accuracy : 0.7404
'Positive' Class : 1
