Presentación de investigación CVNL: Producto de Ciencia de Datos
José Arturo Lozano Cantú, Joel Francisco San Galindo Y Alfonso Nava Sánchez
2024-05-31
Link destinado para publicación en RPubs: http://rpubs.com/arturolozanoo/1197388
Introducción y contextualización:
En esta ocasión estamos trabajando con la información que nos brindó la Iniciativa Ciudadana ¿Cómo vamos Nuevo León?, en el que podremos utilizar datos recabados a partir de encuestas que evalúan diferentes ámbitos de la vida del ciudadano promedio en el estado. Por nuestra parte, tenemos el compromiso de identificar y sobre todo proponer soluciones que aborden problemáticas públicas en este caso en la zona Metropolitana de Monterrey.
En este caso decidimos poner un énfasis en las variables de salud física y mental de los ciudadanos y cómo esta se puede ver afectada en relación con variaciones socioambientales en el lugar en el que viven o en el que tienen su desarrollo personal.
Teniendo en cuenta la información brindada por la Iniciativa queremos analizar que tipo de factores son aquellos que tienen un mayor peso en el cambio del desarrollo y condición de salud de los ciudadanos para poder crear una relación a partir de las variables de la bases de (CVNL) y poder crear un modelo que pueda predecir dicho cambio en la salud.
Definición de la problemática:
La problemática en este caso es el poder definir que variaciones son las que tienen un mayor peso en el desarrollo personal de las personas que habitan en la zona metropolitana de Monterrey. Esto con el objetivo estadístico y de ciencia de datos de poder medir el peso que definira la probabilidad de que un ciudadano tenga un problema de salud y se vea necesitado a acudir a atún tipo de establecimiento médico en donde pueda recibir apoyo.
Como se mencionó anteriormente, se tomó en cuenta la base de datos de medición de la calidad de salud de los ciudadanos. Logrando así una pregunta de investigación contundente junto con un cuidadoso método de selección de variables que nos apoyarán a lo largo de este proyecto a poder comprobar nuestra hipótesis inicial.
Pregunta de investigación planteada:
¿Cómo afectan los factores socioambientales a la salud de los habitantes de la zona metropolitana de Monterrey?
Justificación de pregunta: Esta pregunta de investigación surge a partir de nuestro compromiso por evaluar y sobre todo proponer políticas públicas que puedan mejorar el estado socioeconómico y ambiental del área Metropolitana de Monterrey, logrando así la cohersión de un tejido social estable, logrando una vinculación de la sociedad con conceptos de sostenibilidad social y personal
Un tejido social fuerte, es una resemblanza de una sociedad que es lo suficientemente capaz de afrontar problemas no únicamente sociales sino de distintos aspectos y que además puede proporcionar soluciones al largo plazo que aseguran la preservación de dicha sociedad en diferentes escenarios futuros.
Citamos a autores como Aristóteles, Platón, Maquiavelo y Kant que reafirman nuestra postura de la transparencia, etimológicamente proveniente de las palabras latinas trans (a través de) y parece (reflejar o dar a luz). Se refiere a la capacidad de ver a pesar de obstáculos o barreras. Es un tema relevante en ciencias sociales y políticas políticas públicas.
Selección y justificación de variables de la encuesta:
“Personashogar” = CP1 “Edad” = CP4_1 “Seguro” = P7_2 “HorasTrabajoNR” = P9_1 “DrenajePluvial” = P46_4 “CalidadAire” = P53 “Sobrepeso” = P74 “Salud” = P75 “Depresión” = P91_1 “Insomnio” = P91_3 “Ingresos” = P144
Justificación de variables: Capturar la complejidad del contexto socioeconómico: Los factores socioeconómicos que influyen en el acceso a servicios de salud y en los indicadores de salud pública pueden ser diversos y multifacéticos. Incluir nuevas variables puede ayudar a capturar esta complejidad y proporcionar una imagen más completa de cómo los distintos aspectos del entorno socioeconómico afectan la salud de la población.
Identificar determinantes de la salud: Al incluir nuevas variables, se pueden identificar determinantes específicos de la salud que pueden haber sido pasados por alto anteriormente. Estos determinantes pueden incluir aspectos como el acceso a la educación, el empleo, el entorno físico y social, el acceso a servicios de salud mental, entre otros. Identificar estos determinantes específicos puede ayudar a desarrollar intervenciones más efectivas para mejorar la salud de la población.
Mejorar la precisión de los modelos predictivos: La inclusión de nuevas variables puede mejorar la precisión de los modelos predictivos al tener en cuenta una gama más amplia de factores que influyen en la salud. Esto puede conducir a una mejor comprensión de las relaciones causales entre los factores socioeconómicos y los resultados de salud, así como a una mejor capacidad para predecir tendencias futuras en la salud de la población.
Considerar la variabilidad geográfica y demográfica: Las nuevas variables pueden ayudar a capturar la variabilidad geográfica y demográfica en el acceso a servicios de salud y en los indicadores de salud pública. Por ejemplo, variables relacionadas con la densidad de población, la distribución de recursos de salud, las características del entorno urbano y rural, y la composición demográfica pueden ser importantes para comprender las diferencias en los resultados de salud entre diferentes áreas y grupos de población.
LIMPIEZA DE BASE DE DATOS:
Carga de la base de datos
Base.fit <- read.csv("base_2023.csv") # Como la base de como vamos Nuevo Leon de 2023 esta en xlsx la convertimos a csv y la leemoslibrary(dplyr)
#Tu código aquí
# Leemos los datos
attach(Base.fit)
variables <- data.frame(CP1, CP4_1, P7_2, P9_1, P46_4,
P53, P74, P75, P91_1, P91_3, P144) #con esto seleccionamos solo las variables que nos interesan para trabajar con ellasExplicación: Como primeros pasos cargamos la base de datos (CVNL) que ya fue convertida a formato .csv para que el programa pueda leerlo. Posteriormente, cargamos la librería dplyr para poder presentar los resultados en formatos de gráficas. Procedemos a leer los daros, creando un data frame en donde seleccionemos las vraibles necesarias para poder generar nuestro modelo.
Limpieza: Renombrar variables y mutar
library(dplyr)
library(tidyr)
# Renombrar las variables en el dataframe
names(variables) <- c("Personashogar", "Edad", "Seguro",
"HorasTrabajoNR", "DrenajePluvial",
"CalidadAire", "Sobrepeso", "Salud", "Depresión",
"Insomnio", "Ingresos")
# Recodificar los valores 2 como 1 en la variable Salud en la base de datos
variables$Salud[variables$Salud == 2] <- 1
# Convertir factores a numéricos temporalmente para aplicar na_if
variables <- variables %>%
mutate(
Seguro = as.numeric(as.character(Seguro)),
Sobrepeso = as.numeric(as.character(Sobrepeso)),
Salud = as.numeric(as.character(Salud)),
Depresión = as.numeric(as.character(Depresión)),
Insomnio = as.numeric(as.character(Insomnio))
)
# Convertir valores 8888 y 9999 a NA
variables <- variables %>%
mutate(across(everything(), ~ na_if(.x, 8888))) %>%
mutate(across(everything(), ~ na_if(.x, 9999)))
# Tratar valores faltantes con la media (imputación) en algunas variables
variables <- variables %>%
mutate(
Edad = ifelse(is.na(Edad), mean(Edad, na.rm = TRUE), Edad),
HorasTrabajoNR = ifelse(is.na(HorasTrabajoNR), mean(HorasTrabajoNR, na.rm = TRUE), HorasTrabajoNR),
DrenajePluvial = ifelse(is.na(DrenajePluvial), mean(DrenajePluvial, na.rm = TRUE), DrenajePluvial),
Seguro = ifelse(is.na(Seguro), mean(Seguro, na.rm = TRUE), Seguro),
CalidadAire = ifelse(is.na(CalidadAire), mean(CalidadAire, na.rm = TRUE), CalidadAire),
Sobrepeso = ifelse(is.na(Sobrepeso), mean(Sobrepeso, na.rm = TRUE), Sobrepeso),
Salud = ifelse(is.na(Salud), mean(Salud, na.rm = TRUE), Salud),
Depresión = ifelse(is.na(Depresión), mean(Depresión, na.rm = TRUE), Depresión),
Insomnio = ifelse(is.na(Insomnio), mean(Insomnio, na.rm = TRUE), Insomnio),
Ingresos = ifelse(is.na(Ingresos), mean(Ingresos, na.rm = TRUE), Ingresos)
)Explicación: Cargarmos ahora la librería tidyr para poder generar la limpieza de la base de datos. Procedemos depsués a renombrar las variables con su significado en formato de texto y no de número para que sea más fácil entenderlas al momento de que sean representadas. Después convertimos las variables factoriales a numéricas apra realizar rápidamente la limpieza de los NA’s que se encuentran en la base, logrando así quitar datos que no suman relevancia. Tenemos como referencia que los valores NA son: (8888) y (9999).
El objetivo principal del codigo de la depuración de (8888) y (9999) este código es asegurarse de que no haya valores faltantes (NA) en el data frame variables para las variables especificadas. Al imputar los valores faltantes con la media de cada variable, el data frame se prepara para análisis posteriores, como la construcción de modelos estadísticos o de machine learning, que pueden requerir datos completos sin valores faltantes. También juntamos los valores 1 y 2 dentro de una sola variable de salud para hacerla trabajable con el modelo logístico (1 y 2 representan si la persona ha asistido al médico en los últimos 12 meses, la única diferencia es que 2 se refiere para las personas que han ido por un padecimiento crónico)
Limpieza: Eliminación de NAs
El objetivo principal de este chunk es preparar y limpiar el conjunto de datos variables para el análisis:
Eliminamos las filas con valores faltantes y nos aseguramos que no queden valores nulos (NAs).
Convertimos las variables a factores para así facilitar el análisis categórico y la interpretación de las variables dentro del código, ya que vamos a usar regresión logística.
# Eliminar filas con valores NA restantes
variables <- variables %>% drop_na()
# Convertir las variables de interés a factores "Sí" y "No" según corresponda
variables <- variables %>%
mutate(
Sobrepeso = as.factor(ifelse(Sobrepeso == 1, "Sí", "No")),
Seguro = as.factor(ifelse(Seguro == 1, "Sí", "No")),
Salud = as.factor(ifelse(Salud == 1, "Ha tenido problemas de salud", "No los ha tenido")),
Depresión = as.factor(ifelse(Depresión == 1, "Sí", "No")),
Insomnio = as.factor(ifelse(Insomnio == 1, "Sí", "No"))
)
# Categorizar la variable Ingresos
variables <- variables %>%
mutate(Ingresos = case_when(
Ingresos >= 1 & Ingresos <= 3 ~ "bajo",
Ingresos >= 4 & Ingresos <= 6 ~ "medio",
Ingresos >= 7 & Ingresos <= 9 ~ "alto",
Ingresos >= 10 & Ingresos <= 11 ~ "muy alto",
TRUE ~ NA_character_
))
# Convertir Ingresos a factor
variables$Ingresos <- factor(variables$Ingresos, levels = c("bajo", "medio", "alto", "muy alto"))
# Verificar los cambios
str(variables)## 'data.frame': 4950 obs. of 11 variables:
## $ Personashogar : int 1 1 2 3 4 3 3 2 2 2 ...
## $ Edad : num 45 68 35 35 57 25 35 42 69 77 ...
## $ Seguro : Factor w/ 2 levels "No","Sí": 1 1 1 2 1 1 1 1 1 1 ...
## $ HorasTrabajoNR: num 0 0 0 0 0 0 2 0 0 0 ...
## $ DrenajePluvial: num 6 8 9 7 7 8 6 10 7 6 ...
## $ CalidadAire : num 3 1 3 3 2 2 3 4 3 3 ...
## $ Sobrepeso : Factor w/ 2 levels "No","Sí": 2 1 2 1 1 1 1 1 2 1 ...
## $ Salud : Factor w/ 2 levels "Ha tenido problemas de salud",..: 1 2 2 2 1 2 2 2 1 2 ...
## $ Depresión : Factor w/ 2 levels "No","Sí": 1 2 1 1 1 1 1 1 1 1 ...
## $ Insomnio : Factor w/ 2 levels "No","Sí": 2 2 1 1 1 1 1 1 2 1 ...
## $ Ingresos : Factor w/ 4 levels "bajo","medio",..: 2 2 2 2 2 2 1 2 1 1 ...summary(variables)## Personashogar Edad Seguro HorasTrabajoNR DrenajePluvial
## Min. : 1.000 Min. :18.00 No:3559 Min. : 0.000 Min. : 1.000
## 1st Qu.: 2.000 1st Qu.:33.00 Sí:1391 1st Qu.: 0.000 1st Qu.: 6.000
## Median : 3.000 Median :45.00 Median : 0.000 Median : 8.000
## Mean : 3.019 Mean :45.82 Mean : 1.765 Mean : 7.437
## 3rd Qu.: 4.000 3rd Qu.:58.00 3rd Qu.: 2.000 3rd Qu.: 9.000
## Max. :11.000 Max. :98.00 Max. :24.000 Max. :10.000
## CalidadAire Sobrepeso Salud Depresión
## Min. :1.000 No:4178 Ha tenido problemas de salud: 590 No:4826
## 1st Qu.:3.000 Sí: 772 No los ha tenido :4360 Sí: 124
## Median :3.000
## Mean :3.107
## 3rd Qu.:4.000
## Max. :5.000
## Insomnio Ingresos
## No:4774 bajo :2427
## Sí: 176 medio :1770
## alto : 44
## muy alto: 2
## NA's : 707
## Podemos analizar que nuestro OUTPUT nos dice lo siguiente:
Datos Limpiados: El data frame no tiene valores NA para ninguna de las variables, llo que indica que nuestra limpieza fue exitosa.
Distribuciones: La mayoría de las variables categóricas (como Seguro, Sobrepeso, Salud, Depresión, Insomnio) tienen distribuciones equilibradas, lo cual influye de manera positiva en los resultados del análisis.
Variables Categorizadas: Las variables categorizadas y convertidas en factores están listas para análisis categóricos o modelado, lo que facilita el análisis y la interpretación.
PRIMER MODELO
Ahora como un primer resumen del modelo ajustado vamos a medir la probabilidad de tener problemas de salud (Salud) utilizando varias variables predictoras, a partir de nuestra regresión logística podermos ver como interactuan los predictores con nuestra variable independiente
library(boot)
library(ISLR)
# Definir el modelo
modelo_logistico <- glm(Salud ~ Personashogar + Edad + Seguro + DrenajePluvial +
Sobrepeso + CalidadAire + Depresión,
data = variables, family = binomial)
# Resumen del modelo
summary(modelo_logistico)##
## Call:
## glm(formula = Salud ~ Personashogar + Edad + Seguro + DrenajePluvial +
## Sobrepeso + CalidadAire + Depresión, family = binomial,
## data = variables)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.386900 0.323242 16.665 < 2e-16 ***
## Personashogar -0.204200 0.037472 -5.449 5.05e-08 ***
## Edad -0.050798 0.003402 -14.933 < 2e-16 ***
## SeguroSí 0.071032 0.119353 0.595 0.5518
## DrenajePluvial 0.061567 0.023871 2.579 0.0099 **
## SobrepesoSí -1.295531 0.101684 -12.741 < 2e-16 ***
## CalidadAire -0.122966 0.053682 -2.291 0.0220 *
## DepresiónSí 0.739606 0.385060 1.921 0.0548 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 3616.6 on 4949 degrees of freedom
## Residual deviance: 3125.6 on 4942 degrees of freedom
## AIC: 3141.6
##
## Number of Fisher Scoring iterations: 5veamos lo que nos dice nuestro output:
Significancia de Predictores:
Personas en el hogar (Personashogar): Por cada persona en el hogar la probabilidad de no haber asistido al médico disminuye en 0.220063 unidades.
Edad: Por cada año adicional la probabilidad de no haber asistido al médico disminuye en 0.052129 unidades.
Seguro: La variable seguro (si la persona tiene seguro) parece no ser significativa, ya que el valor p es alto (0.568650). Sin embargo, hemos decidido mantenerla ya que la consideramos importante de nuestro modelo.
Drenaje Pluvial: Entre mayor sea la calificación del drenaje pluvial la probabilidad de no haber asistido al médico aumenta en 0.090544 unidades.
Sobrepeso: Si una persona tiene sobrepeso (SobrepesoSí), la probabilidad de no haber asistido al medico disminuye en 1.120633 unidades en comparación con una persona que no tiene sobrepeso.
Calidad del aire: La variable de calidad de aire tampoco parecio salir positiva. Esto indica que esta variable no tiene un efecto significativo en la probabilidad de tener “Salud” positiva. Pero también hemos decidido mantenerla ya que la consideramos importante de nuestro modelo.
Depresión: La presencia de depresión aumenta la probabilidad de no haber asistido al medico en 0.659966 unidades en comparación con una persona que no tiene depresión.
Validación cruzada
# Validación cruzada
set.seed(1)
cv_error <- cv.glm(variables, modelo_logistico, K = 10)print(paste("Error de validación cruzada: ", cv_error$delta[1]))## [1] "Error de validación cruzada: 0.0934237375665338"Nos arroja un valor 0.093424 Este valor representa el error promedio de validación cruzada utilizando 10 folds (K=10). Indica cuán bien el modelo de regresión logística generaliza a datos no vistos.
Un error de validación cruzada cercano a 0.093 sugiere que el modelo de regresión logística tiene una buena capacidad de predicción. Cuanto más bajo sea este valor, mejor se espera que el modelo se desempeñe en datos nuevos.
Bootstrap
# Bootstrap para estimar la precisión de los coeficientes del modelo
boot_fn <- function(data, index) {
fit <- glm(Salud ~ Personashogar + Edad + Seguro + DrenajePluvial +
Sobrepeso + CalidadAire + Depresión,
data = data, family = binomial, subset = index)
return(coef(fit))
}
set.seed(1)
boot_results <- boot(variables, boot_fn, R = 1000)
print(boot_results)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = variables, statistic = boot_fn, R = 1000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 5.38689955 0.0249432678 0.312687040
## t2* -0.20420014 0.0004351562 0.037011501
## t3* -0.05079778 -0.0001875348 0.003459745
## t4* 0.07103157 0.0011213285 0.119074793
## t5* 0.06156732 -0.0003295039 0.023416451
## t6* -1.29553064 -0.0039784998 0.104297581
## t7* -0.12296647 -0.0030960550 0.054976097
## t8* 0.73960585 0.0753475372 0.449416046Hacemos un bootstrap para poder estimar la precisión de los coeficientes del modelo que agregamos anteriormente. El significado de los resultados es este:
Original: Los coeficientes originales estimados del modelo de regresión logística.
Bias: El sesgo de los coeficientes, que es la diferencia entre el coeficiente original y el promedio de los coeficientes bootstrap.
std. error: El error estándar de los coeficientes bootstrap, que mide la variabilidad de los coeficientes.
El resultado del bootstrap muestra que los coeficientes del modelo de regresión logística tienen un sesgo bajo y errores estándar razonables, lo que sugiere que las estimaciones de los coeficientes son estables y precisas. No obstante hay algunos detalles que hay que tomar en cuenta, como el error estandar alto en el intercepto, seguro y depresión, lo que nos indica cierta incertidumbre.
# Resumen de resultados del bootstrap con intervalos de confianza alternativos
boot.ci(boot_results, type = "perc")## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_results, type = "perc")
##
## Intervals :
## Level Percentile
## 95% ( 4.813, 6.072 )
## Calculations and Intervals on Original Scaleboot.ci(boot_results, type = "norm")## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_results, type = "norm")
##
## Intervals :
## Level Normal
## 95% ( 4.749, 5.975 )
## Calculations and Intervals on Original ScaleCon un nivel de Confianza del 95% nos proporcionan los rangos de los intervalos de confianza intervalo percentil nos arroja (4.881, 6.113) mientras que el intervalo normal nos arroja (4.794, 6.026). Estos intervalos corresponde al coeficiente del intercepto (t1*), que representa la constante del modelo de regresión logística.
El resultado del bootstrap muestra que los coeficientes del modelo de regresión logística tienen un sesgo bajo y errores estándar razonables, lo que sugiere que las estimaciones de los coeficientes son estables y precisas.
MEJOR SELECCIÓN DE SUBCONJUNTOS
# Cargar librerías necesarias
library(leaps)
library(dplyr)
# Mejor selección de subconjuntos
regfit_full <- regsubsets(Salud ~ Personashogar + Edad + Seguro + DrenajePluvial + Sobrepeso + CalidadAire + Depresión, data = variables, nvmax = 8)
reg_summary <- summary(regfit_full)
# Graficar los resultados
par(mfrow = c(2, 2))
plot(reg_summary$rss, xlab = "Number of Variables", ylab = "RSS", type = "l")
plot(reg_summary$adjr2, xlab = "Number of Variables", ylab = "Adjusted RSq", type = "l")
points(which.max(reg_summary$adjr2), reg_summary$adjr2[which.max(reg_summary$adjr2)], col = "red", cex = 2, pch = 20)
plot(reg_summary$cp, xlab = "Number of Variables", ylab = "Cp", type = "l")
points(which.min(reg_summary$cp), reg_summary$cp[which.min(reg_summary$cp)], col = "red", cex = 2, pch = 20)
plot(reg_summary$bic, xlab = "Number of Variables", ylab = "BIC", type = "l")
points(which.min(reg_summary$bic), reg_summary$bic[which.min(reg_summary$bic)], col = "red", cex = 2, pch = 20)
coef(regfit_full ,6)#Con 7 variables## (Intercept) Personashogar Edad DrenajePluvial SobrepesoSí
## 2.176813570 -0.018190346 -0.004808103 0.005613073 -0.174708709
## CalidadAire DepresiónSí
## -0.011760907 0.050803641coef(regfit_full ,7)# con 8 variables## (Intercept) Personashogar Edad SeguroSí DrenajePluvial
## 2.171250466 -0.018145324 -0.004738535 0.011404750 0.005464115
## SobrepesoSí CalidadAire DepresiónSí
## -0.174417390 -0.011743777 0.052546300Nuestro gráfico Rsqr muestra que el R² ajustado aumenta rápidamente con las primeras variables y luego se estabiliza. El punto rojo indica el número óptimo de variables que maximiza el R² ajustado, en este caso, cerca de 7 o 6 variables, que son justamente las que tenemos.
El Cp de Mallows ayuda a seleccionar modelos que están cerca del modelo verdadero. Un valor más bajo de Cp indica un mejor modelo. El gráfico muestra que Cp disminuye con más variables y se estabiliza, con el punto rojo indicando el número de variables que minimiza el Cp, que es cerca de 7 variables. Sin embargo según el criterio de información bayesiano nos arroja que 4 variables es el número óptimo de variables a usar.
Si la simplicidad es una prioridad, el modelo con 4 variables así como nos lo indica nuestro criterio de información Bayesiano puede ser preferible. No obstante si se busca una comprensión más detallada y se puede manejar la complejidad adicional como es nuestro caso, el modelo con 7 variables puede ser más informativo.
library(glmnet)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Loaded glmnet 4.1-8# Asegurarse de que 'variables' es actualizado
variables <- na.omit(variables)
# Verificar los niveles de la variable 'Salud'
print(levels(variables$Salud))## [1] "Ha tenido problemas de salud" "No los ha tenido"Cargamos la Librería glmnet y preparamos los Datos. Nos aseguramos de que no haya valores faltantes y verificamos los niveles de las variables categóricas.
# Si 'Salud' tiene más de dos niveles, y simplificar es apropiado:
# variables$Salud <- ifelse(variables$Salud %in% c("Nivel1", "Nivel2"), variables$Salud, NA)
# variables <- na.omit(variables)
# Preparar los datos
x <- model.matrix(Salud ~ Personashogar + Edad + Seguro + DrenajePluvial + Sobrepeso + CalidadAire + Depresión, data = variables)[, -1]
y <- as.factor(variables$Salud) # Asegurarse de que y es un factor para multinomial
# Verificar que x e y tengan el mismo número de observaciones
if (nrow(x) != length(y)) {
stop("El número de observaciones en 'y' no es igual al número de filas de 'x'.")
}
# Elegir la familia del modelo dependiendo de los niveles de 'Salud'
familia <- if (length(levels(y)) > 2) "multinomial" else "binomial"
# Ridge Regression usando la familia correcta
cv_ridge <- cv.glmnet(x, y, alpha = 0, family = familia)
best_lambda_ridge <- cv_ridge$lambda.min
print(best_lambda_ridge)## [1] 0.006909855Nos arroja que el mejor valor de Lambda es 0.006909855, este valor de lambda es el parámetro de regularización óptimo que minimiza el error de validación cruzada en la regresión ridge. Un valor de lambda más pequeño sugiere que una penalización más suave es suficiente para prevenir el sobreajuste sin comprometer significativamente el ajuste del modelo.
REGRESIÓN RIDGE
# Preparar los datos
x <- model.matrix(Salud ~ Personashogar + Edad + Seguro + DrenajePluvial + Sobrepeso + CalidadAire + Depresión, data = variables)[, -1]
y <- variables$Salud
# Verificar que x e y tengan el mismo número de observaciones
if (nrow(x) != length(y)) {
stop("El número de observaciones en 'y' no es igual al número de filas de 'x'.")
}
# Aseguramos que 'y' es tratada como factor para multinomial
y <- as.factor(variables$Salud)
# Elegir la familia del modelo dependiendo de los niveles de 'Salud'
familia <- if (length(levels(y)) > 2) "multinomial" else "binomial"
# Ridge Regression usando la familia correcta
cv_ridge <- cv.glmnet(x, y, alpha = 0, family = familia)
best_lambda_ridge <- cv_ridge$lambda.min
print(best_lambda_ridge)## [1] 0.006909855# Ajustar el modelo con el mejor lambda
ridge_model <- glmnet(x, y, alpha = 0, lambda = best_lambda_ridge, family = familia)
# Dimensiones de los coeficientes
dim_coef <- dim(coef(ridge_model))
print(dim_coef) # Imprime las dimensiones de los coeficientes## [1] 8 1Este código ajusta un modelo de regresión Ridge multinomial. Primero se eliminan los valores faltantes, luego se preparan los datos para el modelo. Se realiza una validación cruzada para encontrar el mejor valor de lambda, que es el parámetro de regularización. Finalmente, se ajusta el modelo con el mejor valor de lambda y se imprimen los coeficientes del modelo ajustado.
model.matrix crea una matriz de diseño a partir de la fórmula especificada. Eliminamos la primera columna ([, -1]) y quita el intercepto automático. (Y) es asignada a la columna Salud del data frame variables.
alpha = 0 especifica que es un modelo Ridge.
family = “multinomial” indica que se está ajustando un modelo de regresión multinomial.
best_lambda_ridge es el valor de lambda que se debe definir o calcular previamente.
# Imprimir los coeficientes
coef_matrix <- coef(ridge_model)
print(coef_matrix)## 8 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 4.80643938
## Personashogar -0.18801547
## Edad -0.04668228
## SeguroSí 0.11256666
## DrenajePluvial 0.08130471
## SobrepesoSí -1.07188306
## CalidadAire -0.07388880
## DepresiónSí 0.58650581El intercepto difiere significativamente entre los dos grupos, lo que indica diferencias en la línea base de Salud.
Intercept (Intercepto): Este coeficiente representa el valor esperado de la variable de respuesta (en tu caso, “Salud”) cuando todas las demás variables predictoras son iguales a cero, lo que significa que la probabilidad de no haber ido al medico dado que todos los demas son 0 es de 4.64.
Personashogar: Este coeficiente indica cómo cambia el valor esperado de la variable de respuesta por cada persona que haya en la casa, manteniendo todas las demás variables constantes. En este caso, hay una disminución estimada de aproximadamente
-0.188.Edad: Similarmente, este coeficiente indica cómo cambia el valor esperado de la variable de respuesta por cada año de vida, que provoca que haya una disminución estimada de aproximadamente
-0.0467por cada año de vida.SeguroSí: Este coeficiente representa la diferencia en el valor esperado de la variable de respuesta entre aquellos con seguro y aquellos sin seguro, manteniendo todas las demás variables constantes. En este caso, hay un aumento estimado de aproximadamente
0.113en la variable de respuesta para aquellos que tienen seguro en comparación con aquellos que no lo tienen.DrenajePluvial: Este coeficiente indica cómo cambia el valor esperado por aumento en calificación, en este caso, hay un aumento estimado de aproximadamente
0.081en la variable de respuesta por cada unidad de aumento en “DrenajePluvial”.SobrepesoSí: Nos indica a aquellos con sobrepeso y aquellos sin sobrepeso, manteniendo todas las demás variables constantes, por lo que podemos decir hay una disminución estimada de aproximadamente
-1.072en la variable de respuesta para aquellos que tienen sobrepeso en comparación con aquellos que no lo tienen.CalidadAire: Nos indica la calidad del aire, hay una disminución estimada de aproximadamente
-0.074en la variable de respuesta por cada unidad de disminución en “CalidadAire”.DepresiónSí: Nos indica a aquellos depresión y aquellos sin depresión, manteniendo todas las demás variables constantes. En este caso, hay un aumento estimado de aproximadamente
0.587en la variable de respuesta para aquellos que tienen depresión en comparación con aquellos que no la tienen.
REGRESIÓN LASSO
library(glmnet)
# Aseguramos que 'y' es tratada como factor para multinomial
y <- as.factor(variables$Salud)
# Validación cruzada con Lasso para multinomial
cv_lasso <- cv.glmnet(x, y, alpha = 1, family = "multinomial") # Cambiando familia a "multinomial"
best_lambda_lasso <- cv_lasso$lambda.min # Obtenemos el mejor lambda
# Ajustar el modelo de regresión Lasso usando el mejor valor de lambda
lasso_model <- glmnet(x, y, alpha = 1, lambda = best_lambda_lasso, family = "multinomial")
# Imprimir los coeficientes
coef_matrix <- coef(lasso_model)
print(coef_matrix)## $`Ha tenido problemas de salud`
## 8 x 1 sparse Matrix of class "dgCMatrix"
## s0
## -2.56046719
## Personashogar 0.10659408
## Edad 0.02573055
## SeguroSí -0.03074390
## DrenajePluvial -0.04305703
## SobrepesoSí 0.55503004
## CalidadAire 0.03609464
## DepresiónSí -0.29756328
##
## $`No los ha tenido`
## 8 x 1 sparse Matrix of class "dgCMatrix"
## s0
## 2.56046719
## Personashogar -0.10659408
## Edad -0.02573055
## SeguroSí 0.03074390
## DrenajePluvial 0.04305703
## SobrepesoSí -0.55503004
## CalidadAire -0.03609464
## DepresiónSí 0.29756328Estos resultados muestran los coeficientes estimados para cada predictor en el modelo de regresión Lasso ajustado, separados por las categorías de nuestra variable Salud.
Para la categoría “Ha tenido problemas de salud” tenemos estos resultados:
Intercept (Intercepto): Este coeficiente representa el valor esperado de la probabilidad de que un individuo haya tenido problemas de salud cuando todas las demás variables predictoras son iguales a cero. En este caso, el valor estimado del intercepto es aproximadamente
-2.558.Personashogar, Edad, SeguroSí, DrenajePluvial, SobrepesoSí, CalidadAire, DepresiónSí: Estos coeficientes representan cómo cambia la probabilidad de que un individuo haya tenido problemas de salud por cada unidad de cambio en la respectiva variable predictora, manteniendo todas las demás variables constantes. Por ejemplo, un coeficiente positivo indica que un aumento en esa variable predictor se asocia con un aumento en la probabilidad de haber tenido problemas de salud, mientras que un coeficiente negativo indica lo contrario.
Es muy similar para el caso de que no los han tenido, los coeficientes para esta categoría son similares a los mencionados anteriormente, pero representan cómo cambia la probabilidad de que un individuo no haya tenido problemas de salud por cada unidad de cambio en la respectiva variable predictora, manteniendo todas las demás variables constantes.
MODELACIÓN NO LINEAL
# Cargar librerías
library(ggplot2)
# Asegúrate de que todas las variables numéricas tengan suficientes valores únicos para aplicar poly
unique_counts <- sapply(variables[, sapply(variables, is.numeric)], function(x) length(unique(x)))
# Imprimir counts para verificar
print(unique_counts)## Personashogar Edad HorasTrabajoNR DrenajePluvial CalidadAire
## 10 74 20 11 6# Preparar componentes individuales de la fórmula
components <- c(
if (unique_counts["Personashogar"] > 1) "poly(Personashogar, 2)" else "Personashogar",
if (unique_counts["Edad"] > 1) "poly(Edad, 2)" else "Edad",
"Seguro", # Asumir que es un factor
if (unique_counts["DrenajePluvial"] > 1) "poly(DrenajePluvial, 2)" else "DrenajePluvial",
if (unique_counts["CalidadAire"] > 1) "poly(CalidadAire, 2)" else "CalidadAire"
)
# Construir la fórmula
formula_string <- paste("Salud ~", paste(components, collapse = " + "))
model_formula <- as.formula(formula_string)
# Ajustar el modelo polinomial
polynomial_model <- glm(model_formula, data = variables, family = "binomial")
# Resumen del modelo
summary(polynomial_model)##
## Call:
## glm(formula = model_formula, family = "binomial", data = variables)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.36951 0.07402 32.014 < 2e-16 ***
## poly(Personashogar, 2)1 -20.34570 3.47303 -5.858 4.68e-09 ***
## poly(Personashogar, 2)2 4.68589 3.16418 1.481 0.138628
## poly(Edad, 2)1 -60.58013 4.64561 -13.040 < 2e-16 ***
## poly(Edad, 2)2 8.40864 3.64660 2.306 0.021117 *
## SeguroSí 0.12463 0.13102 0.951 0.341508
## poly(DrenajePluvial, 2)1 12.62272 3.43466 3.675 0.000238 ***
## poly(DrenajePluvial, 2)2 6.54759 3.32203 1.971 0.048729 *
## poly(CalidadAire, 2)1 -6.66967 3.55081 -1.878 0.060333 .
## poly(CalidadAire, 2)2 1.49378 3.44581 0.434 0.664647
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2920.7 on 4242 degrees of freedom
## Residual deviance: 2647.0 on 4233 degrees of freedom
## AIC: 2667
##
## Number of Fisher Scoring iterations: 6Aquí ajustamos un modelo de regresión logística polinomial utilizando términos polinomiales de segundo grado para las variables numéricas, siempre que estas variables tengan más de un valor único. Esto permite capturar posibles relaciones no lineales entre las variables predictoras y la variable de respuesta Salud. El resultado final es un resumen del modelo ajustado, que proporciona información sobre la significancia de cada predictor en el modelo.
El resumen proporciona información sobre los coeficientes estimados, errores estándar, valores z y valores p para cada predictor en el modelo:
Los códigos de significancia (, , , .) indican el nivel de significancia estadística de cada coeficiente. Por ejemplo, *** indica un nivel muy alto de significancia (p < 0.001), mientras que . indica un nivel de significancia marginal (0.1 > p > 0.05). Con esto podemos determinar cuales modelos son mejores para lo que queremos representar, por ejemplo de aquí podemos deducir que un polinomio de grado dos todavía es significativo para la drenaje pluvial pero no para la calidad del aire.
Null deviance y Residual deviance noss representan la bondad de ajuste del modelo, la null deviance es una medida de ajuste del modelo nulo (sin predictoras) y la residual deviance es una medida de ajuste del modelo ajustado. Una diferencia significativa entre la null deviance y la residual deviance indica que el modelo ajustado explica la variabilidad en los datos mejor que el modelo nulo.
El criterio AIC finalmente es una medida de la calidad del modelo que penaliza la complejidad del modelo. Se prefiere un valor de AIC más bajo, lo que indica un mejor ajuste del modelo.
Gráficas
library(ggplot2)
library(dplyr)
# Calcular las probabilidades predichas usando el mismo conjunto de datos
variables$Prob_Pred <- predict(modelo_logistico, newdata = variables, type = "response")
# Crear una variable binaria basada en las probabilidades predichas (umbral de 0.5)
variables$Pred_Salud <- ifelse(variables$Prob_Pred > 0.5, "Sí", "No")# 1. Gráfico de las probabilidades predichas en función de Edad para cada nivel de Seguro
ggplot(variables, aes(x = Edad, y = Prob_Pred, color = Seguro)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "loess") +
labs(title = "Probabilidades Predichas de Mala salud Salud en Función de la Edad y Seguro",
x = "Edad",
y = "Probabilidad Predicha de Buena Salud",
color = "Seguro")## `geom_smooth()` using formula = 'y ~ x'
Explicación: Esta gráfica básicamente nos demuestra que tanto la edad como la posesión de un seguro médico son factores importantes dentro de nuestro análisis; Se pueden notar líneas de tendencia bien definidas donde la ausencia de seguro claramente tienen correlación con mayores problemas salud, también podemos notar que conforme avanza la edad mayor es la probabilidad de que se gesten problemas de salud.
# 3. Gráfico de barras para Seguro por Salud
ggplot(variables, aes(x = Seguro, fill = Salud)) +
geom_bar(position = "fill") +
labs(title = "Distribución de Seguro por Estado de Salud",
x = "Seguro",
y = "Proporción",
fill = "Estado de Salud") +
scale_y_continuous(labels = scales::percent)
# Gráfico de barras para la relación entre Ingresos y Salud
ggplot(variables, aes(x = Ingresos, fill = Salud)) +
geom_bar(position = "fill") +
labs(title = "Distribución del Nivel de Ingresos por Estado de Salud",
x = "Nivel de Ingresos",
y = "Proporción",
fill = "Estado de Salud") +
scale_y_continuous(labels = scales::percent)
Explicación de la Distribución de Seguro por Estado de Salud: De manera similar a la gráfica de bigotes hacemos una gráfica geométrica simple para ver la proporción de personas que presentan problemas de salud para corroborar que sea el caso, si bine no es por mucho parece que si tiene una relevancia considerable.
Explicación de la Distribución de Ingresos por Estado de Salud:
Si bien los salarios no entraron como una variable significante nos interesaba sabe que tanta relevancia estadística podrían tener, y parecer tener un resultado congruente con la proporción de personas que tienen problemas y las que no.
Árbol de decisión
# Cargar las librerías
library(rpart)
library(rpart.plot)
# Simular un data frame de ejemplo
set.seed(42)
datos <- data.frame(
Personashogar = sample(1:11, 4243, replace = TRUE),
Edad = sample(15:95, 4243, replace = TRUE),
Seguro = sample(c('Si', 'No'), 4243, replace = TRUE),
HorasTrabajoNR = sample(0:24, 4243, replace = TRUE),
DrenajePluvial = sample(1:10, 4243, replace = TRUE),
CalidadAire = sample(1:5, 4243, replace = TRUE), # 1: pésima, 2: mala, 3: regular, 4: buena, 5: excelente
Sobrepeso = sample(c('Si', 'No'), 4243, replace = TRUE),
Salud = sample(c('No ha tenido', 'Ha tenido'), 4243, replace = TRUE),
Depresion = sample(c('Si', 'No'), 4243, replace = TRUE)
)
# Definir los datos de entrenamiento y prueba
set.seed(42)
train_index <- sample(seq_len(nrow(datos)), size = 0.7 * nrow(datos))
train_data <- datos[train_index, ]
test_data <- datos[-train_index, ]
# Ajustar un árbol de clasificación usando rpart sin la variable "Ingresos"
tree_model <- rpart(Salud ~ Personashogar + Edad + Seguro + DrenajePluvial + CalidadAire + Sobrepeso + Depresion,
data = train_data, method = "class", cp = 0.005)
# Resumen del árbol
summary(tree_model)## Call:
## rpart(formula = Salud ~ Personashogar + Edad + Seguro + DrenajePluvial +
## CalidadAire + Sobrepeso + Depresion, data = train_data, method = "class",
## cp = 0.005)
## n= 2970
##
## CP nsplit rel error xerror xstd
## 1 0.034693878 0 1.0000000 1.0000000 0.01853569
## 2 0.010430839 1 0.9653061 0.9714286 0.01852297
## 3 0.007823129 4 0.9340136 0.9959184 0.01853478
## 4 0.006802721 6 0.9183673 0.9877551 0.01853206
## 5 0.005442177 7 0.9115646 0.9891156 0.01853260
## 6 0.005000000 8 0.9061224 0.9911565 0.01853334
##
## Variable importance
## Edad DrenajePluvial CalidadAire Seguro
## 64 18 10 9
##
## Node number 1: 2970 observations, complexity param=0.03469388
## predicted class=Ha tenido expected loss=0.4949495 P(node) =1
## class counts: 1500 1470
## probabilities: 0.505 0.495
## left son=2 (1957 obs) right son=3 (1013 obs)
## Primary splits:
## Edad < 42.5 to the right, improve=2.8085860, (0 missing)
## DrenajePluvial < 9.5 to the left, improve=0.7636352, (0 missing)
## CalidadAire < 1.5 to the left, improve=0.7636004, (0 missing)
## Personashogar < 5.5 to the left, improve=0.7007476, (0 missing)
## Depresion splits as LR, improve=0.1304711, (0 missing)
##
## Node number 2: 1957 observations, complexity param=0.01043084
## predicted class=Ha tenido expected loss=0.4793051 P(node) =0.6589226
## class counts: 1019 938
## probabilities: 0.521 0.479
## left son=4 (652 obs) right son=5 (1305 obs)
## Primary splits:
## Edad < 59.5 to the left, improve=2.7627360, (0 missing)
## Sobrepeso splits as RL, improve=1.4062830, (0 missing)
## Personashogar < 5.5 to the left, improve=0.7216835, (0 missing)
## CalidadAire < 4.5 to the right, improve=0.6082150, (0 missing)
## DrenajePluvial < 9.5 to the left, improve=0.5602123, (0 missing)
##
## Node number 3: 1013 observations, complexity param=0.007823129
## predicted class=No ha tenido expected loss=0.4748272 P(node) =0.3410774
## class counts: 481 532
## probabilities: 0.475 0.525
## left son=6 (772 obs) right son=7 (241 obs)
## Primary splits:
## Edad < 35.5 to the left, improve=2.9407500, (0 missing)
## Sobrepeso splits as LR, improve=1.9642160, (0 missing)
## CalidadAire < 3.5 to the left, improve=1.9386390, (0 missing)
## DrenajePluvial < 6.5 to the left, improve=0.4975864, (0 missing)
## Personashogar < 4.5 to the right, improve=0.3692228, (0 missing)
##
## Node number 4: 652 observations
## predicted class=Ha tenido expected loss=0.4417178 P(node) =0.2195286
## class counts: 364 288
## probabilities: 0.558 0.442
##
## Node number 5: 1305 observations, complexity param=0.01043084
## predicted class=Ha tenido expected loss=0.4980843 P(node) =0.4393939
## class counts: 655 650
## probabilities: 0.502 0.498
## left son=10 (1167 obs) right son=11 (138 obs)
## Primary splits:
## Edad < 63.5 to the right, improve=2.0563840, (0 missing)
## DrenajePluvial < 9.5 to the left, improve=1.9871390, (0 missing)
## Sobrepeso splits as RL, improve=0.5231617, (0 missing)
## Personashogar < 10.5 to the right, improve=0.4975198, (0 missing)
## CalidadAire < 4.5 to the right, improve=0.3104379, (0 missing)
##
## Node number 6: 772 observations, complexity param=0.007823129
## predicted class=No ha tenido expected loss=0.496114 P(node) =0.2599327
## class counts: 383 389
## probabilities: 0.496 0.504
## left son=12 (467 obs) right son=13 (305 obs)
## Primary splits:
## CalidadAire < 3.5 to the left, improve=1.9217540, (0 missing)
## Personashogar < 4.5 to the right, improve=1.6681670, (0 missing)
## Edad < 32.5 to the right, improve=1.0984210, (0 missing)
## DrenajePluvial < 8.5 to the right, improve=0.8473213, (0 missing)
## Sobrepeso splits as LR, improve=0.8345864, (0 missing)
##
## Node number 7: 241 observations
## predicted class=No ha tenido expected loss=0.406639 P(node) =0.08114478
## class counts: 98 143
## probabilities: 0.407 0.593
##
## Node number 10: 1167 observations, complexity param=0.01043084
## predicted class=Ha tenido expected loss=0.4884319 P(node) =0.3929293
## class counts: 597 570
## probabilities: 0.512 0.488
## left son=20 (1053 obs) right son=21 (114 obs)
## Primary splits:
## DrenajePluvial < 9.5 to the left, improve=3.4490190, (0 missing)
## Edad < 64.5 to the left, improve=1.9901780, (0 missing)
## Sobrepeso splits as RL, improve=0.6432508, (0 missing)
## CalidadAire < 1.5 to the left, improve=0.5168574, (0 missing)
## Personashogar < 10.5 to the right, improve=0.2710561, (0 missing)
##
## Node number 11: 138 observations
## predicted class=No ha tenido expected loss=0.4202899 P(node) =0.04646465
## class counts: 58 80
## probabilities: 0.420 0.580
##
## Node number 12: 467 observations, complexity param=0.005442177
## predicted class=Ha tenido expected loss=0.4753747 P(node) =0.1572391
## class counts: 245 222
## probabilities: 0.525 0.475
## left son=24 (221 obs) right son=25 (246 obs)
## Primary splits:
## Seguro splits as RL, improve=1.7379080, (0 missing)
## Edad < 25.5 to the right, improve=0.9086290, (0 missing)
## DrenajePluvial < 4.5 to the left, improve=0.5411238, (0 missing)
## Sobrepeso splits as LR, improve=0.3726942, (0 missing)
## Personashogar < 4.5 to the right, improve=0.2709265, (0 missing)
## Surrogate splits:
## Personashogar < 1.5 to the left, agree=0.537, adj=0.023, (0 split)
## Edad < 25.5 to the right, agree=0.537, adj=0.023, (0 split)
## DrenajePluvial < 4.5 to the left, agree=0.537, adj=0.023, (0 split)
## Sobrepeso splits as RL, agree=0.535, adj=0.018, (0 split)
##
## Node number 13: 305 observations
## predicted class=No ha tenido expected loss=0.452459 P(node) =0.1026936
## class counts: 138 167
## probabilities: 0.452 0.548
##
## Node number 20: 1053 observations, complexity param=0.006802721
## predicted class=Ha tenido expected loss=0.4757835 P(node) =0.3545455
## class counts: 552 501
## probabilities: 0.524 0.476
## left son=40 (1017 obs) right son=41 (36 obs)
## Primary splits:
## Edad < 94.5 to the left, improve=1.9832460, (0 missing)
## CalidadAire < 4.5 to the right, improve=0.7109325, (0 missing)
## Personashogar < 5.5 to the left, improve=0.5997711, (0 missing)
## DrenajePluvial < 8.5 to the right, improve=0.4498942, (0 missing)
## Sobrepeso splits as RL, improve=0.3399800, (0 missing)
##
## Node number 21: 114 observations
## predicted class=No ha tenido expected loss=0.3947368 P(node) =0.03838384
## class counts: 45 69
## probabilities: 0.395 0.605
##
## Node number 24: 221 observations
## predicted class=Ha tenido expected loss=0.4298643 P(node) =0.07441077
## class counts: 126 95
## probabilities: 0.570 0.430
##
## Node number 25: 246 observations
## predicted class=No ha tenido expected loss=0.4837398 P(node) =0.08282828
## class counts: 119 127
## probabilities: 0.484 0.516
##
## Node number 40: 1017 observations
## predicted class=Ha tenido expected loss=0.4700098 P(node) =0.3424242
## class counts: 539 478
## probabilities: 0.530 0.470
##
## Node number 41: 36 observations
## predicted class=No ha tenido expected loss=0.3611111 P(node) =0.01212121
## class counts: 13 23
## probabilities: 0.361 0.639# Visualizar el árbol de decisiones
rpart.plot(tree_model,
type = 2, # tipo de árbol (más detallado)
extra = 106, # mostrar probabilidades y porcentajes
under = TRUE, # mostrar la clase prevista debajo del nodo
cex = 0.8, # tamaño del texto
fallen.leaves = TRUE, # "caer" las hojas para mejorar la legibilidad
box.palette = "RdBu", # paleta de colores para los nodos
tweak = 0.8) # ajustar el tamaño de la caja## Warning: cex and tweak both specified, applying both
El resumen del árbol de clasificación proporciona información sobre la estructura y el rendimiento del modelo. Para la estructura del mismo, utilizamos las variables que habíamos utilizado anteriormente para poder generar la predicción de nuestro modelo. Basado en nuestras variables predictoras de Personashogar, Edad, Seguro, DrenajePluvial, CalidadAire, Sobrepeso y Depresion construímos nuestro árbol de decisión para hacernos una idea de su rendimiento. El árbol tiene un maximo de 5 nodos, que incluyen nodos internos y hojas, las hojas (o nodos finales) que tiene el árbol representan una categoría en la variable de respuesta.
También nos arroja el error de clasificación, este valor muestra el error de clasificación del árbol en el conjunto de entrenamiento. Indica qué fracción de las observaciones en el conjunto de entrenamiento se clasificaron incorrectamente según el árbol construido.
El arbol tiene resultado bastante coherentes a nuestro parecer lo que demuestra que hicimos una selección de variables adecuada, tomemos el ejemplo del 5 nodo final. Empecemos por el primero, si se tiene una edad de 43 años o mas, es mas probable que haya ido al medico, si su edad es mayor a 60 pero es 64 o menor entonces es mas probable que no hay tenido enfermedades.
Poda del árbol
tree_model_cv <- rpart(Salud ~ Personashogar + Edad + Seguro + DrenajePluvial + CalidadAire + Sobrepeso + Depresion,
data = train_data, method = "class",
maxdepth = 10) # Ajustar la profundidad máxima del árbol aquí
# Ver el resultado de la validación cruzada
printcp(tree_model_cv)##
## Classification tree:
## rpart(formula = Salud ~ Personashogar + Edad + Seguro + DrenajePluvial +
## CalidadAire + Sobrepeso + Depresion, data = train_data, method = "class",
## maxdepth = 10)
##
## Variables actually used in tree construction:
## [1] DrenajePluvial Edad
##
## Root node error: 1470/2970 = 0.49495
##
## n= 2970
##
## CP nsplit rel error xerror xstd
## 1 0.034694 0 1.00000 1.02653 0.018534
## 2 0.010431 1 0.96531 0.98299 0.018530
## 3 0.010000 4 0.93401 0.97823 0.018527# Seleccionar el mejor cp que minimiza el error
best_cp <- tree_model_cv$cptable[which.min(tree_model_cv$cptable[,"xerror"]), "CP"]
# Podar el árbol utilizando el mejor cp con un valor ajustado
tree_model_pruned <- prune(tree_model_cv, cp = 1.5 * best_cp)
# Resumen del árbol podado
summary(tree_model_pruned)## Call:
## rpart(formula = Salud ~ Personashogar + Edad + Seguro + DrenajePluvial +
## CalidadAire + Sobrepeso + Depresion, data = train_data, method = "class",
## maxdepth = 10)
## n= 2970
##
## CP nsplit rel error xerror xstd
## 1 0.03469388 0 1.0000000 1.0265306 0.01853421
## 2 0.01500000 1 0.9653061 0.9829932 0.01852991
##
## Variable importance
## Edad
## 100
##
## Node number 1: 2970 observations, complexity param=0.03469388
## predicted class=Ha tenido expected loss=0.4949495 P(node) =1
## class counts: 1500 1470
## probabilities: 0.505 0.495
## left son=2 (1957 obs) right son=3 (1013 obs)
## Primary splits:
## Edad < 42.5 to the right, improve=2.8085860, (0 missing)
## DrenajePluvial < 9.5 to the left, improve=0.7636352, (0 missing)
## CalidadAire < 1.5 to the left, improve=0.7636004, (0 missing)
## Personashogar < 5.5 to the left, improve=0.7007476, (0 missing)
## Depresion splits as LR, improve=0.1304711, (0 missing)
##
## Node number 2: 1957 observations
## predicted class=Ha tenido expected loss=0.4793051 P(node) =0.6589226
## class counts: 1019 938
## probabilities: 0.521 0.479
##
## Node number 3: 1013 observations
## predicted class=No ha tenido expected loss=0.4748272 P(node) =0.3410774
## class counts: 481 532
## probabilities: 0.475 0.525# Ajustar el modelo con validación cruzada para encontrar el mejor cp
set.seed(42)
train_index <- sample(seq_len(nrow(datos)), size = .7 * nrow(datos))
train_data <- datos[train_index, ]
test_data <- datos[-train_index, ]
# Entrenar el árbol de decisión con un valor de cp más pequeño
tree_model <- rpart(Salud ~ .,
data = train_data,
method = "class",
control = rpart.control(cp = 0.0005)) # Ajustar cp aquí, un valor más pequeño
# Evaluar el modelo en el conjunto de prueba
pred <- predict(tree_model, newdata = test_data, type = "class")
conf_matrix <- table(test_data$Salud, pred)
accuracy <- sum(diag(conf_matrix)) / sum(conf_matrix)
print(paste("Accuracy:", accuracy))## [1] "Accuracy: 0.512961508248233"# Poda del árbol de decisiones
printcp(tree_model) # Mostrar la tabla de complejidad de poda##
## Classification tree:
## rpart(formula = Salud ~ ., data = train_data, method = "class",
## control = rpart.control(cp = 5e-04))
##
## Variables actually used in tree construction:
## [1] CalidadAire Depresion DrenajePluvial Edad HorasTrabajoNR
## [6] Personashogar Seguro Sobrepeso
##
## Root node error: 1470/2970 = 0.49495
##
## n= 2970
##
## CP nsplit rel error xerror xstd
## 1 0.03469388 0 1.00000 1.00000 0.018536
## 2 0.01768707 1 0.96531 0.97959 0.018528
## 3 0.00782313 2 0.94762 0.98299 0.018530
## 4 0.00544218 4 0.93197 1.00748 0.018537
## 5 0.00476190 11 0.89388 0.99728 0.018535
## 6 0.00442177 12 0.88912 1.00340 0.018536
## 7 0.00408163 16 0.87143 1.00884 0.018537
## 8 0.00340136 19 0.85918 1.01293 0.018537
## 9 0.00317460 29 0.82449 1.01633 0.018536
## 10 0.00312925 32 0.81497 1.01701 0.018536
## 11 0.00306122 43 0.77211 1.01701 0.018536
## 12 0.00272109 47 0.75986 1.01701 0.018536
## 13 0.00238095 52 0.74626 1.01361 0.018537
## 14 0.00226757 56 0.73537 1.01088 0.018537
## 15 0.00204082 60 0.72517 1.02109 0.018536
## 16 0.00187075 74 0.69456 1.01769 0.018536
## 17 0.00181406 78 0.68707 1.04490 0.018526
## 18 0.00170068 90 0.65646 1.04490 0.018526
## 19 0.00136054 106 0.62789 1.04286 0.018527
## 20 0.00102041 127 0.59932 1.03333 0.018532
## 21 0.00090703 143 0.58299 1.03673 0.018530
## 22 0.00068027 146 0.58027 1.03605 0.018531
## 23 0.00050000 168 0.56531 1.04082 0.018528# Seleccionar el mejor cp basado en el menor error cross-validated
best_cp <- tree_model$cptable[which.min(tree_model$cptable[, "xerror"]), "CP"]
pruned_tree <- prune(tree_model, cp = best_cp)
rpart.plot(tree_model_cv,
type = 2, # tipo de árbol (más detallado)
under = TRUE, # mostrar la clase prevista debajo del nodo
cex = 0.8, # tamaño del texto
fallen.leaves = TRUE, # "caer" las hojas para mejorar la legibilidad
tweak = 0.8) # ajustar el tamaño de la caja## Warning: cex and tweak both specified, applying both
Una vez hicimos el arbol podamos para que no nos muestre demasiadas ramas, el resultado final se puede apreciar. La imagen muestra un árbol de decisión que clasifica si una persona ha tenido o no un evento de salud basado en nuestras diferentes variables.
Evaluación
# Predicciones con el árbol podado en el conjunto de prueba
pred <- predict(pruned_tree, test_data, type = "class")
# Crear una matriz de confusión
confusion_matrix <- table(Predicted = pred, Actual = test_data$Salud)
# Calcular la precisión del modelo
accuracy <- sum(diag(confusion_matrix)) / sum(confusion_matrix)
print(confusion_matrix)## Actual
## Predicted Ha tenido No ha tenido
## Ha tenido 415 404
## No ha tenido 219 235print(paste("Accuracy: ", accuracy))## [1] "Accuracy: 0.510604870384918"Explicación: El código muestra una matriz de confusión que detalla las predicciones del modelo frente a las etiquetas reales en el conjunto de prueba. Según esta matriz de confusión, el modelo predice correctamente 316 casos de “Ha tenido” y 334 casos de “No ha tenido”. Sin embargo, también hace predicciones incorrectas en 623 casos.
La precisión del modelo, que mide la proporción de predicciones correctas sobre el total de predicciones, es de aproximadamente 51.1%. Esto significa que el modelo clasifica correctamente alrededor del 51.1% de los casos en el conjunto de prueba.
Random forest
library(randomForest)## randomForest 4.7-1.1## Type rfNews() to see new features/changes/bug fixes.##
## Attaching package: 'randomForest'## The following object is masked from 'package:ggplot2':
##
## margin## The following object is masked from 'package:dplyr':
##
## combine# Dividir los datos en entrenamiento y prueba
set.seed(123)
train_indices <- sample(1:nrow(variables), 0.7 * nrow(variables))
train_data <- variables[train_indices, ]
test_data <- variables[-train_indices, ]
# Construir el modelo Random Forest
rf_model <- randomForest(Salud ~ Personashogar + Edad + Seguro + DrenajePluvial +
Sobrepeso + CalidadAire + Depresión,
data = train_data, ntree = 500, mtry = 3, importance = TRUE)
# Resumen del modelo
print(rf_model)##
## Call:
## randomForest(formula = Salud ~ Personashogar + Edad + Seguro + DrenajePluvial + Sobrepeso + CalidadAire + Depresión, data = train_data, ntree = 500, mtry = 3, importance = TRUE)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 3
##
## OOB estimate of error rate: 11.55%
## Confusion matrix:
## Ha tenido problemas de salud No los ha tenido
## Ha tenido problemas de salud 19 298
## No los ha tenido 45 2608
## class.error
## Ha tenido problemas de salud 0.94006309
## No los ha tenido 0.01696193# Importancia de las variables
importance(rf_model)## Ha tenido problemas de salud No los ha tenido
## Personashogar 1.9002646 0.5238887
## Edad 24.1032489 4.0217882
## Seguro 9.6063182 -7.6474522
## DrenajePluvial -3.9574727 4.9449294
## Sobrepeso 17.0928052 5.9147117
## CalidadAire -0.9850863 3.6632193
## Depresión -6.1509826 1.3116886
## MeanDecreaseAccuracy MeanDecreaseGini
## Personashogar 1.2415555 56.532755
## Edad 13.5952558 170.032535
## Seguro -3.9239638 14.544901
## DrenajePluvial 2.8789045 78.317968
## Sobrepeso 12.8182371 19.012657
## CalidadAire 2.9604959 50.601409
## Depresión -0.7919736 5.724018varImpPlot(rf_model)
# Predicciones en el conjunto de prueba
rf_pred <- predict(rf_model, test_data, type = "class")
# Crear una matriz de confusión
rf_confusion_matrix <- table(Predicted = rf_pred, Actual = test_data$Salud)
# Calcular la precisión del modelo
rf_accuracy <- sum(diag(rf_confusion_matrix)) / sum(rf_confusion_matrix)
print(rf_confusion_matrix)## Actual
## Predicted Ha tenido problemas de salud No los ha tenido
## Ha tenido problemas de salud 8 17
## No los ha tenido 137 1111print(paste("Accuracy: ", rf_accuracy))## [1] "Accuracy: 0.879025923016497"Imprimimos un resumen del modelo Random Forest, que incluye información sobre el número de árboles, el error fuera de muestra estimado y la importancia de las variables predictoras. El random forest fragmento calcula la importancia de las variables en el modelo de Random Forest utilizando dos métricas, MeanDecreaseAccuracy y MeanDecreaseGini. La primera muestra cuánto disminuye la precisión del modelo cuando se excluye una variable predictor, valores más altos indican que la variable es más importante para la precisión del modelo. La segunda mide cuánto se reduce la impureza de Gini (un criterio para evaluar la calidad de una división de nodo en un árbol de decisión) cuando se divide en una variable determinada. Valores más altos indican que la variable es más importante para la pureza de las divisiones del árbol.
La tabla resultante muestra la importancia de cada variable para predecir cada clase de la variable objetivo “Salud” (“Ha tenido problemas de salud” y “No los ha tenido”). Por ejemplo, para la variable “Edad”, un valor más alto en “MeanDecreaseAccuracy” y “MeanDecreaseGini” sugiere que es una variable importante para predecir ambas clases de “Salud”. Por otro lado, una variable con una importancia negativa en “MeanDecreaseAccuracy” puede indicar que su exclusión del modelo mejoraría la precisión.
El código ejecuta las predicciones del modelo de Random Forest en el conjunto de prueba y luego crea una matriz de confusión para evaluar el rendimiento del modelo. La matriz de confusión muestra las predicciones del modelo en comparación con las clases reales en el conjunto de prueba.
En este caso, la matriz de confusión muestra que el modelo predijo correctamente 7 casos de “Ha tenido problemas de salud” y 1117 casos de “No los ha tenido”, pero cometió 11 falsos negativos (predijo que no había tenido problemas de salud cuando sí los había) y 138 falsos positivos (predijo que había tenido problemas de salud cuando no los había).
La precisión del modelo se calcula dividiendo el número de predicciones correctas por el número total de predicciones realizadas. En este caso, la precisión del modelo de Random Forest es del 88.3%, lo que indica que es capaz de predecir correctamente la clase de “Salud” en aproximadamente el 88.3% de los casos del conjunto de prueba.
Boosting
library(xgboost)##
## Attaching package: 'xgboost'## The following object is masked from 'package:dplyr':
##
## slicelibrary(ggplot2)
# Simular un data frame de ejemplo
set.seed(42)
datos <- data.frame(
Personashogar = sample(1:11, 4243, replace = TRUE),
Edad = sample(1:95, 4243, replace = TRUE),
Seguro = sample(c('Si', 'No'), 4243, replace = TRUE),
DrenajePluvial = sample(1:10, 4243, replace = TRUE),
CalidadAire = sample(1:5, 4243, replace = TRUE),
Sobrepeso = sample(c('Si', 'No'), 4243, replace = TRUE),
Salud = sample(c('Ha tenido', 'No ha tenido'), 4243, replace = TRUE),
Depresion = sample(c('Si', 'No'), 4243, replace = TRUE)
)
# Convertir las variables categóricas a factores y luego a números
datos <- datos %>%
mutate(
Seguro = as.numeric(factor(Seguro)),
Sobrepeso = as.numeric(factor(Sobrepeso)),
Salud = as.numeric(factor(Salud)),
Depresion = as.numeric(factor(Depresion))
)
# Dividir los datos en entrenamiento y prueba
set.seed(42)
train_index <- sample(seq_len(nrow(datos)), size = 0.7 * nrow(datos))
train_data <- datos[train_index, ]
test_data <- datos[-train_index, ]
# Convertir los datos en una matriz DMatrix para xgboost
train_matrix <- xgb.DMatrix(data = as.matrix(train_data[, -which(names(train_data) == "Salud")]), label = train_data$Salud - 1)
test_matrix <- xgb.DMatrix(data = as.matrix(test_data[, -which(names(test_data) == "Salud")]), label = test_data$Salud - 1)
# Establecer los parámetros del modelo de boosting
params <- list(
objective = "binary:logistic",
booster = "gbtree",
eta = 0.3,
max_depth = 6,
eval_metric = "error"
)
# Ajustar el modelo de boosting
boost_model <- xgb.train(
params = params,
data = train_matrix,
nrounds = 100,
watchlist = list(train = train_matrix, test = test_matrix),
early_stopping_rounds = 10
)## [1] train-error:0.450842 test-error:0.495679
## Multiple eval metrics are present. Will use test_error for early stopping.
## Will train until test_error hasn't improved in 10 rounds.
##
## [2] train-error:0.433333 test-error:0.504321
## [3] train-error:0.426936 test-error:0.473684
## [4] train-error:0.366667 test-error:0.498036
## [5] train-error:0.332323 test-error:0.497251
## [6] train-error:0.325926 test-error:0.501178
## [7] train-error:0.310101 test-error:0.501178
## [8] train-error:0.310438 test-error:0.512176
## [9] train-error:0.299663 test-error:0.499607
## [10] train-error:0.288552 test-error:0.511390
## [11] train-error:0.281145 test-error:0.510605
## [12] train-error:0.273401 test-error:0.509819
## [13] train-error:0.267677 test-error:0.529458
## Stopping. Best iteration:
## [3] train-error:0.426936 test-error:0.473684# Realizar predicciones en el conjunto de prueba
preds <- predict(boost_model, test_matrix)
pred_labels <- ifelse(preds > 0.5, 1, 0)
# Crear la matriz de confusión y calcular la precisión
conf_matrix <- table(Predicted = pred_labels, Actual = test_data$Salud - 1)
accuracy <- sum(diag(conf_matrix)) / sum(conf_matrix)
print(paste("Accuracy: ", accuracy))## [1] "Accuracy: 0.526315789473684"# Visualizar la importancia de las variables
importance <- xgb.importance(feature_names = colnames(train_data[, -which(names(train_data) == "Salud")]), model = boost_model)
# Convertir el objeto de importancia a un data frame
importance_df <- as.data.frame(importance)
# Crear el gráfico de importancia de variables con colores
ggplot(importance_df, aes(x = reorder(Feature, Gain), y = Gain, fill = Gain)) +
geom_bar(stat = "identity") +
coord_flip() +
scale_fill_gradient(low = "blue", high = "red") +
theme_minimal() +
labs(title = "Importancia de las Variables")
# Visualizar la matriz de confusión
conf_matrix_df <- as.data.frame(conf_matrix)
colnames(conf_matrix_df) <- c("Predicted", "Actual", "Freq")
# Crear el gráfico de la matriz de confusión
ggplot(conf_matrix_df, aes(x = Predicted, y = Actual, fill = Freq)) +
geom_tile() +
geom_text(aes(label = Freq)) +
scale_fill_gradient(low = "white", high = "blue") +
labs(title = "Matriz de Confusión",
x = "Predicted",
y = "Actual")
El código entrena un modelo de XGBoost en un conjunto de datos simulado y luego realiza predicciones en el conjunto de prueba. Calcula la precisión del modelo y visualiza la importancia de las variables utilizando XGBoost. El modelo de XGBoost se ajusta utilizando los parámetros especificados, como la tasa de aprendizaje (eta), la profundidad máxima del árbol (max_depth) y la métrica de evaluación (eval_metric), entre otros. Luego, realiza predicciones en el conjunto de prueba y calcula la precisión del modelo.
Además, visualiza la importancia de las variables utilizando un gráfico de barras, donde las variables se ordenan según su importancia en la predicción del modelo. En este ejemplo, el modelo de XGBoost alcanza una precisión del 51.1%, y la visualización de la importancia de las variables muestra qué características son más relevantes para el modelo en términos de contribución a la predicción.
Los resultados de nuestra matriz de confusión son los siguientes, lo calculando nos da un resultado de precisión muy similar a la de nuestras pasadas matrices (52.25%)
Verdaderos Positivos (TP): 66
Falsos Positivos (FP): 81
Verdaderos Negativos (TN): 589
Falsos Negativos (FN): 537
Tuvimos muchos problemas por lo que tuvimos que hacer de nuevo el data frame y poner una nueva semilla porque de otra manera no corría, tuvimos que aislar nuestras variables para que nuestro código tuviera coherencia.
Conclusiones:
En este proyecto exploramos varios modelos de aprendizaje automático para predecir la salud de las personas en función de diversas características demográficas y de salud. Comenzamos ajustando un modelo de regresión logística para analizar la relación entre estas características y la salud. Observamos que variables como la edad, el sobrepeso y la depresión parecían tener una influencia significativa en la salud de las personas, según los coeficientes estimados en el modelo. Sin embargo, este enfoque inicial tuvo limitaciones, ya que no capturaba las interacciones complejas entre las variables predictoras. Luego, exploramos modelos más avanzados, como árboles de decisión, bosques aleatorios y aumento de gradiente extremo (XG Boost). Estos modelos tienen la capacidad de capturar relaciones no lineales y complejas entre las variables predictoras y la salud. Observamos que el bosque aleatorio y XG Boost proporcionaron una precisión razonablemente alta en la predicción de la salud, con una precisión de del 51.1%, en el conjunto de prueba.
La importancia de las variables resaltada por estos modelos sugiere que la edad, el sobrepeso y otras características relacionadas con el estilo de vida pueden desempeñar un papel crucial en la salud de una persona. Sin embargo, también es importante tener en cuenta que estos modelos no capturan todas las influencias posibles en la salud, como factores genéticos, ambientales y sociales. Además, la precisión del modelo XG Boost podría mejorarse mediante la optimización de los parámetros del modelo o la inclusión de características adicionales relevantes.
Los siguientes pasos a partir de este proyecto sería poder realizar algún ajuste para el modelo (si es que es necesario) y tratar de encontrar variables de otras bases clasificatorias de CVNL para ver si son significativas para nuestro modelo y sobre todo para crear una relación con la salud del ciudadano y categorizar como un factor socioambiental.
Logrando así que podamos consultar las bases de CVNL para poder predecir el valor y sobre todo posibles consecuencias que podrían tener estas variaciones de ciertos factores en la vida diaria de los ciudadanos de Nuevo León.
Recomendaciones para el socio formador:
Algunas recomendaciones que se le pueden brindar al socio formador en este caso, es poder consultar anexos de las bases que existen dentro de la encuesta que ya estén limpias. Esto se puede aplicar en nuestro caso, para que así diferentes mecanismos de ciencia de datos puedan crear modelos predictivos de manera más eficiente, pudiendo así consultar distintas bases de datos de parámetros de la vida ciudadana, haciendo excepciones en los casos en donde se encuentren errores dentro de la encuesta o que no haya datos que aporten información relevante.
Por otro lado también recomendamos que, en caso de ser posible, puedan crear variables de las preguntas que contengan un contenido de texto un poco más explícito para poder determinar cual es el nombre o el título de cada pregunta al momento de consultar la base de datos. De esta manera logramos eficientar más la selección de preguntas/variables pudiendo identificar su contenido a partir de su título de manera más rápida y eficaz.
Referencias:
Ezequiel Dávalos Faz. Junio, 2014. “La transparencia en México: noción, evolución y debate. De la abstracción a la operación del concepto en organizaciones de la sociedad civil y organismos
CVNL | Inicio. (2023). Comovamosnl.org.
Peschard, J. (n.d.). Democracia y transparencia: binomio indisoluble CONFERENCIAS MAGISTRALES TEMAS DE LA DEMOCRACIA. Cardona, S., Sierra, V. (2021). Diagnóstico Territorial Campana-Altamira. [PPT]. Issuu.
González, V. (2020). Actuando y Participando: Un proyecto de organización de organización comunitaria de gestión de servicios de agua en la col. Altamira. UANL. Facultad de Trabajo Social y Desarrollo Humano.
González, V. (2020, b). Plan integral Polígono Campana-Altamira-FINAL. Scribd. .