Dataset: heart failure
Eduardo González Constán
28/6/2021
Para este ejercicio voy a utilizar uno de los numerosos datasets que ofrece la magnífica web de Kagle. Comienzo por uno relativamente sencillo: heart failure prediction. La estructura del dataset es la siguiente:
## 'data.frame': 299 obs. of 13 variables:
## $ age : num 75 55 65 50 65 90 75 60 65 80 ...
## $ anaemia : int 0 0 0 1 1 1 1 1 0 1 ...
## $ creatinine_phosphokinase: int 582 7861 146 111 160 47 246 315 157 123 ...
## $ diabetes : int 0 0 0 0 1 0 0 1 0 0 ...
## $ ejection_fraction : int 20 38 20 20 20 40 15 60 65 35 ...
## $ high_blood_pressure : int 1 0 0 0 0 1 0 0 0 1 ...
## $ platelets : num 265000 263358 162000 210000 327000 ...
## $ serum_creatinine : num 1.9 1.1 1.3 1.9 2.7 2.1 1.2 1.1 1.5 9.4 ...
## $ serum_sodium : int 130 136 129 137 116 132 137 131 138 133 ...
## $ sex : int 1 1 1 1 0 1 1 1 0 1 ...
## $ smoking : int 0 0 1 0 0 1 0 1 0 1 ...
## $ time : int 4 6 7 7 8 8 10 10 10 10 ...
## $ DEATH_EVENT : int 1 1 1 1 1 1 1 1 1 1 ...En primer lugar, tranformo en factores aquellas variables que presentan categorías (anaemia, diabetes, etc).
El tipo de variables queda como sigue:
Vemos como no tenemos datos perdidos. El 46% de las variables son categóricas, y el 54% son variables contínuas.
1. EXPLORACIÓN Y ESTADÍSTICA DESCRIPTIVA DE LA BASE DE DATOS .
Respecto a las variables categóricas, en el siguiente gráfico de barras vemos cómo se distribuyen:
Y la relación de las variables categóricas con la variable respuesta de este estudio (‘Death_event’) se muestra a continuación:
En principio no se aprecian diferencias importantes entre las variables tipo factor y el evento muerte. Aplico un test Chi2, los p-value son los siguientes:
## [1] "anaemia ..... 0.307"
## [1] "diabetes ..... 1"
## [1] "high_blood_pressure ..... 0.214"
## [1] "sex ..... 1"
## [1] "smoking ..... 0.932"A continuación muestro la relación entre las variables contínuas y la variable respuesta:
En este caso se aprecia que sí que existen diferencias importantes en la distribución de las variables cuantitativas en relación a la variable respuesta. El p valor del test t de Student lo muestro en la siguiente tabla:
| variable | p.value |
|---|---|
| age | 0.0000474 |
| creatinine_phosphokinase | 0.3692160 |
| ejection_fraction | 0.0000096 |
| platelets | 0.3993231 |
| serum_creatinine | 0.0000640 |
| serum_sodium | 0.0018723 |
| time | 0.0000000 |
Vemos que todas ellas, excepto la CPK y las plaquetas, presentan diferencias significativas entre las categorías de la variable respuesta.
Por otra parte la correlación entre las diferentes variables contínuas vemos que es muy baja. Así pues podemos descartar problemas de multicolinealidad.
2. MODELO DE REGRESIÓN.
Ahora ya puedo construir un modelo de regresión logística, tomando como variable respuesta la probabilidad de muerte tras un evento cardiaco, y como variables predictoras el resto. Divido el dataset en un grupo de entrenamiento (80% de las observaciones, n=239), donde desarrolaré el modelo, y un grupo de test (20% de las observaciones, n=60), donde lo validaré.
En primer lugar calculo la significación global del modelo (test de la razón de verosimilitud) comparando el modelo completo vs modelo nulo:
| Resid. Df | Resid. Dev | Df | Deviance | Pr(>Chi) |
|---|---|---|---|---|
| 226 | 176.8909 | NA | NA | NA |
| 238 | 303.3212 | -12 | -126.4303 | 0 |
Comprobamos que efectivamente el modelo completo es significativo . A continuación vemos el summary del modelo completo:
##
## Call:
## glm(formula = DEATH_EVENT ~ ., family = binomial, data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3873 -0.5558 -0.2288 0.4220 2.6183
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.604e+01 6.929e+00 2.315 0.020614 *
## age 4.267e-02 1.858e-02 2.296 0.021685 *
## anaemia1 -8.625e-02 4.059e-01 -0.212 0.831739
## creatinine_phosphokinase 4.401e-04 3.039e-04 1.448 0.147565
## diabetes1 1.494e-01 3.894e-01 0.384 0.701179
## ejection_fraction -7.910e-02 1.833e-02 -4.315 1.60e-05 ***
## high_blood_pressure1 -9.075e-02 3.976e-01 -0.228 0.819442
## platelets -1.034e-06 2.510e-06 -0.412 0.680329
## serum_creatinine 7.122e-01 1.965e-01 3.625 0.000289 ***
## serum_sodium -1.108e-01 4.970e-02 -2.230 0.025748 *
## sex1 -3.200e-01 4.747e-01 -0.674 0.500263
## smoking1 1.226e-01 4.732e-01 0.259 0.795582
## time -1.911e-02 3.187e-03 -5.995 2.04e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 303.32 on 238 degrees of freedom
## Residual deviance: 176.89 on 226 degrees of freedom
## AIC: 202.89
##
## Number of Fisher Scoring iterations: 6Seguidamente aplico el método stepwise para seleccionar las mejores variables regresoras. El summary lo muestro a continuación:
##
## Call:
## glm(formula = DEATH_EVENT ~ age + creatinine_phosphokinase +
## ejection_fraction + serum_creatinine + serum_sodium + time,
## family = binomial, data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4293 -0.5359 -0.2322 0.4333 2.6992
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 15.7590144 6.7319069 2.341 0.01924 *
## age 0.0401335 0.0179025 2.242 0.02498 *
## creatinine_phosphokinase 0.0004069 0.0002851 1.427 0.15355
## ejection_fraction -0.0767124 0.0178937 -4.287 1.81e-05 ***
## serum_creatinine 0.7195121 0.1913735 3.760 0.00017 ***
## serum_sodium -0.1115868 0.0483534 -2.308 0.02101 *
## time -0.0189662 0.0030763 -6.165 7.04e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 303.32 on 238 degrees of freedom
## Residual deviance: 177.71 on 232 degrees of freedom
## AIC: 191.71
##
## Number of Fisher Scoring iterations: 5Vemos que la variable creatinine_phosphokinase es no significativa. A continuación vamos a contrastar el modelo con y sin dicha variable:
| Resid. Df | Resid. Dev | Df | Deviance | Pr(>Chi) |
|---|---|---|---|---|
| 233 | 180.4331 | NA | NA | NA |
| 232 | 177.7095 | 1 | 2.723583 | 0.0988758 |
Comprobamos que no existen diferencias significativas entre ambos modelos (p=0.09). La inclusión de esta variable sí que permite una ligera disminución del AIC (192.43 vs 191.70). No obstante, puesto que la disminución del error (medido como AIC) es mínima decido eliminar dicha variable del modelo y quedarme con el modelo reducido. Así pues el summary del modelo reducido es:
##
## Call:
## glm(formula = DEATH_EVENT ~ age + serum_sodium + serum_creatinine +
## ejection_fraction + time, family = binomial, data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3361 -0.5587 -0.2418 0.5083 2.6998
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 15.13649 6.59974 2.293 0.0218 *
## age 0.03648 0.01728 2.111 0.0348 *
## serum_sodium -0.10397 0.04719 -2.203 0.0276 *
## serum_creatinine 0.72641 0.18395 3.949 7.85e-05 ***
## ejection_fraction -0.07696 0.01773 -4.341 1.42e-05 ***
## time -0.01881 0.00305 -6.167 6.98e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 303.32 on 238 degrees of freedom
## Residual deviance: 180.43 on 233 degrees of freedom
## AIC: 192.43
##
## Number of Fisher Scoring iterations: 5El test de Hosmer-Lemeshow es significativo, lo que nos indica que la bondad de ajuste no es tan buena como nos gustaría:
## Warning in Ops.factor(1, y): '-' not meaningful for factors##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: train$DEATH_EVENT, fitted(lr.model.2)
## X-squared = 239, df = 7, p-value < 2.2e-16El porcentaje de incertidumbre explicada por el modelo reducido lo calculo mediante la R2 de McFadden, en este caso comprobamos que es del 40%:
## fitting null model for pseudo-r2## McFadden
## 0.4051419A continuación aplico el modelo al grupo test para su validación. El porcentaje de individuos correctamente clasificados es del 85%:
## [1] 0.85La tabla de confusión la muestro a continuación:
| Clasificacion/Muerte | 0 | 1 |
|---|---|---|
| 0 | 40 | 6 |
| 1 | 3 | 11 |
Los valores de sensibilidad (S), especificidad (E), valor predictivo positivo y valor predictivo negativo son:
## Sensitivity Specificity Pos Pred Value Neg Pred Value
## 0.6470588 0.9302326 0.7857143 0.8695652Vemos cómo la capacidad de clasificar a los que no fallecen (especificidad) es muy superior a la capacidad del modelo de clasificar a los que fallecen (sensibilidad).
Para evaluar todos los posibles puntos de corte que definen la S y E construyo la curva ROC y muestro el mejor punto de corte:
## Setting levels: control = 0, case = 1## Setting direction: controls < cases
Observamos cómo el punto de corte óptimo lo obtenemos para una prob=0.4, con lo que conseguimos una S=82% y una E=88%. Puesto que este estudio trata de predecir las fallecidos tras un fallo cardiaco, podría interesarnos mejorar la S para no fallar al predecir las muertes. Ello implicaría un aumento considerable de los falsos positivos (disminución de la E). En la siguiente curva ROC muestro un supuesto punto de corte que presenta una elevada S pero menor E:
En este caso vemos que tomando como punto de corte una probabilidad de 0.10, conseguimos un 94% de S, pero a costa de un 56% de falsos positivos.