En esta primera entrada voy a centrarme en una de las múltiples bases de datos del paquete Datarium. Este paquete constituye un excelente recurso para aplicar las técnicas estadísticas que voy aprendiendo. He elegido el dataset titanic.raw. A continuación vemos la estructura del dataset:
## 'data.frame': 2201 obs. of 4 variables:
## $ Class : Factor w/ 4 levels "1st","2nd","3rd",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ Sex : Factor w/ 2 levels "Male","Female": 1 1 1 1 1 1 1 1 1 1 ...
## $ Age : Factor w/ 2 levels "Child","Adult": 1 1 1 1 1 1 1 1 1 1 ...
## $ Survived: Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
Estamos ante un dataset con 2201 observaciones y 4 variables, todas ellas factores. Dichas variables recogen información básica sobre los 2201 pasajeros del Titánic. A diferencia de otros datasets sobre el hundimiento del Titánic en el que hay más variables cuantitativas, éste dataset unicamente presenta variables de tipo categórico.
Veamos en primer lugar si el dataset está completo o tenemos casos perdidos:
El gráfico muestra que los datos están muy bien estructurados, sin missing data. Esto nos va a facilitar parte del trabajo.
Comenzamos pues analizando cómo se distribuyen los niveles de los diferentes factores:
Destaca que la mayoría de los pasajeros formaban parte de la tripulación (Crew), seguido por los de tercera clase. En cuanto al sexo, claramente había mayoría de varones. Vemos también que casi toda la tripulación era adulta (tan solo aparecen 109 niños). En cuanto a los supervivientes, más de la mitad (67%) no sobrevivieron al naugrafio:
| No | Yes |
|---|---|
| 0.676965 | 0.323035 |
Bien, vamos a sacar un poco más de punta a estos datos. Nos interesa saber la relación de las variables Class, Sex y Age con la variable Survived. El siguiente gráfico de barras apiladas nos puede ayudar:
Así, en la variable Class, comprobamos que tanto la tripulación como los pasajeros de tercera clase presentaron una mayor mortalidad. De hecho es llamativo que entre los de primera clase sobrevivieron más que murieron (203 vs 122). En cuanto al sexo, la mayoría de los no supervivientes son varones. De la misma forma, la mortalidad entre los adultos fue muy superior a los niños.
A continuación muestro las tablas de frecuencias con el test Chi2 correspondiente:
| Class/Survived | No | Yes |
|---|---|---|
| 1st | 122 | 203 |
| 2nd | 167 | 118 |
| 3rd | 528 | 178 |
| Crew | 673 | 212 |
##
## Pearson's Chi-squared test
##
## data: t1
## X-squared = 190.4, df = 3, p-value < 2.2e-16
| Sex/Survived | No | Yes |
|---|---|---|
| Male | 1364 | 367 |
| Female | 126 | 344 |
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: t2
## X-squared = 454.5, df = 1, p-value < 2.2e-16
| Age/Survived | No | Yes |
|---|---|---|
| Child | 52 | 57 |
| Adult | 1438 | 654 |
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: t3
## X-squared = 20.005, df = 1, p-value = 7.725e-06
Vemos, efectivamente, que las diferencias entre las columnas y filas de las tablas son significativas.
Ahora podemos plantear un modelo predictivo mediante una regresión logística. De esta forma intentaremos predecir con nuestro modelo si un determinado pasajero es clasificado correctamente como superviviente o como fallecido en el hundimiento. Para empezar selecciono al azar un subgrupo train con el 80% de los pasajeros que constituirá el grupo de entrenamiento del modelo. Posteriormente validaremos el modelo en el subgrupo test formado por el 20% restante de los pasajeros. A continuación muestro el test de la razón de verosimilitud, al comparar el modelo frente al modelo nulo. Comprobamos que nuestro modelo es significativo:
| Resid. Df | Resid. Dev | Df | Deviance | Pr(>Chi) |
|---|---|---|---|---|
| 1760 | 2199.433 | NA | NA | NA |
| 1755 | 1721.238 | 5 | 478.1948 | 0 |
Ahora paso a validar el modelo. Analizo la bondad de ajuste mediante el test de Hosmer-Lemeshow:
## Warning in Ops.factor(1, y): '-' not meaningful for factors
##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: train$Survived, fitted(model)
## X-squared = 1761, df = 7, p-value < 2.2e-16
Comprobamos que el resultado es significativo, así pues el ajuste que conseguimos con nuestro modelo no es demasiado bueno. Voy a analizar a continuación el porcentaje de la incertidumbre que nuestro modelo explica en el grupo train. Para ello calculo la R2 de McFadden:
## fitting null model for pseudo-r2
## McFadden
## 0.2174173
De nuevo se confirma las limitaciones del modelo, ya que tan solo es capaz de explicar un 22% de la variabilidad.
A continuación muestro el porcentaje de pasajeros clasificados correctamente por el modelo, en el grupo test:
## [1] 0.7386364
Vemos que nuestro modelo es capaz de clasificar correctamente al 74% de los pasajeros.
La tabla de confusión en el grupo test sería:
| Asignacion/Sobrevive | No | Yes |
|---|---|---|
| No | 258 | 86 |
| Yes | 29 | 67 |
Podemos comprobar cómo el modelo funciona mejor para clasificar a los que no sobreviven, pero presenta problemas para clasificar correctamente a los que sí sobreviven. La sensibilidad y especificidad del modelo (tomando como punto de corte una probabilidad=0.5) la muestro a continuación (consideramos como grupo de referencia a los que no sobreviven):
## Sensitivity
## 0.4379085
## Specificity
## 0.8989547
La capacidad para clasificar bien a los que no sobreviven (especificidad) es mayor que la capacidad para clasificar correctamente a los que sí sobreviven (sensibilidad).
La curva ROC del modelo se muestra en el siguiente gráfico. Además se añade el mejor punto de corte. Para este modelo sería una prob=0.36, lo que nos da una S=0.53 y una E=0.8.
## Setting levels: control = No, case = Yes
## Setting direction: controls < cases
En conclusión, nuestro modelo de predicción presenta un débil ajuste, sobre todo en la clasificación de los supervivientes. En el grupo test hay un error de un 20% de falsos positivos lo que puede ser excesivo para las intenciones predictivas del modelo. Podríamos incrementar el punto de corte (treshold) a una probabilidad de 0.6. Con ello mejoramos la Especificidad y los falsos positivos (tan solo 2.4%) pero a costa de una reducción de la sensibilidad al 36%, tal y como muestro en el siguinte gráfico:
En el análisis previo he considerado a todos los pasajeros de forma global. Sin embargo ahora me voy a centrar en lo que ocurrió, en base a los datos que tengo, con los niños de viajaron en el Titánic. Los siguientes gráficos de barras nos muestran cómo se repartían en función de la Clase en la que viajaben, en función del sexo y en función de los que sobrevivieron:
Si bien la mortalidad fue en torno al 50% en los niños, ésta se reparte de forma muy desigual en función del sexo y de la clase en la que viajaban. Me llama la atención que entre los niños de primera y segunda clase no hubo fallecidos, pero sí en los de tercera clase (un 66%). De la misma forma la mortalidad entre los niños varones fue muy superior a la de las niñas (54% vs 37%).
Así pues, podemos plantearnos el realizar una regresión logística y comprobar como se comporta el modelo en los niños.
El test de Razón de Verosimilitud, calculado a continuación, es significativo cuando comparo el modelo con el modelo nulo:
| Resid. Df | Resid. Dev | Df | Deviance | Pr(>Chi) |
|---|---|---|---|---|
| 86 | 118.65779 | NA | NA | NA |
| 83 | 80.20553 | 3 | 38.45225 | 0 |
La bondad de ajuste del modelo, mediante el test de Hosmer-Lemeshow, de nuevo vemos que es significativa, lo que nos indica que el modelo no se ajusta muy bien a los datos observados en los niños:
## Warning in Ops.factor(1, y): '-' not meaningful for factors
##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: train$Survived, fitted(model)
## X-squared = 87, df = 7, p-value = 5.551e-16
El porcentaje de variabilidad explicado, medido mediante la R2 de McFadden, es del 32%, superior al modelo global calculado inicialmente:
## fitting null model for pseudo-r2
## McFadden
## 0.3240601
Ahora paso a validar el modelo en el grupo test. A continuación muestro el porcentaje de concordancia en la clasificación de los supervivientes con el modelo. En este caso el porcentaje de aciertos en la clasificación es del 68%:
## [1] 0.6818182
La tabla de confusión la muestro a continuación:
| Asignacion/Sobrevive | No | Yes |
|---|---|---|
| No | 10 | 2 |
| Yes | 5 | 5 |
Y la sensibilidad y especificidad del modelo en el grupo test es:
## Sensitivity
## 0.7142857
## Specificity
## 0.6666667
Vemos que la S mejora la del modelo global, pero la E disminuye. A continuación construyo la curva ROC y calculo el AUC, señalando el punto óptimo de corte:
## Setting levels: control = No, case = Yes
## Setting direction: controls < cases
Vemos que para un punto de corte de 0.77 conseguimos una especificidad del 100%, con una sensibilidad del 71%, y con un AUC del 81%. Tomando este punto de corte, el porcentaje de aciertos es de:
## [1] 0.9090909
Y la tabla de confusión:
| Asignacion/Sobrevive | No | Yes |
|---|---|---|
| No | 15 | 2 |
| Yes | 0 | 5 |
Así pues, el modelo en niños consigue un mejor ajuste que el modelo considerado de forma global.