Como primer paso, se deberá cargar la base con la cual se va a trabajar, al mismo tiempo de cargar las librerias a utilizar
library(readr)
misdatos <- read_csv("C:\\Users\\DELL\\Desktop\\Curso de especializacion\\bases modulo 3\\heart_failure_clinical_records_dataset.csv")## Rows: 299 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (13): age, anaemia, creatinine_phosphokinase, diabetes, ejection_fractio...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Una vez ya activada la base, corresponde hacerle un análisis exploratorio.
misdatos=as.data.frame(misdatos)
str(misdatos)## 'data.frame': 299 obs. of 13 variables:
## $ age : num 75 55 65 50 65 90 75 60 65 80 ...
## $ anaemia : num 0 0 0 1 1 1 1 1 0 1 ...
## $ creatinine_phosphokinase: num 582 7861 146 111 160 ...
## $ diabetes : num 0 0 0 0 1 0 0 1 0 0 ...
## $ ejection_fraction : num 20 38 20 20 20 40 15 60 65 35 ...
## $ high_blood_pressure : num 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 : num 130 136 129 137 116 132 137 131 138 133 ...
## $ sex : num 1 1 1 1 0 1 1 1 0 1 ...
## $ smoking : num 0 0 1 0 0 1 0 1 0 1 ...
## $ time : num 4 6 7 7 8 8 10 10 10 10 ...
## $ DEATH_EVENT : num 1 1 1 1 1 1 1 1 1 1 ...la base contiene 299 observaciones con 13 variables las cuales están:
age= Edad, variable cuantitativa.
anaemia= si posee anemia, variable categórica (si/no).
creatinine_phosphokinase: Nivel de creatinina Variable numérica.
Diabetes: diabetes, variable categórica (si/no).
ejection_fraction: fracción de eyección, variable numérica.
high_blood_pressure: Hipertensión, categórica (si/no).
platelets: Plaquetas, variable numérica.
serum_creatinine:Suero de creatinina, varibale numérica.
serum_sodium: suero de sodio variable numérica.
sex: sexo categórica (femenino/masculino).
smoking: Fumador, categórica (si/no).
time: Dias de ingreso, variable numérica.
DEATH_EVENT: evento de muerte, Variable categórica(muere/vive).
head(misdatos)## age anaemia creatinine_phosphokinase diabetes ejection_fraction
## 1 75 0 582 0 20
## 2 55 0 7861 0 38
## 3 65 0 146 0 20
## 4 50 1 111 0 20
## 5 65 1 160 1 20
## 6 90 1 47 0 40
## high_blood_pressure platelets serum_creatinine serum_sodium sex smoking time
## 1 1 265000 1.9 130 1 0 4
## 2 0 263358 1.1 136 1 0 6
## 3 0 162000 1.3 129 1 1 7
## 4 0 210000 1.9 137 1 0 7
## 5 0 327000 2.7 116 0 0 8
## 6 1 204000 2.1 132 1 1 8
## DEATH_EVENT
## 1 1
## 2 1
## 3 1
## 4 1
## 5 1
## 6 1Se observan los datos de las primeras seis columnas para visualizar como se encuentran en la base se procede a la activación de una paqueteria y se hacen los siguientes análisis.
set.seed(2021)Número de observaciones que tiene la base
n=nrow(misdatos)
n## [1] 299En este apartado se hace una sección de las bases, las cuales se tomarán para la base de entrenamiento el 70% que corresponde a 209 datos y la base de prueba el 30%, conteniendo 90 datos.
R ofrece la función estándar sample() para tomar una muestra de los conjuntos de datos, en donde se guardara en una nueva varibale, y se le indica que de N datos, tomará una muestra del porcentaje establecido de manera aletatoria.
d_ind=sample(n,n*0.70) #0.7*299 Para construir las bases de manera individual, se crea otra nueva variable y se guarda el paso anterior, extrayendo la muestra de las base original “misdatos” dejando libre el segundo campo para que pase todos los datos de cada observación de la muestra, la primer base tiene un 70% de los datos por tanto se tratará como base de entrenamiento.
records=misdatos[d_ind,] #trainingY para la base de prueba solo construimos el complemento de la base.
d_test=misdatos[-d_ind,] #testVerificando que todas las varibales esten en la base de entrenamiento, se verifican los nombres de las varibles con el comando “name”
names(records)## [1] "age" "anaemia"
## [3] "creatinine_phosphokinase" "diabetes"
## [5] "ejection_fraction" "high_blood_pressure"
## [7] "platelets" "serum_creatinine"
## [9] "serum_sodium" "sex"
## [11] "smoking" "time"
## [13] "DEATH_EVENT"La matriz de correlación ayuda a determinar cómo se relacionan o dependen entre sí dos o más variables. Se muestra en formato de tabla, lo que facilita su lectura y comprensión, este sentido solo se toman las varibles de la base de entreanamiento de edad,creatinina, eyección, plaquetas, suero de creatinina y suero de sodio, debido a que son las variables cuantitativas que tiene la base.
cor(records[c(1,3,5,7,8,9)])## age creatinine_phosphokinase ejection_fraction
## age 1.00000000 -0.018732424 0.02564676
## creatinine_phosphokinase -0.01873242 1.000000000 -0.00801319
## ejection_fraction 0.02564676 -0.008013190 1.00000000
## platelets -0.10971729 -0.002011334 0.04137405
## serum_creatinine 0.13568621 0.054162850 0.05174527
## serum_sodium -0.07268846 0.072730561 0.10135513
## platelets serum_creatinine serum_sodium
## age -0.109717290 0.13568621 -0.07268846
## creatinine_phosphokinase -0.002011334 0.05416285 0.07273056
## ejection_fraction 0.041374045 0.05174527 0.10135513
## platelets 1.000000000 -0.08736598 0.07051430
## serum_creatinine -0.087365981 1.00000000 -0.14122197
## serum_sodium 0.070514305 -0.14122197 1.00000000En este caso se puede ver que los valores estan relacionados y por tanto podemos generar un modelo a partir de ello.
Para la generación del modelo se introducirá la biblioteca Mathematical Acceleration Subsystem (MASS) consta de funciones incorporadas escalares y vectoriales matemáticas ajustadas específicamente para un rendimiento óptimo.
Y utilizaremos la función gml que corresponde a General Lineal Model, en donde almacenaremos la base de datos de entrenamiento y la variable predictora de \(Y\) y agregandole la parte de binomial, para que el modelo sea establecido como un modelo logístico binario múltiple.
library(MASS)
my_model=glm(data = records,DEATH_EVENT~.,
family = binomial())
summary(my_model)##
## Call:
## glm(formula = DEATH_EVENT ~ ., family = binomial(), data = records)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.102e+01 7.542e+00 1.461 0.14410
## age 5.624e-02 1.927e-02 2.918 0.00352 **
## anaemia 4.698e-01 4.160e-01 1.129 0.25874
## creatinine_phosphokinase 4.584e-04 3.483e-04 1.316 0.18807
## diabetes 1.068e-01 4.252e-01 0.251 0.80174
## ejection_fraction -7.459e-02 1.905e-02 -3.916 9.01e-05 ***
## high_blood_pressure -2.212e-01 4.511e-01 -0.490 0.62391
## platelets -1.896e-06 2.294e-06 -0.826 0.40870
## serum_creatinine 5.489e-01 1.926e-01 2.849 0.00438 **
## serum_sodium -7.284e-02 5.314e-02 -1.371 0.17043
## sex -1.338e+00 5.321e-01 -2.514 0.01195 *
## smoking 2.293e-01 5.223e-01 0.439 0.66061
## time -2.250e-02 3.911e-03 -5.754 8.73e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 259.12 on 208 degrees of freedom
## Residual deviance: 153.24 on 196 degrees of freedom
## AIC: 179.24
##
## Number of Fisher Scoring iterations: 6Como podemos observar de todas las varibales, las que son significativas para el modelo son la eyección, la edad y el los dias de ingreso de la persona, con un indicador de AIC de 167.02 que sera comparado posteriromente con otro modelo.
Unos de los problemas mas comunes en una regresión lineal múltiple es el de multicolinealidad , ya que lo que se quieres es que las variables respuestas sean ortogonales y esto se da raramente en casos reales, por lo tanto se aplicará la funcion vif, con la librería \(car\) y \(carData\).
library(carData)
library(car)El VIF para cada término del modelo mide el efecto combinado que tienen las dependencias entre los regresores sobre la varianza de ese término. Si hay uno o más VIF grandes, hay multicolinealidad. Lo que indica que si cualquiera de los VIF es mayor que 5 o 10, es indicio de que los coeficientes asociados de regresión están mal estimados debido a la multicolinealidad, sin emabrgo observamos que en este caso los coeficientes asociados no estan mal estimados ya que a todos los valores estan en el rango de 1 a 1.5
vif(my_model)## age anaemia creatinine_phosphokinase
## 1.129135 1.046886 1.150944
## diabetes ejection_fraction high_blood_pressure
## 1.077959 1.220016 1.076420
## platelets serum_creatinine serum_sodium
## 1.059796 1.116415 1.096600
## sex smoking time
## 1.624757 1.451073 1.269458Como hay variables que no son significativas para el modelo, se modificarán las variables con mayor p valor y se quitarán del modelo dos variables categóricas que seria la presión y si la persona es fumadora, y volvemos a ejetutar un nuevo modelo con gml.
my_model2=glm(data=records,DEATH_EVENT~.-high_blood_pressure-smoking,
family = binomial(link = "logit"))
summary(my_model2)##
## Call:
## glm(formula = DEATH_EVENT ~ . - high_blood_pressure - smoking,
## family = binomial(link = "logit"), data = records)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.043e+01 7.403e+00 1.409 0.15877
## age 5.516e-02 1.918e-02 2.876 0.00403 **
## anaemia 4.404e-01 4.110e-01 1.072 0.28387
## creatinine_phosphokinase 4.491e-04 3.389e-04 1.325 0.18509
## diabetes 8.655e-02 4.237e-01 0.204 0.83815
## ejection_fraction -7.481e-02 1.905e-02 -3.928 8.57e-05 ***
## platelets -1.913e-06 2.283e-06 -0.838 0.40205
## serum_creatinine 5.381e-01 1.888e-01 2.850 0.00437 **
## serum_sodium -6.842e-02 5.201e-02 -1.316 0.18831
## sex -1.197e+00 4.562e-01 -2.624 0.00870 **
## time -2.232e-02 3.887e-03 -5.742 9.35e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 259.12 on 208 degrees of freedom
## Residual deviance: 153.64 on 198 degrees of freedom
## AIC: 175.64
##
## Number of Fisher Scoring iterations: 6Se aprecia la existencia de más variables significativas, lo cual es la edad, la eyección, la suero de creatinina, el sexo y el tiempo de ingreso, también podemos comparar el AIC del el segundo modelo.
Para comparar modelos mediante el AIC, es necesario calcular el AIC de cada modelo. podemos considerar Si un modelo es más de 2 unidades AIC inferior a otro, es significativamente mejor que ese modelo en otras palabras, cuanto menor sea el AIC, mejor será el modelo, porque significa que tiene una mayor probabilidad y una menor complejidad, siguiendo esta teoría analizaremos el AIC del segundo modelo.
my_model2$aic## [1] 175.6393Comparando con el primer modelo, el segundo modelo tiene más robustez por tanto seguiremos trabajando bajo ese modelo y a este se le excluiran otras dos variables que son el análisis de las plaquetas y la diabetes,guardandolo en un nuevo modelo.
my_model3=glm(data=records,DEATH_EVENT~.-high_blood_pressure-smoking-platelets-diabetes,
family = binomial(link = "logit"))
summary(my_model3)##
## Call:
## glm(formula = DEATH_EVENT ~ . - high_blood_pressure - smoking -
## platelets - diabetes, family = binomial(link = "logit"),
## data = records)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 10.1176501 7.2122510 1.403 0.16066
## age 0.0543893 0.0187074 2.907 0.00364 **
## anaemia 0.4637488 0.4101194 1.131 0.25815
## creatinine_phosphokinase 0.0004687 0.0003342 1.403 0.16075
## ejection_fraction -0.0738926 0.0188889 -3.912 9.16e-05 ***
## serum_creatinine 0.5413446 0.1876921 2.884 0.00392 **
## serum_sodium -0.0699922 0.0511068 -1.370 0.17083
## sex -1.1717596 0.4515627 -2.595 0.00946 **
## time -0.0219631 0.0038254 -5.741 9.39e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 259.12 on 208 degrees of freedom
## Residual deviance: 154.39 on 200 degrees of freedom
## AIC: 172.39
##
## Number of Fisher Scoring iterations: 6my_model3$aic## [1] 172.388Como en el análisis anterior siguen siendo significativas las mismas variables y el AIC ha mejorado a partir del primero modelo, para el siguiente modelo se suprimirán las variables de nivel de creatinina y el suero de sodio.
my_model4=glm(data=records,DEATH_EVENT~.-high_blood_pressure-smoking-platelets-diabetes-creatinine_phosphokinase-serum_sodium ,
family = binomial(link = "logit"))
summary(my_model4)##
## Call:
## glm(formula = DEATH_EVENT ~ . - high_blood_pressure - smoking -
## platelets - diabetes - creatinine_phosphokinase - serum_sodium,
## family = binomial(link = "logit"), data = records)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.708705 1.385597 0.511 0.60902
## age 0.053327 0.017994 2.964 0.00304 **
## anaemia 0.401068 0.403233 0.995 0.31992
## ejection_fraction -0.073797 0.018512 -3.986 6.71e-05 ***
## serum_creatinine 0.569549 0.176790 3.222 0.00127 **
## sex -1.077485 0.441182 -2.442 0.01460 *
## time -0.021078 0.003659 -5.761 8.37e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 259.12 on 208 degrees of freedom
## Residual deviance: 158.21 on 202 degrees of freedom
## AIC: 172.21
##
## Number of Fisher Scoring iterations: 5EL AIC sigue mejorando y Siguen siendo signifativas las mismas variables, por tanto también se excluirá la variable anemia y también la variable sexo, y se ejecutará el nuevo modelo
my_model5=glm(data=records,DEATH_EVENT~.-high_blood_pressure-smoking-platelets-diabetes-creatinine_phosphokinase-serum_sodium-anaemia-sex-1 ,family = binomial(link = "logit"))
summary(my_model5)##
## Call:
## glm(formula = DEATH_EVENT ~ . - high_blood_pressure - smoking -
## platelets - diabetes - creatinine_phosphokinase - serum_sodium -
## anaemia - sex - 1, family = binomial(link = "logit"), data = records)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## age 0.047627 0.010969 4.342 1.41e-05 ***
## ejection_fraction -0.062928 0.015835 -3.974 7.07e-05 ***
## serum_creatinine 0.612431 0.178637 3.428 0.000607 ***
## time -0.020307 0.003297 -6.159 7.31e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 289.74 on 209 degrees of freedom
## Residual deviance: 165.67 on 205 degrees of freedom
## AIC: 173.67
##
## Number of Fisher Scoring iterations: 5Observamos que el AIc del modelo 5 es el más bajo de todos los creados, y que las variables que están incluidas son las variables numéricas, por tanto se procede a realizar una bondad de ajuste.
Perimero se importan las librerías \(zoo\) que es un estructura para series temporales regulares e irregulares (observaciones ordenadas de Z) y la librería \(lmtest\), que es una librería de prueba de modelos de regresión lineal.
library(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(lmtest)La prueba de Wald es una prueba de hipótesis estadística que se utiliza para evaluar si los parámetros de un modelo estadístico son significativamente diferentes de los valores hipotéticos. Se utiliza ampliamente en el contexto del análisis de regresión, donde se prueba la significancia de coeficientes individuales o grupos de coeficientes en un modelo de regresión.
waldtest(my_model,my_model5)## Wald test
##
## Model 1: DEATH_EVENT ~ age + anaemia + creatinine_phosphokinase + diabetes +
## ejection_fraction + high_blood_pressure + platelets + serum_creatinine +
## serum_sodium + sex + smoking + time
## Model 2: DEATH_EVENT ~ (age + anaemia + creatinine_phosphokinase + diabetes +
## ejection_fraction + high_blood_pressure + platelets + serum_creatinine +
## serum_sodium + sex + smoking + time) - high_blood_pressure -
## smoking - platelets - diabetes - creatinine_phosphokinase -
## serum_sodium - anaemia - sex - 1
## Res.Df Df F Pr(>F)
## 1 196
## 2 205 -9 1.1655 0.3191Como se puede observar el valor de la probabilidad es 0.3191 que es mayor a 0.05,podemos eliminar estos términos del modelo ya que no mejoran significativamente desde el punto de vista estadístico el ajuste general del modelo.
Esto lo reconfirmamos con la siguente prueba de Hosmer-Lemeshow es una prueba estadística para medir la bondad de ajuste y calibración de los modelos de regresión logística . Se utiliza con frecuencia en los modelos de predicción de riesgos .
Para este caso se Calcula la prueba de bondad de ajuste de Hosmer-Lemeshow para un modelo lineal generalizado ajustado a respuestas binarias.
library(ResourceSelection)## Warning: package 'ResourceSelection' was built under R version 4.3.3## ResourceSelection 0.3-6 2023-06-27hoslem.test(records$DEATH_EVENT,fitted(my_model5))##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: records$DEATH_EVENT, fitted(my_model5)
## X-squared = 8.2596, df = 8, p-value = 0.4085La hipótesis nula del test de Hosmer-Lemeshow es que no hay diferencias entre los valores observados y los valores pronosticados (el rechazo este test indicaría que el modelo no está bien ajustado).Como el P valor es mayor que 0.05, nos lleva a no rechazar la hipótesis nula entonces podemos decir que el modelo se ajusta bien a los datos.
A partir de ello podemos sacar los coeficientes de la regresion lineal \(a_1,a_2,a_3,a_4\) para escribir el modelo en su forma general. A continuación, tenemos varias formas de extraerlos.
coef(my_model5)## age ejection_fraction serum_creatinine time
## 0.04762671 -0.06292788 0.61243102 -0.02030746my_model5$coefficients## age ejection_fraction serum_creatinine time
## 0.04762671 -0.06292788 0.61243102 -0.02030746Para saber si nuestros coeficientes son significativos para el modelo, procedemos a sacar los intervalos de confianza
confint(my_model5)## Waiting for profiling to be done...## 2.5 % 97.5 %
## age 0.02712785 0.07051571
## ejection_fraction -0.09572358 -0.03334280
## serum_creatinine 0.26755137 0.98855160
## time -0.02731683 -0.01431199Y como se puede ver en los coeficientes que fueron seleccionados para el modelo en ningún intervalo esta incluido el cero, por tanto la variable aporta información al modelo.
Para probar que el modelo tiene una eficiencia, predeciremos unas observaciones con la función predict()
y=predict(my_model5,type='response')
head(y)## 166 231 70 192 251 102
## 0.39533733 0.10040396 0.78290708 0.03414446 0.02802915 0.42462119Luego, convertiremos a la predicción de si vive o muere
pred=ifelse(as.double(y)>0.5,1,0)
pred## [1] 0 0 1 0 0 0 0 1 1 0 0 0 1 1 0 1 0 0 1 0 1 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 1
## [38] 1 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1 1 1 1 1 1 1 0 1 1 0 1 0 0 0 0 1 0 0 0 0 1
## [75] 1 0 1 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0
## [112] 0 1 0 0 0 0 0 0 1 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 1 1 1 0 0 0
## [149] 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 0 0 0 0 1 1 0 0 0 0 0 1 0 1 1 0 0
## [186] 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 1 1 0 0Se carga la libreria MLmetrisc que es el paquete de Métricas de evaluación del aprendizaje automático Una colección de métricas de evaluación, que incluyen funciones de pérdida, puntuación y utilidad, que miden el rendimiento de la regresión, la clasificación y el ranking
library(MLmetrics)## Warning: package 'MLmetrics' was built under R version 4.3.3##
## Attaching package: 'MLmetrics'## The following object is masked from 'package:base':
##
## RecallUtilizamos la función accurancy en donde nos dara el porcentaje de aciertos del modelo, coloamos las predicciones hechas con el modelo 5 y las comparamos con la base de entrenamiento y se tiene una presicion de 82.77%
Accuracy(pred,records$DEATH_EVENT)## [1] 0.8277512Hacemos una segunda predicción con el primer modelo obtenido y se codifican como si vive o mueren
pred2=ifelse(my_model$fitted.values>0.5,1,0)
head(pred2)## 166 231 70 192 251 102
## 1 0 1 0 0 0Realizamos nuevamente la fucnion acurracy, en donde colocamos las prediciciones del primer modelo y la base de estrenamiento con el evento muerte y obtenemos un acierto de 85.64%
Accuracy(pred2,records$DEATH_EVENT)## [1] 0.8564593