Logistic Regression
Cristian Ramos Lorenzo
29/5/2019
1. Introducción.
En esta sección trabajaremos con el método de clasificación de Regresión Logística sobre el conocido dataset del Titanic de Kaggle. El objetivo es trabajar sobre una variable discreta, normalmente se relaciona con grupos binarios (es decir, 0 y 1).
Si usásemos el modelo de regresión lineal para intentar predecir en grupos binarios obtendríamos un modelo erróneo, es por ello que lo mejor es usar la regresión logística, que junto a la función sigmoide nos dará una curva continua de nuestros datos para predecir la probabilidad. La función sigmoide viene dada por la siguiente fórmula:
\[P(t)=\dfrac{1}{1+e^{-t}}\] Podemos apreciar la diferencia entre los dos modelos comentados previamente en la siguiente imagen:
One
Supongamos que el eje de las Y representa la probabilidad de la devolución de un préstamo. Con el modelo lineal \[y=b_{0}+b_{1}·x\] no logramos obtener los resultados correctamente, sin embargo, si pasamos el modelo lineal como la variable t en la curva sigmoide, obtenemos una transformación del modelo en la curva de la derecha.
\[p=\dfrac{1}{1+e^{-(b_{0}+b{1}·x)}}\]
Se toma el valor 0.5 como punto de corte para determinar si un valor pertenece a las clase 0 o a la clase 1.
2. Descarga y estructura del dataset.
setwd("~/R/R-Course-HTML-Notes/R-for-Data-Science-and-Machine-Learning/Machine Learning with R")
df.train <- read.csv("titanic_train.csv")
str(df.train)## 'data.frame': 891 obs. of 12 variables:
## $ PassengerId: int 1 2 3 4 5 6 7 8 9 10 ...
## $ Survived : int 0 1 1 1 0 0 0 0 1 1 ...
## $ Pclass : int 3 1 3 1 3 3 1 3 3 2 ...
## $ Name : Factor w/ 891 levels "Abbing, Mr. Anthony",..: 109 191 358 277 16 559 520 629 417 581 ...
## $ Sex : Factor w/ 2 levels "female","male": 2 1 1 1 2 2 2 2 1 1 ...
## $ Age : num 22 38 26 35 35 NA 54 2 27 14 ...
## $ SibSp : int 1 1 0 1 0 0 0 3 0 1 ...
## $ Parch : int 0 0 0 0 0 0 0 1 2 0 ...
## $ Ticket : Factor w/ 681 levels "110152","110413",..: 524 597 670 50 473 276 86 396 345 133 ...
## $ Fare : num 7.25 71.28 7.92 53.1 8.05 ...
## $ Cabin : Factor w/ 148 levels "","A10","A14",..: 1 83 1 57 1 1 131 1 1 1 ...
## $ Embarked : Factor w/ 4 levels "","C","Q","S": 4 2 4 4 4 3 4 4 4 2 ...3. Visualización de NA’s y del dataset.
library(Amelia)## Loading required package: Rcpp## ##
## ## Amelia II: Multiple Imputation
## ## (Version 1.7.5, built: 2018-05-07)
## ## Copyright (C) 2005-2019 James Honaker, Gary King and Matthew Blackwell
## ## Refer to http://gking.harvard.edu/amelia/ for more information
## ##missmap(df.train, main='Missing Map', col=c('yellow','black'), legend=FALSE)
library(ggplot2)## Registered S3 methods overwritten by 'ggplot2':
## method from
## [.quosures rlang
## c.quosures rlang
## print.quosures rlangVeamos la proporción de pasajeros que sobrevivieron frente a los que no.
ggplot(df.train, aes(Survived)) + geom_bar()
Ahora, veamos la distribución de la clase según el ticket.
ggplot(df.train, aes(Pclass)) + geom_bar(aes(fill=factor(Pclass))) + scale_fill_discrete(name="Ticket class")
Veamos la distribución de pasajeros según su sexo.
ggplot(df.train, aes(Sex)) + geom_bar(aes(fill=factor(Sex))) + scale_fill_discrete(name="Sex")
Veamos ahora en función de su edad.
ggplot(df.train,aes(Age)) + geom_histogram(bins=20, alpha=0.5, fill="blue")## Warning: Removed 177 rows containing non-finite values (stat_bin).
La mayoría de gente a bordo era gente entre los 20 y 40 años de edad.
ggplot(df.train, aes(SibSp)) + geom_bar()
Los pasajeros que más frecuentaban el barco en el viaje eran gente que iba sola a bordo.
ggplot(df.train, aes(Fare)) + geom_histogram(fill="green", col="black", alpha=0.5)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Tal y como era esperado, la gente había comprado el ticket más barato, lo que tiene sentido pues el mayor número de pasajeros iban en tercera clase.
Veamos ahora la cantidad de gente que no sobrevivió en función de la clase en la que viajaban.
ggplot(df.train, aes(Survived, fill=factor(Pclass))) + geom_bar(position="dodge") + scale_fill_discrete(name="Passenger's Class")
Podemos observar que el mayor número de gente que no sobrevivió al hundimiento eran viajeros de la tercera clase.
Por último veamos en función de su sexo.
ggplot(df.train, aes(Survived, fill=factor(Sex))) + geom_bar(position="dodge") + scale_fill_discrete(name="Passenger's Sex")
Un mayor número de hombres fallecieron en el hundimiento.
4. Limpieza de valores NA.
En muchas ocasiones lo mejor ante este tipo de situaciones es eliminar las observaciones con NA’s, pero en este caso hay una cantidad muy alta de observaciones con esa condición.
table(is.na(df.train))##
## FALSE TRUE
## 10515 177Una de las soluciones típicas es usar la media de las edades de los pasajeros, pero en este caso, vamos a usar la media de los pasajeros en función de su clase. Podemos hacerlo numéricamente:
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, uniondf.train %>% filter(Age!='Na') %>% group_by(Pclass) %>% summarize(media=median(Age))## # A tibble: 3 x 2
## Pclass media
## <int> <dbl>
## 1 1 37
## 2 2 29
## 3 3 24O verlo gráficamente:
ggplot(df.train, aes(Pclass, Age)) + geom_boxplot(aes(group=Pclass, fill=factor(Pclass)), alpha=0.4) + scale_fill_discrete(name="Passenger's class") + scale_y_continuous(breaks= seq(0,70, by=2)) + theme_bw()## Warning: Removed 177 rows containing non-finite values (stat_boxplot).
Vamos a modificar la columna edad para esas observaciones:
fix_age <- function(age, class){
out <- age
for (i in 1:length(age)){
if (is.na(age[i])){
if (class[i]==1){
out[i] <- 37
}
else if (class[i]==2){
out[i] <- 29
}
else{
out[i] <- 24
}
}
else{
out[i] <- age[i]
}
}
return(out)
}
fixed.ages <- fix_age(df.train$Age, df.train$Pclass)
df.train$Age <- fixed.ages5. Preparación del dataset para aplicar el modelo.
Comprobamos que no existan NA’s en el dataset.
missmap(df.train, main="Missing Map", col=c("yellow","black"), legend=FALSE)
Eliminamos las variables con demasiados niveles pues no las vamos a usar.
df.train <- df.train %>% select(-PassengerId, -Ticket,-Cabin, -Name)
head(df.train)## Survived Pclass Sex Age SibSp Parch Fare Embarked
## 1 0 3 male 22 1 0 7.2500 S
## 2 1 1 female 38 1 0 71.2833 C
## 3 1 3 female 26 0 0 7.9250 S
## 4 1 1 female 35 1 0 53.1000 S
## 5 0 3 male 35 0 0 8.0500 S
## 6 0 3 male 24 0 0 8.4583 QCambiamos el formato de las columnas, aunque la columna SibSp se puede interpretar como una variable continua, en realidad el número de hermanos y esposos va a estar “limitado”. Es por ello que mejor es tratarla como una variable discreta.
df.train$Survived <- factor(df.train$Survived)
df.train$Pclass <- factor(df.train$Pclass)
df.train$SibSp <- factor(df.train$SibSp)
df.train$Parch <- factor(df.train$Parch)
str(df.train)## 'data.frame': 891 obs. of 8 variables:
## $ 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 24 54 2 27 14 ...
## $ SibSp : Factor w/ 7 levels "0","1","2","3",..: 2 2 1 2 1 1 1 4 1 2 ...
## $ Parch : Factor w/ 7 levels "0","1","2","3",..: 1 1 1 1 1 1 1 2 3 1 ...
## $ Fare : num 7.25 71.28 7.92 53.1 8.05 ...
## $ Embarked: Factor w/ 4 levels "","C","Q","S": 4 2 4 4 4 3 4 4 4 2 ...6. Entrenamiento del modelo
log.model <- glm(Survived~., family=binomial(link='logit'), data=df.train)
summary(log.model)##
## Call:
## glm(formula = Survived ~ ., family = binomial(link = "logit"),
## data = df.train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8158 -0.6134 -0.4138 0.5808 2.4896
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.845e+01 1.660e+03 0.011 0.991134
## Pclass2 -1.079e+00 3.092e-01 -3.490 0.000484 ***
## Pclass3 -2.191e+00 3.161e-01 -6.930 4.20e-12 ***
## Sexmale -2.677e+00 2.040e-01 -13.123 < 2e-16 ***
## Age -3.971e-02 8.758e-03 -4.534 5.79e-06 ***
## SibSp1 8.135e-02 2.245e-01 0.362 0.717133
## SibSp2 -2.897e-01 5.368e-01 -0.540 0.589361
## SibSp3 -2.241e+00 7.202e-01 -3.111 0.001862 **
## SibSp4 -1.675e+00 7.620e-01 -2.198 0.027954 *
## SibSp5 -1.595e+01 9.588e+02 -0.017 0.986731
## SibSp8 -1.607e+01 7.578e+02 -0.021 0.983077
## Parch1 3.741e-01 2.895e-01 1.292 0.196213
## Parch2 3.862e-02 3.824e-01 0.101 0.919560
## Parch3 3.655e-01 1.056e+00 0.346 0.729318
## Parch4 -1.586e+01 1.055e+03 -0.015 0.988007
## Parch5 -1.152e+00 1.172e+00 -0.983 0.325771
## Parch6 -1.643e+01 2.400e+03 -0.007 0.994536
## Fare 2.109e-03 2.490e-03 0.847 0.397036
## EmbarkedC -1.458e+01 1.660e+03 -0.009 0.992995
## EmbarkedQ -1.456e+01 1.660e+03 -0.009 0.993001
## EmbarkedS -1.486e+01 1.660e+03 -0.009 0.992857
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1186.66 on 890 degrees of freedom
## Residual deviance: 763.41 on 870 degrees of freedom
## AIC: 805.41
##
## Number of Fisher Scoring iterations: 15Abrimos el dataset para testear el modelo.
setwd("~/R/R-Course-HTML-Notes/R-for-Data-Science-and-Machine-Learning/Machine Learning with R")
df.test <- read.csv('titanic_test.csv')7. Preparación dataset para testear el modelo
Al igual que ocurría con el dataset para entrenar el modelo, tenemos que modificar los valores NA’s que se encuentran en la columna Age.
ggplot(df.test, aes(Pclass, Age)) + geom_boxplot(aes(group=Pclass, fill=factor(Pclass)), alpha=0.4) + scale_fill_discrete(name="Passenger's class") + scale_y_continuous(breaks= seq(0,70, by=2)) + theme_bw()## Warning: Removed 86 rows containing non-finite values (stat_boxplot).
df.test %>% filter(Age!='Na') %>% group_by(Pclass) %>% summarize(media=median(Age))## # A tibble: 3 x 2
## Pclass media
## <int> <dbl>
## 1 1 42
## 2 2 26.5
## 3 3 24fix_age_test <- function(age, class){
out <- age
for (i in 1:length(age)){
if (is.na(age[i])){
if (class[i]==1){
out[i] <- 42
}
else if (class[i]==2){
out[i] <- 26.5
}
else{
out[i] <- 24
}
}
else{
out[i] <- age[i]
}
}
return(out)
}
fixed.ages.test <- fix_age_test(df.test$Age, df.test$Pclass)
df.test$Age <- fixed.ages.testdf.test$Pclass <- factor(df.test$Pclass)
df.test$SibSp <- factor(df.test$SibSp)
df.test$Parch <- factor(df.test$Parch)
str(df.test)## 'data.frame': 418 obs. of 11 variables:
## $ PassengerId: int 892 893 894 895 896 897 898 899 900 901 ...
## $ Pclass : Factor w/ 3 levels "1","2","3": 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 : Factor w/ 7 levels "0","1","2","3",..: 1 2 1 1 2 1 1 2 1 3 ...
## $ Parch : Factor w/ 8 levels "0","1","2","3",..: 1 1 1 1 2 1 1 2 1 1 ...
## $ 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/ 77 levels "","A11","A18",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Embarked : Factor w/ 3 levels "C","Q","S": 2 3 2 3 3 3 2 3 1 3 ...library(plyr)## -------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## -------------------------------------------------------------------------##
## Attaching package: 'plyr'## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarizedf.test$Parch <- revalue(df.test$Parch,c("9" = "6"))8. Ejecución del modelo sobre el dataset
fitted.probabilities <- predict(log.model, df.test, type='response')A continuación tenemos que dividir esas probabilidades en aquellas mayores que 0.5y las que son menores, para incluír al pasajero en el grupo de supervivientes o en otro.
fitted.results <- ifelse(fitted.probabilities>0.5, 1,0)
df.test$Survived_predict <- fitted.resultsComprobación de la eficacia del modelo
Esta columna sería la correspondiente para saber si los pasajeros sobrevivirían o no aplicando el método de regresión logística. No podemos saber la eficacia del modelo pues en el dataset de test no existe una columna denominada Survived. Sin embargo, para realizar esa comprobación podemos usar el dataset df.train y tomar el 70% de sus observaciones para ver la eficacia del modelo.
library(caTools)
sample <- sample.split(df.train, 0.7)
train_final <- subset(df.train, sample=TRUE)
test_final <- subset(df.train, sample=FALSE)Entrenamos el modelo sobre el dataset de entrenamiento:
final_log_model <- glm(Survived ~., family=binomial(link='logit'),train_final)
summary(final_log_model)##
## Call:
## glm(formula = Survived ~ ., family = binomial(link = "logit"),
## data = train_final)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8158 -0.6134 -0.4138 0.5808 2.4896
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.845e+01 1.660e+03 0.011 0.991134
## Pclass2 -1.079e+00 3.092e-01 -3.490 0.000484 ***
## Pclass3 -2.191e+00 3.161e-01 -6.930 4.20e-12 ***
## Sexmale -2.677e+00 2.040e-01 -13.123 < 2e-16 ***
## Age -3.971e-02 8.758e-03 -4.534 5.79e-06 ***
## SibSp1 8.135e-02 2.245e-01 0.362 0.717133
## SibSp2 -2.897e-01 5.368e-01 -0.540 0.589361
## SibSp3 -2.241e+00 7.202e-01 -3.111 0.001862 **
## SibSp4 -1.675e+00 7.620e-01 -2.198 0.027954 *
## SibSp5 -1.595e+01 9.588e+02 -0.017 0.986731
## SibSp8 -1.607e+01 7.578e+02 -0.021 0.983077
## Parch1 3.741e-01 2.895e-01 1.292 0.196213
## Parch2 3.862e-02 3.824e-01 0.101 0.919560
## Parch3 3.655e-01 1.056e+00 0.346 0.729318
## Parch4 -1.586e+01 1.055e+03 -0.015 0.988007
## Parch5 -1.152e+00 1.172e+00 -0.983 0.325771
## Parch6 -1.643e+01 2.400e+03 -0.007 0.994536
## Fare 2.109e-03 2.490e-03 0.847 0.397036
## EmbarkedC -1.458e+01 1.660e+03 -0.009 0.992995
## EmbarkedQ -1.456e+01 1.660e+03 -0.009 0.993001
## EmbarkedS -1.486e+01 1.660e+03 -0.009 0.992857
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1186.66 on 890 degrees of freedom
## Residual deviance: 763.41 on 870 degrees of freedom
## AIC: 805.41
##
## Number of Fisher Scoring iterations: 15Igual que antes podemos observar que tanto pertenecer a la clase 2 como a la clase 3 parece estar relacionado con la variable survived, así como ser hombre. Esto tiene sentido viendo los gráficos de la sección 3.
Aplicamos el método al dataset para testear.
fitted.probabilities.final <- predict(final_log_model, test_final, type='response')
fitted.results.final <- ifelse(fitted.probabilities.final>0.5, 1,0)Calculamos ahora la eficacia del modelo.
1-mean(fitted.results.final != test_final$Survived)## [1] 0.8181818Ha acertado en un 81% de los datos!
Creamos, finalmente, una matriz de confusión.
table(test_final$Survived, fitted.results.final)## fitted.results.final
## 0 1
## 0 482 67
## 1 95 247