eea2020_tp3_cisco_nicolas
Nicolas Cisco
set.seed(3717)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(GGally)## Loading required package: ggplot2## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2library(ggplot2)
library(rsample)
library(tidyr)
library(purrr)
library(broom)
library(pROC)## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, varlibrary(caret)## Loading required package: lattice##
## Attaching package: 'caret'## The following object is masked from 'package:purrr':
##
## liftEntrenamiento
Preparación de datos
dataset.train <- read.csv('./titanic_complete_train.csv')
glimpse(dataset.train)## Rows: 891
## Columns: 12
## $ PassengerId <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…
## $ Survived <int> 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, …
## $ Pclass <int> 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 (F…
## $ Sex <chr> "male", "female", "female", "female", "male", "male", "ma…
## $ Age <dbl> 22.00000, 38.00000, 26.00000, 35.00000, 35.00000, 26.5075…
## $ SibSp <int> 1, 1, 0, 1, 0, 0, 0, 3, 0, 1, 1, 0, 0, 1, 0, 0, 4, 0, 1, …
## $ Parch <int> 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", "3…
## $ Fare <dbl> 7.2500, 71.2833, 7.9250, 53.1000, 8.0500, 8.4583, 51.8625…
## $ 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", "S…El dataset tiene doce columnas:
PassengerId: entero que identifica al pasajero (único por cada fila)Survived: booleano que representa si el pasajero sobrevivió:0: No sobrevivió1: Sobrevivió
Pclass: categórico que representa la clase del ticket:1: Primera clase2: Segunda clase3: Tercera clase
Name: cadena de texto que representa el Nombre de la personaSex: categórico que representa el sexo de la persona:male: Masculinofemale: Femenino
Age: entero que representa la edad de la personaSibSp: entero que representa la cantidad de hermanos / esposas abordo del TitanicParch: entero que representa la cantidad de padres / hijos abordo del TitanicTicket: cadena de texto que identifica al ticketFare: double que representa la tarifa del ticketCabin: cadena de texto que identifica la cabinaEmbarked: categórico que represeta el puerto en donde se embarcó:C: CherbourgQ: CherbourgS: Southampton
columns <- c(
'PassengerId',
'Survived',
'Pclass',
'Sex',
'Age',
'SibSp',
'Parch',
'Fare',
'Embarked'
)
dataset.train.clean <- dataset.train %>%
select(one_of(columns)) %>%
mutate(
Survived = as.factor(Survived),
Pclass = as.factor(Pclass),
Embarked = as.factor(Embarked)
)
ggpairs(
data = dataset.train.clean,
columns = c(
'Survived',
'Pclass',
'Sex',
'Age',
'Fare'
),
upper = list(continuous = ggally_cor_v1_5)
)## `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`.
Se obvserva que en la población general la mayoría pasajeros no sobrevivieron. Esta tendencia no esta presente si se divide a la población por Pclasso por Sex. Se puede notar que los individuos con Pclass=1 o Sex=female tuvieron una mayor tendencia a sobrevivir. Hay que notar, que aunque las distrubciones de Age agrupadas por Survived parecen ser normales con un pico entre los 20 y 30 años, el grupo que sobrevivió tiene un pico más marcado, que es muy similar al gráfico de la distribución de Age agrupado por Sex para el conjunto de Sex=female. Esto último puede estar relacionado con la mayor tendencia de sobrevivir de que los individuos de Sex=female antes mencionada.
Si se miran a los inviduos agrupados por Pclass, se puede notar diferencias entre las medias de la edad de cada clase, siendo Pclass=1 la clase con un Age más alto, seguido por Pclass=2 y por último Pclass=3 con Age más bajo.
Al observar a los indivudos agrupados por Sex, se nota que hay una mayor proporcion de mujeres en el conjunto de datos. Y aunque las medias de las edades parecieran ser bastante similares, la de los hombres sería un poco menor a la de las mujeres.
En resumen, a simple vista pareciera que las variables más importantes para determinar si un individuo sobrevivió podrían ser Sex y Pclass. Siendo los invididuos con Pclass=1 y Sex=female los que más sobrevivieron.
ggplot(
dataset.train.clean,
aes(x = Survived)
) +
geom_bar()
dataset.train.clean %>%
group_by(Survived) %>%
summarise(n = n()) %>%
mutate(frec = n / sum(n))## `summarise()` ungrouping output (override with `.groups` argument)| Survived | n | frec |
|---|---|---|
| 0 | 549 | 0.6161616 |
| 1 | 342 | 0.3838384 |
Dividir al dataset en conjunto de entrenamiento (70% de los datos) y validación (30% de los datos).
dataset.train.split <- initial_split(dataset.train.clean, prop = 0.7)
dataset.train.train <- training(dataset.train.split)
dataset.train.validation <- testing(dataset.train.split)
ggplot(
dataset.train.train,
aes(x = Survived)
) +
geom_bar()
dataset.train.train %>%
group_by(Survived) %>%
summarise(n = n()) %>%
mutate(frec = n / sum(n))## `summarise()` ungrouping output (override with `.groups` argument)| Survived | n | frec |
|---|---|---|
| 0 | 385 | 0.6169872 |
| 1 | 239 | 0.3830128 |
ggplot(
dataset.train.validation,
aes(x = Survived)
) +
geom_bar()
dataset.train.validation %>%
group_by(Survived) %>%
summarise(n = n()) %>%
mutate(frec = n / sum(n))## `summarise()` ungrouping output (override with `.groups` argument)| Survived | n | frec |
|---|---|---|
| 0 | 164 | 0.6142322 |
| 1 | 103 | 0.3857678 |
Se observa que la distribución de la clase Survived es aproximadamente similar en el conjunto original y en los generados ( train y validation ). Siendo aproximadamente un 61% para el grupo que no sobrevivió y 39% para el que sí.
Predicciones
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.
model.glm <- glm(Survived ~ Pclass + Sex + Age, data = dataset.train.train, family = binomial)
tidy(model.glm)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 3.7204552 | 0.4520272 | 8.230601 | 0.0000000 |
| Pclass2 | -1.2013231 | 0.3146754 | -3.817658 | 0.0001347 |
| Pclass3 | -2.5174266 | 0.3101330 | -8.117249 | 0.0000000 |
| Sexmale | -2.6480600 | 0.2278113 | -11.623917 | 0.0000000 |
| Age | -0.0338531 | 0.0091405 | -3.703621 | 0.0002125 |
Se observa que todos los coeficientes son estadísticamente significativos.
Hay que destacar en el modelo de regresión logística las variables no tienen una relación lineal con la probabilidad. En este modelo un coeficiente positivo indica que frente a aumentos de dicha variable la probabilidad aumenta, mientras que un coeficiente negativo nos indica lo contrario.
Particularmente:
- Dadas el resto de las variables constantes, el coeficiente
Pclass2(-1.20132307) indica que cuando la clase del invidivuo es dos disminuye la probabilidad de sobrevivir respecto a la categoría basal. - Dadas el resto de las variables constantes, el coeficiente
Pclass3(-2.51742659) indica que cuando la clase del invidivuo es tres disminuye la probabilidad de sobrevivir respecto a la categoría basal. - Dadas el resto de las variables constantes, el coeficiente
Sexmale(-2.64805998) indica que cuando el sexo del invidivuo es hombre disminuye la probabilidad de sobrevivir respecto a la categoría basal. - Dadas el resto de las variables constantes, El coeficiente
Age(-0.03385308) indica que cuando la edad del invidivuo aumenta la probabilidad de sobrevivir disminuye.
¿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.
people <- data.frame(
Name = c('Rose', 'Jack'),
Pclass = as.factor(c(1, 3)),
Sex = as.factor(c('female', 'male')),
Age = c(17, 20)
)
data.frame(
Nombre = people$Name,
`Probabilidad Supervivencia` = predict(model.glm, newdata = people, type = 'response')
)| Nombre | Probabilidad.Supervivencia |
|---|---|
| Rose | 0.9587094 |
| Jack | 0.1069680 |
La probabilidad de sobrevivir de Rose es mayor que la de Jack (0.96 > 0.11)
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
El primer modelo se generó con todas las covariables seleccionadas del dataset:
model.1 <- glm(Survived ~ Pclass + Sex + Age + SibSp + Parch + Fare + Embarked, data = dataset.train.train, family = binomial)
tidy(model.1)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 4.4656675 | 0.5902406 | 7.5658421 | 0.0000000 |
| Pclass2 | -1.1665168 | 0.3617334 | -3.2247970 | 0.0012606 |
| Pclass3 | -2.5488556 | 0.3775372 | -6.7512701 | 0.0000000 |
| Sexmale | -2.7550186 | 0.2419246 | -11.3879245 | 0.0000000 |
| Age | -0.0400542 | 0.0098323 | -4.0737592 | 0.0000463 |
| SibSp | -0.3401253 | 0.1433072 | -2.3734001 | 0.0176252 |
| Parch | -0.0242890 | 0.1696520 | -0.1431698 | 0.8861561 |
| Fare | -0.0006014 | 0.0029689 | -0.2025682 | 0.8394726 |
| EmbarkedQ | -0.0022472 | 0.4745475 | -0.0047355 | 0.9962216 |
| EmbarkedS | -0.4095513 | 0.2870198 | -1.4269097 | 0.1536059 |
Se tiene que resaltar que los coeficientes Parch, Fare, EmbarkedQ y EmbarkedS no resultan estadísticamente significativos. Todos los coeficientes resultan negativos, por lo tanto, si variamos las covariables de a una, significa que al estar presente respecto a la categoría basal (en el caso de las covariables categóricas: Pclass2, Pclass3, Sexmale, EmbarkedQ y EmbarkedS) o de aumentar (en el caso de las covariables númericas: Age, SibSp y Fare) disminuye la probabilidad de sobrevivir.
El segundo modelo se realizó con todas las variables que resultaron estadisticamente significativas del anterior modelo.
model.2 <- glm(Survived ~ Pclass + Sex + Age + SibSp, data = dataset.train.train, family = binomial)
tidy(model.2)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 4.2331189 | 0.4991022 | 8.481467 | 0.0000000 |
| Pclass2 | -1.2529205 | 0.3187579 | -3.930634 | 0.0000847 |
| Pclass3 | -2.5524949 | 0.3139647 | -8.129878 | 0.0000000 |
| Sexmale | -2.7727293 | 0.2377029 | -11.664683 | 0.0000000 |
| Age | -0.0412525 | 0.0097198 | -4.244152 | 0.0000219 |
| SibSp | -0.3782431 | 0.1324183 | -2.856427 | 0.0042844 |
Todos los coeficientes resultan estadísticamente significativos. Al ser todos negativos, si mantenemos el resto de las covariables constantes, disminuyen la probabilidad si la covariable esta presente respecto a la categoría basal (en el caso de los categóricos) o aumenta (en el caso de las numéricas).
Por último se realiza un modelo utilizando solo las variables que se supusieron como más importantes en la sección de procesamiento de datos: Pclass y Sex:
model.3 <- glm(Survived ~ Pclass + Sex, data = dataset.train.train, family = binomial)
tidy(model.3)| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 2.4604148 | 0.2700217 | 9.111917 | 0.0000000 |
| Pclass2 | -0.9121735 | 0.2990050 | -3.050696 | 0.0022831 |
| Pclass3 | -1.9892349 | 0.2628211 | -7.568778 | 0.0000000 |
| Sexmale | -2.7430724 | 0.2253309 | -12.173530 | 0.0000000 |
Todas los coeficientes resultan estadísticamente significativos. Al ser todos negativos, si mantenemos el resto de las covariables constantes, disminuye la probabilidad de sobrevivir respecto a la categoría basal del individuo si la covariable esta presente.
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.
models <- data.frame(
name = c(
'Punto 2.a',
'1',
'2',
'3'
),
model = I(list(
model.glm,
model.1,
model.2,
model.3
))
) %>%
mutate(glance = map(model, glance)) %>%
unnest(glance) %>%
mutate(perc_explained_dev = 1-deviance/null.deviance) %>%
select(-c(model, df.null, AIC, BIC)) %>%
arrange(deviance)
models| name | null.deviance | logLik | deviance | df.residual | nobs | perc_explained_dev |
|---|---|---|---|---|---|---|
| 1 | 830.5687 | -265.2517 | 530.5035 | 614 | 624 | 0.3612768 |
| 2 | 830.5687 | -266.5935 | 533.1870 | 618 | 624 | 0.3580458 |
| Punto 2.a | 830.5687 | -271.4448 | 542.8895 | 619 | 624 | 0.3463641 |
| 3 | 830.5687 | -278.7085 | 557.4171 | 620 | 624 | 0.3288730 |
Se selecciona el modelo 1, el cual utiliza todas las covariables, ya que es el mejor en términos de la deviance explicada. Aunque hay que destacar que el modelo 2, el que remueve las variables que no parecen ser estadísticamente significativas, es muy similar en términos de deviance.
Evaluación del modelo
Realizar el gráfico de curva ROC y obtener el AUC para el modelo elegido. Interpretar el gráfico.
model.1.roc <- roc(response = dataset.train.train$Survived, predictor = model.1$fitted.values)## Setting levels: control = 0, case = 1## Setting direction: controls < casesprint(model.1.roc)##
## Call:
## roc.default(response = dataset.train.train$Survived, predictor = model.1$fitted.values)
##
## Data: model.1$fitted.values in 385 controls (dataset.train.train$Survived 0) < 239 cases (dataset.train.train$Survived 1).
## Area under the curve: 0.8647ggroc(model.1.roc) +
geom_abline(slope = 1, intercept = 1, linetype='dashed') +
theme_bw() +
labs(title='Curva ROC')
El gráfico de curva ROC es una representación gráfica de la sensibilidad frente a la especificidad para un clasificador binario, es decir, del porcentaje de individuos clasificados correctamente positivos dentro de los individuos que son realmente positivos frente al porcentaje de individuos correctamente clasificados negativos dentro de los individuos clasificados negativamente, según se varía el umbral de clasificación (o cutoff). Un clasificador al azar supondría la recta x = y (línea punteada en el gráfico), ese clasificador tendría un area bajo la curva (auc, por su siglas en inglés) de 0.5.
En este caso, el area bajo la curva es de 0.8647 lo que nos dice que es un buen modelo.
Realizar un violin plot e interpretar.
ggplot(
augment(model.1, type.predict = 'response'),
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')
Los gráficos de violin muestran la distribución de probabilidades para cada una de las clases, en este caso, los pasajeros que sobrevivieron y los que no. En rojo esta la distribución de los individuos que no sobrevivieron, se puede apreciar como hay una mayor concentración cerca del cero, más especificamente menor a 0.25. Por otro lado, en verde se observa la distribución de probabilidades de los individuos que sobrevieron, en dicha distribución se nota que hay una mayor concentración cercana al 1.
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.
predictions_validation <- augment(model.1, newdata = dataset.train.validation, type.predict = 'response')
get_prediction_metrics <- function (cutoff, predictions = predictions_validation) {
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')) %>%
mutate(cutoff = cutoff)
}
cutoffs <- seq(0.01,0.95,0.01)
logit_pred <- map_dfr(cutoffs, get_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 y Precision', subtitle= 'Modelo completo', color="")
Elegir un punto de corte y explicar su decisión.
Para decidir el punto de corte de la probabilidad, hay que entender el significado de cada métrica. En términos técnicos, cada una se calcula con las siguientes fórmulas:
\(accuracy = \frac{TP+TN}{TP+FP+FN+TN}\)
\(sensitivity = recall = \frac{TP}{TP+FN}\)
\(specificity = \frac{TN}{TN+FP}\)
\(precision = \frac{TP}{TP+FP}\)
En términos practicos, se puede decir que el accuracy o la exactitud es el porcentaje de individuos correctamente clasificados. La sensibility o recall (sensibilidad) es el porcentaje de individuos clasificados correctamente positivos dentro de los individuos que son realmente positivos, esta métrica penaliza cuando un individuo fue mal clasificado negativo (es decir, fue clasificado negativo, pero realmente es positivo). La specificity o especificidad es el porcentaje de individuos correctamente clasificados negativos dentro de los individuos clasificados negativamente, esta métrica penaliza cuando un indivudo fue mal clasificado positivamente (es decir, fue clasificado como positivo, pero realmente es negativo). La precision o presición es el porcentaje de individuos clasificado correctamente positivos dentro de los individuos que se clasificaron como positivos, esta métrica penaliza a un individuo mal clasificado positivamente (es decir, fue clasificado como positivo, pero realmente es negativo).
La decisión acerca de cual métrica es más importante frente a las demás es una decisión que esta dada completamente por el problema de negocio. Por ejemplo, si estamos en el caso de método de testeo de una enfermedad muy contagiosa, los falsos negativos son muy costosos frente a los falsos positivos, por ende, se eligirá optimizar una métrica que penalice eso, como lo hace el recall. Por otro lado, si se esta midiendo un algoritmo de clasificación de correo basura, podría ser muy costoso clasificar un mail como spam cuando no lo es, en este caso se intentará optimizar una metrica que penalicé eso, como lo es la especificidad. No hay que dejar de lado que algunas métricas, como la exactitud, pueden presentar problemas si el dataset llega a estar desbalanceado.
En el caso particular de este trabajo, ya que el dataset se encuentra balanceado y nuestro objetivo de negocio es sencillamente obtener la clasificación más acertada, se buscará el punto de corte que optimice la exactitud.
logit_pred[logit_pred$term == 'accuracy', ] %>%
arrange(-estimate) %>%
select(c('estimate', 'cutoff')) %>%
head()| estimate | cutoff |
|---|---|
| 0.8052434 | 0.65 |
| 0.8014981 | 0.57 |
| 0.8014981 | 0.59 |
| 0.8014981 | 0.66 |
| 0.8014981 | 0.67 |
| 0.7977528 | 0.56 |
Se eligirá el punto de corte 0.65.
Obtener la matriz de confusión con el modelo y punto de corte elegidos. Interpretarla.
cutoff_point <- 0.65
get_confusion_matrix <- function (predictions, cutoff = cutoff_point) {
table <- predictions %>%
mutate(
predicted_class = if_else(
.fitted > cutoff,
1,
0
) %>%
as.factor(),
Survived = factor(Survived)
)
return(confusionMatrix(table$predicted_class, table$Survived, positive = "1"))
}
get_confusion_matrix(predictions_validation)## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 153 41
## 1 11 62
##
## Accuracy : 0.8052
## 95% CI : (0.7526, 0.851)
## No Information Rate : 0.6142
## P-Value [Acc > NIR] : 1.406e-11
##
## Kappa : 0.5655
##
## Mcnemar's Test P-Value : 5.781e-05
##
## Sensitivity : 0.6019
## Specificity : 0.9329
## Pos Pred Value : 0.8493
## Neg Pred Value : 0.7887
## Prevalence : 0.3858
## Detection Rate : 0.2322
## Detection Prevalence : 0.2734
## Balanced Accuracy : 0.7674
##
## 'Positive' Class : 1
## La matriz de confusión muestra en las columnas la clasificación real de los individuos, mientras que en las filas los individuos clasificados según el modelo. En este caso particular, se observa que de los 164 pasajeros que no sobrevivieron solo 153 fueron clasificados correctamente por el modelo, mientras que de los 103 pasajeros que sobrevivieron solo 62 fueron clasificados correctamente. Esta información se ve resumida en las métricas anteriormente descriptas (accuracy, sensitivity, specificity, precision)
Testeo
Evaluación del modelo
Leer el archivo titanic_complete_test.csv y transformar las variables Survived, Pclass y Embarked a factor.
dataset.test <- read.csv('./titanic_complete_test.csv')
dataset.test.clean <- dataset.train %>%
select(one_of(columns)) %>%
mutate(
Survived = as.factor(Survived),
Pclass = as.factor(Pclass),
Embarked = as.factor(Embarked)
)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).
predictions_test <- augment(model.1, newdata = dataset.test.clean, type.predict = 'response')
get_confusion_matrix(predictions_test)## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 515 128
## 1 34 214
##
## Accuracy : 0.8182
## 95% CI : (0.7913, 0.843)
## No Information Rate : 0.6162
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.5946
##
## Mcnemar's Test P-Value : 2.736e-13
##
## Sensitivity : 0.6257
## Specificity : 0.9381
## Pos Pred Value : 0.8629
## Neg Pred Value : 0.8009
## Prevalence : 0.3838
## Detection Rate : 0.2402
## Detection Prevalence : 0.2783
## Balanced Accuracy : 0.7819
##
## 'Positive' Class : 1
## El modelo presentó resultados similares en el set de validación y en el de testeo. Ambos set de datos obtuvieron un accuracy aproximada de 0.8, sensitivity aproximada de 0.6 y una specificity aproximada de 0.9.