Que se mueran los pobres - Un modelo de regresión logística para predecir la supervivencia en el Titanic
En este post, vamos a crear un modelo de regresión logística para calcular la probabilidad de supervivencia en el Titanic.
Las explicaciones y el código en R están debajo, pero me gustaría comentar algunas cosas previamente.
SOBRE LA IMPUTACIÓN DE DATOS FALTANTES
Una de las principales dificultades del modelo es la cantidad de datos con NA que tenemos. La imputación, o eliminación, de datos faltantes es una de las tareas más importantes en el pre-procesamiento en Machine Learning.
En este caso, la cantidad de veces que aparece Age = NA es demasiado grande para borrarlos. Lo que tenemos que hacer es imputarlos, que puede hacerse de varias formas. La más sencilla es imputalos a la media (datos no sesgados) o a la mediana (datos sesgados). Y una forma más compleja, pero que puede dar mejores resultados es usar regresión logística en otras variables para estimar su valor. Esta última, en una versión sencilla, es la que haremos.
¿CÓMO SE ELIGEN LAS VARIABLES?
De las muchas formas que puede modelizarse el problema, hay que ir eligiendo sobre la marcha qué variables queremos incluir y cuáles excluiremos. Dado que estamos interesados en un modelo predictivo, deberemos tener cuidado de no suprimir variables que resten capacidad de explicar los resultados.
Pese a todo, no toda la información del modelo es útil. El número de billete, por ejemplo, no nos ofrece nada interesante. Como mucho, el lugar donde fue expedido y el precio, pero esta información está contenida en las variables Embark y Pclass respectivamente. De hecho, más adelante veremos que Fare, el precio del billete, no ofrece nada relevante.
LA VALIDACIÓN DEL MODELO
Para poder realizar una predicción con seguridad, deberemos previamente adoptar un protocolo para la validación de los datos. El elegido en este problema se llama hold-out, y consiste en reservar una parte de los datos para entrenar el modelo y otra para validarlos.
Seguiremos esta misma idea tanto en la imputación de Age, para todos esos NA, como en el cálculo de las predicciones para la probabilidad de supervivencia.
EVALUANDO LAS PREDICCIONES
De las formas que hay para evaluar las predicciones hemos elegido la basada en la matriz de confusión. Al final del modelo, calcularemos la matriz de predicción (falsos positivos, verdaderos positivos, falsos negativos y verdaderos negativos) sobre los datos de test y con ello mediremos la capacidad predictiva esperada de nuestro modelo.
Con ayuda de la curva ROC, podremos configurar un valor de corte, que en nuestro caso resultará ser p=0.5, para la probabilidad de que alguien viva o muera. Si p>0.5 determinaremos que esa persona muere y si es menor o igual que sobrevive.
LOS RESULTADOS
El sentido común, la historia y la triste lógica de la vida llevan a pensar que el modelo lo único que hace es confirmar lo obvio: los pobres mueren a una tasa mucho mayor que los ricos. Los hombres mueren más que las mujeres.
Datos del problema
CARGA DE LIBRERÍAS y DATOS
Cargamos las librerias y los datos de entrenamiento
library(corrplot)## corrplot 0.92 loadedlibrary(caret)## Loading required package: ggplot2## Loading required package: latticelibrary(e1071) #normalidad
library(nortest) #coef.asimetría de Fisher
library(dplyr) # manipulacion de datos##
## 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(broom) # resumenes de los modelos
library(visreg) # graficos para regresion logistica
library(margins) # efecto marginal medio
library(rcompanion) # "coeficiente de determinacion"
library(ROCR) # curvas ROC##
## Attaching package: 'ROCR'## The following object is masked from 'package:margins':
##
## predictionlibrary(caTools)
train <- read.csv("train.csv")
test <- read.csv("test.csv")
head(train)## PassengerId Survived Pclass
## 1 1 0 3
## 2 2 1 1
## 3 3 1 3
## 4 4 1 1
## 5 5 0 3
## 6 6 0 3
## Name Sex Age SibSp Parch
## 1 Braund, Mr. Owen Harris 1 22 1 0
## 2 Cumings, Mrs. John Bradley (Florence Briggs Thayer) 0 38 1 0
## 3 Heikkinen, Miss. Laina 0 26 0 0
## 4 Futrelle, Mrs. Jacques Heath (Lily May Peel) 0 35 1 0
## 5 Allen, Mr. William Henry 1 35 0 0
## 6 Moran, Mr. James 1 NA 0 0
## Ticket Fare Cabin Embarked
## 1 A/5 21171 7.2500 S
## 2 PC 17599 71.2833 C85 C
## 3 STON/O2. 3101282 7.9250 S
## 4 113803 53.1000 C123 S
## 5 373450 8.0500 S
## 6 330877 8.4583 Qhead(test)## PassengerId Pclass Name Sex Age
## 1 892 3 Kelly, Mr. James male 34.5
## 2 893 3 Wilkes, Mrs. James (Ellen Needs) female 47.0
## 3 894 2 Myles, Mr. Thomas Francis male 62.0
## 4 895 3 Wirz, Mr. Albert male 27.0
## 5 896 3 Hirvonen, Mrs. Alexander (Helga E Lindqvist) female 22.0
## 6 897 3 Svensson, Mr. Johan Cervin male 14.0
## SibSp Parch Ticket Fare Cabin Embarked
## 1 0 0 330911 7.8292 Q
## 2 1 0 363272 7.0000 S
## 3 0 0 240276 9.6875 Q
## 4 0 0 315154 8.6625 S
## 5 1 1 3101298 12.2875 S
## 6 0 0 7538 9.2250 SPRIMER VISTAZO A LOS DATOS
Significado de cada una de las variables
Vamos a hacernos una idea de lo que contienen las tablas. Esta es una breve explicacion de cada variable:
PassengerId is a unique identifying number assigned to each passenger.
Survived is a flag that indicates if a passenger survived or died ( i.e., 0 = No, 1 = Yes).
Pclass is the passenger class (i.e., 1 = 1st class, 2 = 2nd class, 3 = 3rd class).
Name is the name of the passenger.
Sex indicates the gender of the passenger (i.e., Male or female).
Age indicates the age of the passenger.
Sibsp is the number of siblings/spouses aboard.
Parch is the number of parents/children aboard.
Ticket indicates the ticket number issued to the passenger.
Fare indicates the amount of money spent on their ticket.
Cabin indicates the cabin category occupied by the passenger.
Embarked indicates the port where the passenger embarked from (i.e., C = Cherbourg, Q = Queenstown, S = Southampton)
Survived es la variable que queremos predecir.
Exámen preliminar de la información en la tabla
Lo primero es hacer factors todas las variables categoricas. De otra forma, R se empeña en sacarnos medias y medianas y no tiene sentido (así nos mostrará porcentajes y se entiende mejor)
train$Survived<-as.factor(train$Survived)
train$Pclass<-as.factor(train$Pclass)
train$Sex<-as.factor(train$Sex)
train$SibSp<-as.factor(train$SibSp)
train$Parch<-as.factor(train$Parch)
train$Cabin<-as.factor(train$Cabin)
train$Embarked<-as.factor(train$Embarked)
#test$<-as.factor(test$Survived) #esta no existe
test$Pclass<-as.factor(test$Pclass)
#test$Sex<-as.factor(test$Sex) #esto hay que hacerlo mas adelante
test$SibSp<-as.factor(test$SibSp)
test$Parch<-as.factor(test$Parch)
test$Cabin<-as.factor(test$Cabin)
test$Embarked<-as.factor(test$Embarked)Ahora sí, si hacemos summary nos saldrá algo comprensible:
summary(train)## PassengerId Survived Pclass Name Sex Age
## Min. : 1.0 0:549 1:216 Length:891 0:314 Min. : 0.42
## 1st Qu.:223.5 1:342 2:184 Class :character 1:577 1st Qu.:20.12
## Median :446.0 3:491 Mode :character Median :28.00
## Mean :446.0 Mean :29.70
## 3rd Qu.:668.5 3rd Qu.:38.00
## Max. :891.0 Max. :80.00
## NA's :177
## SibSp Parch Ticket Fare Cabin Embarked
## 0:608 0:678 Length:891 Min. : 0.00 :687 : 2
## 1:209 1:118 Class :character 1st Qu.: 7.91 B96 B98 : 4 C:168
## 2: 28 2: 80 Mode :character Median : 14.45 C23 C25 C27: 4 Q: 77
## 3: 16 3: 5 Mean : 32.20 G6 : 4 S:644
## 4: 18 4: 4 3rd Qu.: 31.00 C22 C26 : 3
## 5: 5 5: 5 Max. :512.33 D : 3
## 8: 7 6: 1 (Other) :186summary(test)## PassengerId Pclass Name Sex Age
## Min. : 892.0 1:107 Length:418 Length:418 Min. : 0.17
## 1st Qu.: 996.2 2: 93 Class :character Class :character 1st Qu.:21.00
## Median :1100.5 3:218 Mode :character Mode :character Median :27.00
## Mean :1100.5 Mean :30.27
## 3rd Qu.:1204.8 3rd Qu.:39.00
## Max. :1309.0 Max. :76.00
## NA's :86
## SibSp Parch Ticket Fare
## 0:283 0 :324 Length:418 Min. : 0.000
## 1:110 1 : 52 Class :character 1st Qu.: 7.896
## 2: 14 2 : 33 Mode :character Median : 14.454
## 3: 4 3 : 3 Mean : 35.627
## 4: 4 4 : 2 3rd Qu.: 31.500
## 5: 1 9 : 2 Max. :512.329
## 8: 2 (Other): 2 NA's :1
## Cabin Embarked
## :327 C:102
## B57 B59 B63 B66: 3 Q: 46
## A34 : 2 S:270
## B45 : 2
## C101 : 2
## C116 : 2
## (Other) : 80Vamos a tener que hacer limpieza de los datos. A simple vista tenemos: - 177 valores de age con NA - 2 valores de embark en blancos - Casi todos los valores de Cabin sin informar - Un blanco en Fare, aparte de los ceros que comentábamos antes.
Afortunadamente, no tenemos que preocuparnos de los outliers. También vemos que sexo en una de las tablas está como 0,1 y en la otra como female, male. Esto habrá que modificarlo.
PRIMERAS CORRECCIONES A LOS DATOS
Male y female en la tabla de test
Como comentábamos antes, Sex toma valores female,male en lugar de 0,1. Tenemos que convertirlos adecuadamente. Male=1, Female=0.
for(i in 1:nrow(test))
{
if (test$Sex[i]=='male') test$Sex[i]<-1
if (test$Sex[i]=='female') test$Sex[i]<-0
}
#ahora si lo convertimos a variable categorica. Si lo hacemos antes da error aqui.
test$Sex<-as.factor(test$Sex)Creación de una tabla conjunta para las imputaciones de Age, Fare,…
#creamos la tabla conjunta con toda la info, salvo la columna de supervivencia
#es necesaria para las imputaciones
AGE <- rbind(test, train[, names(test)])Corrigiendo las variables numericas
También salta a la vista que el tipo de Age es double pero debería ser un entero. Esto no lo vamos a tocar porque más adelante, al imputarle valores, aumentaríamos el error por redondeo, dado que los valores que calculemos serán doubles por usar un pequeño modelo de regresión logística para esto.
Imputando valores a embarked
En embarked tenemos dos valores sin definir - cuidado, no son NA, son blancos. Ambos tienen los mismos valores de cabin, survived y ticked.
which(train$Embarked=="")## [1] 62 830Los pasajeros con numero de ticket parecido a ellos, en su mayoria, embarcaron en S, que es donde los vamos a poner y que coincide con la moda de la variable embarked. Por lo tanto es un valor más que razonable
train$Embarked[62]<-"S"
train$Embarked[830]<-"S"
summary(train$Embarked) #comprobamos que este bien.## C Q S
## 0 168 77 646Analogamente en la tabla AGE las mismas imputaciones (con cuidado, porque el numero de fila en la tabla AGE no tiene porque ser el mismo).
which(AGE$Embarked=="")## [1] 480 1248AGE$Embarked[480 ]<-"S"
AGE$Embarked[1248]<-"S"Imputando valores a FARE
La columna Fare, como puede verse a continuación, contiene valores que son cero. Bien puede ser gente que tenga el ticket gratis por algún motivo o un error en las notas. Como no lo sabemos, y ambas cosas tienen sentido, los dejaremos tal y como están. Si bien ponerlas a la mediana tampoco sería aceptable.
summary(train$Fare)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 7.91 14.45 32.20 31.00 512.33Muchas de las decisiones que tomemos en la imputación serán parecidas a esto: varias opciones igualmente razonables para elegir. Siempre nos movemos en el terreno de la ambigüedad.
Hay un unico valor de fare a NA en toda la tabla, concretamente este, que aparece tambien en la tabla test:
which(is.na(AGE$Fare))## [1] 153Realmente, un valor entre tantos no es una preocupación, y más teniendo en cuenta que Fare no va a ser determinante en los resultados. Como FARE es una variable claramente sesgada lo imputaremos a su mediana
#fare es sesgada a derecha (asimetría positiva)
hist(AGE$Fare)
#imputación a la mediana
AGE$Fare[153]<-median(AGE$Fare[-153])
AGE$Fare[153]## [1] 14.4542#no olvidemos actualizar todas las tablas
which(is.na(train$Fare))## integer(0)test$Fare[153]<-AGE$Fare[153] IMPUTANDO VALORES A AGE
Para imputar valores necesitaremos sacar informacion de los que ya estan informados. Tenemos principalmente tres opciones para imputar Age:
- Imputar a la media o la mediana, según los datos sean o no sesgados.
- Usar un modelo de regresión logística para extraer age del resto de las variables e imputar usando este modelo.
- Ver con qué variables varía age, dividir en grupos e imputar a la media o mediana de cada uno de esos grupos.
La primera alternativa es la más sencilla pero la segunda, al margen de su complejidad, es la que mejores resultados, al menos teoricamente, debería ofrecer. La tercera es un punto intermedio.
AGE_YA_IMPUTADOS<-subset(AGE,is.na(AGE$Age)==FALSE)
AGE_POR_IMPUTAR<-subset(AGE,is.na(AGE$Age)==TRUE)Tenemos 263 valores de Age para imputar, más o menos el 20% del total. Una cantidad más que considerable. Notar que sería un error gravísimo borrar las filas que no tienen Age informado, dado que perderíamos mucha información.
Imputación de Age por la mediana
La variable age es sesgada a derecha, por lo que, como primera aproximacion, podríamos imputar a la mediana si no encontramos nada mejor. Si lo hiciéramos a la media, estaríamos cometiendo un error.
skewness(AGE_YA_IMPUTADOS$Age) #coef. de simetria de fisher## [1] 0.4065061mean(AGE_YA_IMPUTADOS$Age) #media## [1] 29.88114median(AGE_YA_IMPUTADOS$Age) #mediana = 28## [1] 28hist(AGE_YA_IMPUTADOS$Age)
Pero vamos a intentar obtener algo más preciso, dado que tenemos un 20% de edades con NA. Veamos con qué variable categórica varía age.
Imputando Age mediante regresión logística
¿Qué variables incluímos en la imputación de Age?
Con un boxplot por pares, podemos ver que cualquiera de SibSp, Pclass y Parch seguramente sean buenos predictores. Fare no lo será, dado que tiene una correlación muy pequeña con age. Puede verse cómo va variando por grupo:
cor(AGE_YA_IMPUTADOS$Age, AGE_YA_IMPUTADOS$Fare) #correlación AGE, FARE## [1] 0.1775283boxplot(AGE_YA_IMPUTADOS$Age~AGE_YA_IMPUTADOS$SibSp)
boxplot(AGE_YA_IMPUTADOS$Age~AGE_YA_IMPUTADOS$Pclass)
boxplot(AGE_YA_IMPUTADOS$Age~AGE_YA_IMPUTADOS$Parch)
Y de aquí que construyamos este modelo para predecir la edad. De todas ellas, la que mejor parece funcionar es Pclass, que es la que usaremos.
Datos de test y de validación para Age
Vamos dividir los datos de los que conocemos la edad en dos partes. Una de ellas para ajustar el modelo y otra para validarlo.
#make this example reproducible
set.seed(123)
#usamos el 70% de los datos para crear el modelo y el 30% para su validación.
sample <- sample.split(AGE_YA_IMPUTADOS$Age, SplitRatio = 0.7)
AGEtrain <- subset(AGE_YA_IMPUTADOS, sample == TRUE)
AGEtest <- subset(AGE_YA_IMPUTADOS, sample == FALSE)Y ahora sí podemos calcular el modelo con los datos de entrenamiento para age. Usaremos un modelo de regresión logística para ajustar las predicciones.
Regresión logísica para imputar Age
modelo_age<-glm(Age~Pclass,data=AGEtrain)
summary(modelo_age)##
## Call:
## glm(formula = Age ~ Pclass, data = AGEtrain)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -37.12 -8.06 -1.06 7.94 48.94
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.1224 0.9699 40.337 < 2e-16 ***
## Pclass2 -8.9921 1.3752 -6.539 1.16e-10 ***
## Pclass3 -14.0622 1.2033 -11.686 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 180.6119)
##
## Null deviance: 157417 on 737 degrees of freedom
## Residual deviance: 132750 on 735 degrees of freedom
## AIC: 5934.3
##
## Number of Fisher Scoring iterations: 2Validando el modelo
Ahora sí tenemos un buen modelo para predecir Age. Veamos su grado de acierto sobre el conjunto de datos de test:
predictions<-predict(modelo_age,newdata=AGEtrain,type="response")
#calculamos la raiz del error cuadratico medio de las predicciones
predictions<-as.double(predictions)
mean(predictions-AGEtrain$Age)## [1] 7.675411e-14sd(predictions-AGEtrain$Age)## [1] 13.42094summary(predictions)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 25.06 25.06 30.13 30.02 39.12 39.12Sin ser excepcional, el resultado no es malo y desde luego será mejor que imputar a la mediana y olvidarnos del asunto. Ya podemos calcular la edad para los datos faltantes.
PREDICTED_AGES<-predict(modelo_age,newdata=AGE_POR_IMPUTAR,type="response")
summary(PREDICTED_AGES)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 25.06 25.06 25.06 27.45 25.06 39.12Cargando la variable Age calculada
Y metemos los datos
for(i in 1:nrow(AGE_POR_IMPUTAR)) AGE_POR_IMPUTAR$Age[i]<-PREDICTED_AGES[i]
head(AGE_POR_IMPUTAR$Age)## [1] 25.06017 39.12240 25.06017 25.06017 25.06017 25.06017La imputación de la edad ha funcionado estupendamente. Volcamos los datos en las tablas del problema original y podemos pasar a la regresión logística para la supervivencia.
No vamos a redondear las edades imputadas, ya que perderíamos información.
for (i in 1:nrow(AGE_POR_IMPUTAR))
{
fila <- AGE_POR_IMPUTAR$PassengerId[i];
if(fila<=891)
train$Age[fila]<-AGE_POR_IMPUTAR$Age[i]
else
test$Age[fila-891]<-AGE_POR_IMPUTAR$Age[i]
}Ahora ya tanto la tabla train como test tienen age imputado.
summary(train$Age)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.42 22.00 26.00 29.31 37.00 80.00summary(test$Age)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.17 23.00 25.06 29.56 36.38 76.00El modelo de regresión logística
¿Que tipo de resultados esperamos?
De un simple vistazo a la tabla vemos que la probabilidad de sobrevivir depende claramente del precio del billete que compres:
TABLA<-table(train$Survived,train$Pclass)
TABLA##
## 1 2 3
## 0 80 97 372
## 1 136 87 119El porcentaje de supervivientes de primera clase es mucho mayor que el de tercera. Y esto no puede ser casualidad. Vamos a confirmar nuestra intuición con un test de independencia chi-cuadrado:
chisq.test(TABLA)##
## Pearson's Chi-squared test
##
## data: TABLA
## X-squared = 102.89, df = 2, p-value < 2.2e-16Por lo tanto está claro que deberemos incluir el precio del billete en nuestro modelo predictivo. Los pasajeros de primera clase sobrevivieron mucho más que los de tercera.
Variables que NO deberíamos incluir
No todas las variables son relevantes en el modelo. El hecho de incluirlas, puede resultar en p-valores pequeños y un aumento artificial del ajuste del modelo, que llevará a malos resultados.
No vamos a incluir cabin, porque tiene muy pocos valores informados, y tampoco el nombre - por no tener ninguna relación con lo demás.
El modelo de regresión logística para la supervivencia
modelo_Survived <- glm(Survived ~ Pclass + Sex + Age + SibSp + Parch + Fare + Embarked, data=train, family="binomial")
summary(modelo_Survived)##
## Call:
## glm(formula = Survived ~ Pclass + Sex + Age + SibSp + Parch +
## Fare + Embarked, family = "binomial", data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8248 -0.6164 -0.4157 0.5810 2.4880
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.920e+00 5.233e-01 7.491 6.84e-14 ***
## Pclass2 -1.105e+00 3.103e-01 -3.561 0.000369 ***
## Pclass3 -2.213e+00 3.173e-01 -6.974 3.07e-12 ***
## Sex1 -2.687e+00 2.038e-01 -13.182 < 2e-16 ***
## Age -4.003e-02 8.765e-03 -4.566 4.96e-06 ***
## SibSp1 7.171e-02 2.247e-01 0.319 0.749659
## SibSp2 -3.008e-01 5.375e-01 -0.560 0.575752
## SibSp3 -2.248e+00 7.216e-01 -3.115 0.001839 **
## SibSp4 -1.687e+00 7.624e-01 -2.213 0.026894 *
## SibSp5 -1.595e+01 9.584e+02 -0.017 0.986720
## SibSp8 -1.604e+01 7.573e+02 -0.021 0.983105
## Parch1 3.650e-01 2.901e-01 1.258 0.208361
## Parch2 2.517e-02 3.830e-01 0.066 0.947602
## Parch3 3.603e-01 1.057e+00 0.341 0.733099
## Parch4 -1.587e+01 1.054e+03 -0.015 0.987987
## Parch5 -1.166e+00 1.172e+00 -0.995 0.319750
## Parch6 -1.644e+01 2.400e+03 -0.007 0.994534
## Fare 2.101e-03 2.499e-03 0.841 0.400411
## EmbarkedQ 3.122e-02 3.869e-01 0.081 0.935673
## EmbarkedS -2.845e-01 2.450e-01 -1.161 0.245507
## ---
## 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.58 on 871 degrees of freedom
## AIC: 803.58
##
## Number of Fisher Scoring iterations: 15Vamos poco a poco. Primero quitaremos embarked
modelo_Survived <- glm(Survived ~ Pclass + Sex + Age + SibSp + Parch + Fare, data=train, family="binomial")
summary(modelo_Survived)##
## Call:
## glm(formula = Survived ~ Pclass + Sex + Age + SibSp + Parch +
## Fare, family = "binomial", data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8825 -0.6188 -0.4193 0.5956 2.4577
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.772e+00 4.957e-01 7.611 2.73e-14 ***
## Pclass2 -1.173e+00 3.060e-01 -3.833 0.000127 ***
## Pclass3 -2.209e+00 3.116e-01 -7.090 1.35e-12 ***
## Sex1 -2.720e+00 2.015e-01 -13.499 < 2e-16 ***
## Age -4.077e-02 8.767e-03 -4.650 3.32e-06 ***
## SibSp1 6.213e-02 2.246e-01 0.277 0.782068
## SibSp2 -3.141e-01 5.362e-01 -0.586 0.558072
## SibSp3 -2.392e+00 7.187e-01 -3.328 0.000875 ***
## SibSp4 -1.766e+00 7.669e-01 -2.302 0.021322 *
## SibSp5 -1.606e+01 9.565e+02 -0.017 0.986601
## SibSp8 -1.615e+01 7.548e+02 -0.021 0.982926
## Parch1 3.693e-01 2.897e-01 1.275 0.202426
## Parch2 -2.153e-03 3.814e-01 -0.006 0.995497
## Parch3 3.311e-01 1.043e+00 0.318 0.750846
## Parch4 -1.599e+01 1.052e+03 -0.015 0.987870
## Parch5 -1.207e+00 1.170e+00 -1.032 0.302273
## Parch6 -1.656e+01 2.400e+03 -0.007 0.994494
## Fare 2.558e-03 2.480e-03 1.032 0.302221
## ---
## 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: 765.54 on 873 degrees of freedom
## AIC: 801.54
##
## Number of Fisher Scoring iterations: 15Y ahora es el turno de Fare
modelo_Survived <- glm(Survived ~ Pclass + Sex + Age + SibSp + Parch, data=train, family="binomial")
summary(modelo_Survived)##
## Call:
## glm(formula = Survived ~ Pclass + Sex + Age + SibSp + Parch,
## family = "binomial", data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8480 -0.5989 -0.4173 0.5896 2.4673
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.975e+00 4.571e-01 8.698 < 2e-16 ***
## Pclass2 -1.318e+00 2.735e-01 -4.821 1.43e-06 ***
## Pclass3 -2.379e+00 2.671e-01 -8.909 < 2e-16 ***
## Sex1 -2.729e+00 2.011e-01 -13.572 < 2e-16 ***
## Age -4.139e-02 8.739e-03 -4.736 2.18e-06 ***
## SibSp1 7.327e-02 2.238e-01 0.327 0.743376
## SibSp2 -2.715e-01 5.383e-01 -0.504 0.613980
## SibSp3 -2.300e+00 6.984e-01 -3.292 0.000994 ***
## SibSp4 -1.781e+00 7.697e-01 -2.314 0.020692 *
## SibSp5 -1.604e+01 9.561e+02 -0.017 0.986616
## SibSp8 -1.606e+01 7.539e+02 -0.021 0.983000
## Parch1 4.090e-01 2.862e-01 1.429 0.152981
## Parch2 7.202e-02 3.746e-01 0.192 0.847532
## Parch3 3.679e-01 1.048e+00 0.351 0.725493
## Parch4 -1.587e+01 1.040e+03 -0.015 0.987827
## Parch5 -1.138e+00 1.168e+00 -0.974 0.330205
## Parch6 -1.646e+01 2.400e+03 -0.007 0.994528
## ---
## 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: 766.71 on 874 degrees of freedom
## AIC: 800.71
##
## Number of Fisher Scoring iterations: 15Claramente Parch debe ser eliminada, dado que no es significativa.
modelo_Survived <- glm(Survived ~ Pclass + Sex + Age + SibSp, data=train, family="binomial")
summary(modelo_Survived)##
## Call:
## glm(formula = Survived ~ Pclass + Sex + Age + SibSp, family = "binomial",
## data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.9022 -0.5985 -0.4180 0.6187 2.5000
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.16757 0.43315 9.621 < 2e-16 ***
## Pclass2 -1.33204 0.27084 -4.918 8.74e-07 ***
## Pclass3 -2.47836 0.26104 -9.494 < 2e-16 ***
## Sex1 -2.70949 0.19615 -13.813 < 2e-16 ***
## Age -0.04577 0.00837 -5.469 4.53e-08 ***
## SibSp1 0.12033 0.21109 0.570 0.56864
## SibSp2 -0.18292 0.52249 -0.350 0.72626
## SibSp3 -2.14936 0.68939 -3.118 0.00182 **
## SibSp4 -1.69904 0.74430 -2.283 0.02245 *
## SibSp5 -16.01036 958.62424 -0.017 0.98667
## SibSp8 -15.97828 754.65794 -0.021 0.98311
## ---
## 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: 774.98 on 880 degrees of freedom
## AIC: 796.98
##
## Number of Fisher Scoring iterations: 15Sisbp aparentemente contiene información relevante. Como estamos en un modelo predictivo, no nos conviene eliminar variables salvo que esté claro que así debe hacerse. Dado que cada variable que eliminamos hace que se reduzca la capacidad de explicar la variación en los datos.
Predicciones del modelo
Veamos que tal funciona el modelo sobre los datos de test. Tomaremos p>0.5 como sobrevivir=TRUE y p<0.5 como sobrevivir=FALSE.
predictions<-predict(modelo_Survived,newdata=train,type="response")
predictions <- ifelse(predictions > 0.5,1,0)
predictions<-as.factor(predictions)Validación del modelo
Dos de los métodos más usados para validar los modelos de regresión logística son la matriz de confusión, con sus indicadores asociados (sensibilidad, especificidad, precisión, F1-score) y las curvas ROC, que proporcionan una forma gráfica de ver qué bien ajusta el modelo.
Matriz de confusión
Creamos la matriz de confusion y las estadisticas para las predicciones.
## Confusion matrix and statistics
library(caret)
confusionMatrix(data=predictions, reference=train$Survived)## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 480 93
## 1 69 249
##
## Accuracy : 0.8182
## 95% CI : (0.7913, 0.843)
## No Information Rate : 0.6162
## P-Value [Acc > NIR] : < 2e-16
##
## Kappa : 0.6105
##
## Mcnemar's Test P-Value : 0.07075
##
## Sensitivity : 0.8743
## Specificity : 0.7281
## Pos Pred Value : 0.8377
## Neg Pred Value : 0.7830
## Prevalence : 0.6162
## Detection Rate : 0.5387
## Detection Prevalence : 0.6431
## Balanced Accuracy : 0.8012
##
## 'Positive' Class : 0
## Podemos ver que la precisión del modelo es bastante alta, un 81%. También los valores de la sensibilidad y especificidad son muy buenos.
Curva ROC y área bajo la curva
Ahora optimizaremos el modelo calculando la curva ROC
## ROC Curve and calculating the area under the curve(AUC)
library(ROCR)
predictions <- predict(modelo_Survived, newdata=train, type="response")
ROCRpred <- prediction(predictions, train$Survived)
ROCRperf <- performance(ROCRpred, measure = "tpr", x.measure = "fpr")
plot(ROCRperf, colorize = TRUE, text.adj = c(-0.2,1.7), print.cutoffs.at = seq(0,1,0.1))
El area bajo la curva es la medida del performance del test
auc <- performance(ROCRpred, measure = "auc")
auc <- auc@y.values[[1]]
auc## [1] 0.8580726La medida del área bajo la curva ROC nos dice que el modelo, con un 0.85, es francamente bueno. A partir de 0.9 se consideraría excepcional.
Observar también que el valor 0.5, como valor de corte, es un buen punto intermedio en cuanto a la sensibilidad/especifidad se refiere.
PREDICCIONES DE SUPERVIVENCIA
Veamos quién vive y quién muere según nuestro modelo:
test<-cbind(test,Survived=0)
#test$Survived<-as.factor(test$Survived)
predictions<-predict(modelo_Survived,newdata=test,type="response")
summary(predictions)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.1073 0.2964 0.4026 0.6627 0.9696predictions <- ifelse(predictions > 0.5,1,0) #metemos el valor de corte calculado con la ROCFinalmente, guardamos en un fichero las predicciones
for(i in 1:nrow(test))
test$Survived[i]<-predictions[i]Y vemos cuántos sobreviven
table(test$Survived)##
## 0 1
## 258 160CONCLUSIONES
En este trabajo, hemos creado un modelo de regresión logística relativamente sencillo para calcular la probabilidad de supervivencia en el hundimiento del Titanic.
La tarea que más trabajo ha sido la imputación a la edad, que podíamos haberla hecho de varias formas posibles.
como suele pasar en la regresión logística, vamos navegando por la ambigüedad de la elección de variables y diversos métodos de imputación. Así pues, según lo que nos pregunten, deberemos elegir una forma de trabajo adecuada.
Comentarios
Publicar un comentario