El conjunto de datos de la Encuesta Nacional de Factores de Riesgo (ENFR) 2018 del repositorio del INDEC (disponible en: https://www.indec.gob.ar/indec/web/Institucional-Indec-BasesDeDatos-2) contiene información sobre los comportamientos y condiciones de vida que pueden afectar la salud de la población de 18 años y más de Argentina, así como algunas características de las viviendas y de los hogares donde residen.
El objetivo de este trabajo es predecir la práctica de mamografía (variable control_mamografía) en mujeres de 50 o más años de edad (población objetivo del control mamográfico), a partir de las características sociodemográficas (edad, estado civil, características del hogar/familia, ingresos, educación, ocupación, cobertura de salud, lugar de residencia), de salud (salud percibida, otras prácticas de cuidado o control de la salud, capacidades funcionales en la vida diaria, etc.) y variables relacionadas a estilo de vida (hábitos relacionados a actividad física y alimentación), entre otras.
En este marco, las variables de dicho dataset incluidas en este trabajo son:
Variable de estratificación:
Variables predictoras:
bhch03: sexo del respondente 1.Varón 2.Mujer
bhch04: edad del respondente (años)
bhch05: Situación conyugal 1 Unido/a 2 Casado/a 3 Separado/a 4 Divorciado/a 5 Viudo/a 6 Soltero/a
nivel_instruccion: Nivel de instrucción
cobertura_salud: Cobertura de salud 1 Con obra social, prepaga o servicio de emergencia médica 2 Sólo cobertura pública
condicion_actividad: Condición de actividad 1 Ocupado 2 Desocupado 3 Inactivo
bisg01: En general, ¿usted diría que su salud es… 1 …excelente? 2 …muy buena? 3 …buena? 4 …regular? 5 …mala?
bisg02: En relación a su movilidad, ¿en el día de hoy… 1 …no tiene problemas para caminar? 2 …tiene algunos problemas para caminar? 3 …no puede caminar?
bisg06: En relación a la ansiedad/depresión, ¿en el día de hoy… 1 …no está ansioso ni deprimido? 2 …está moderadamente ansioso o deprimido? 3 …está muy ansioso o deprimido?
control_pap: Realización de Papanicolaou en los últimos 2 años
region: Código de región INDEC
localidades_150: Tamaño de localidad - 150.000 habitantes o más 1. Localidad de 150.000 o más habitantes 0. Localidad de menos de 150.000 habitantes
cant_componentes: Cantidad de miembros del hogar (Cantidad en cifras)
tipo_hogar: Tipo de hogar 1 Hogar unipersonal 2 Hogar multipersonal conyugal completo sin hijos ni otros miembros 3 Hogar multipersonal conyugal completo sin hijos y con otros miembros 4 Hogar multipersonal conyugal completo con hijos sin otros miembros 5 Hogar multipersonal conyugal completo con hijos y con otros miembros 6 Hogar multipersonal conyugal incompleto sin otros miembros 7 Hogar multipersonal conyugal incompleto con otros miembros 8 Hogar multipersonal no conyugal
quintil_uc: Quintil de hogares según ingreso por unidad consumidora 1. Quintil 1 2. Quintil 2 3. Quintil 3 4. Quintil 4 5. Quintil 5
nivel_actividad_fisica: Nivel de actividad física 1. Alto 2. Medio 3. Bajo 99. Ns/Nc
bial08: Pensando en su alimentación o dieta habitual o de todos los días, ¿considera que es… 1…muy saludable? 2…bastante saludable? 3…poco saludable? 4…nada saludable? 99. Ns/Nc
Llamo a las librerías necesarias:
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 4.0.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(tidymodels)
## ── Attaching packages ────────────────────────────────────── tidymodels 1.4.1 ──
## ✔ broom 1.0.10 ✔ rsample 1.3.1
## ✔ dials 1.4.2 ✔ tailor 0.1.0
## ✔ infer 1.0.9 ✔ tune 2.0.1
## ✔ modeldata 1.5.1 ✔ workflows 1.3.0
## ✔ parsnip 1.3.3 ✔ workflowsets 1.1.1
## ✔ recipes 1.3.1 ✔ yardstick 1.3.2
## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ scales::discard() masks purrr::discard()
## ✖ dplyr::filter() masks stats::filter()
## ✖ recipes::fixed() masks stringr::fixed()
## ✖ dplyr::lag() masks stats::lag()
## ✖ yardstick::spec() masks readr::spec()
## ✖ recipes::step() masks stats::step()
library(haven)
Cargo mi dataset, que lo tengo disponible en formato .dta:
data <- read_dta("data/ENFR2018.dta")
Verifico los nombres de las columnas para podes escribirlas correctamente:
names(data)
## [1] "id" "cod_provincia"
## [3] "region" "tamanio_aglomerado"
## [5] "aglomerado" "localidades_150"
## [7] "submuestra" "bhcv01"
## [9] "bhcv02" "bhcv03"
## [11] "bhcv04" "bhcv05"
## [13] "bhcv06" "bhcv07"
## [15] "bhcv08" "bhcv09"
## [17] "bhcv10" "bhcv11"
## [19] "bhho01" "bhho02"
## [21] "bhho03" "cant_componentes"
## [23] "miembros_18" "tipo_hogar"
## [25] "bhih01" "bhih01_02"
## [27] "rango_ingreso" "quintil_uc"
## [29] "imputado" "bhih03"
## [31] "bhch03_j" "bhch04_j"
## [33] "rango_edad_j" "bhch05_j"
## [35] "nivel_instruccion_j" "nivel_instruccion_agrupado_j"
## [37] "bhch10_01_j" "bhch10_02_j"
## [39] "bhch10_03_j" "bhch10_04_j"
## [41] "bhch10_05_j" "bhch10_06_j"
## [43] "bhch10_99_j" "cobertura_salud_j"
## [45] "bhsl01" "bhsl02"
## [47] "bhsl03" "bhsl04"
## [49] "bhsl05" "bhsl06"
## [51] "condicion_actividad_j" "bhch02"
## [53] "bhch03" "bhch04"
## [55] "rango_edad" "bhch05"
## [57] "nivel_instruccion" "nivel_instruccion_agrupado"
## [59] "bhch10_01" "bhch10_02"
## [61] "bhch10_03" "bhch10_04"
## [63] "bhch10_05" "bhch10_06"
## [65] "bhch10_99" "cobertura_salud"
## [67] "bisl01" "bisl02"
## [69] "bisl03" "bisl04"
## [71] "bisl05" "bisl06"
## [73] "condicion_actividad" "bisg01"
## [75] "bisg02" "bisg03"
## [77] "bisg04" "bisg05"
## [79] "bisg06" "biaf01"
## [81] "biaf02_m" "biaf02_99"
## [83] "biaf03" "biaf04_m"
## [85] "biaf04_99" "biaf05"
## [87] "biaf06_m" "biaf06_99"
## [89] "biaf07_m" "biaf07_99"
## [91] "biaf08" "biaf09"
## [93] "biaf10_01" "biaf10_02"
## [95] "biaf10_03" "biaf10_04"
## [97] "nivel_actividad_fisica" "barreras_actividad_fisica"
## [99] "bita01" "bita02"
## [101] "bita02_99" "bita03"
## [103] "bita04" "bita04_01"
## [105] "bita04_02" "bita05"
## [107] "bita06_a" "bita06_b"
## [109] "bita06_b_99" "bita06_c"
## [111] "bita06_d" "bita07"
## [113] "bita07_99" "bita08"
## [115] "bita09_01" "bita09_02"
## [117] "bita09_03" "bita09_04"
## [119] "bita09_05" "bita09_06"
## [121] "bita10_01" "bita10_02"
## [123] "bita10_03" "bita10_04"
## [125] "bita10_05" "bita10_06"
## [127] "bita11" "bita12"
## [129] "bita13" "bita14"
## [131] "bita15" "bita16"
## [133] "consumo_tabaco_100" "ta_paquete_y_armado"
## [135] "ta_dejar_fumar" "ta_otros_productos"
## [137] "hta_nofumadores" "ta_perc_publicidad"
## [139] "ta_percepcion_riesgo" "imagenes_tabaco"
## [141] "biha01" "biha02"
## [143] "biha03" "biha04"
## [145] "biha05_01" "biha05_02"
## [147] "biha06" "biha06_99"
## [149] "biha07" "biha08"
## [151] "biha09" "biha10"
## [153] "biha11" "biha11_99"
## [155] "biha12" "biha13"
## [157] "biha14" "biha15"
## [159] "control_hipertension" "prevalencia_hipertension"
## [161] "bipc01" "bipc02"
## [163] "bipc03" "bipc04"
## [165] "bipc04_99" "bipc05"
## [167] "bipc05_99" "imc"
## [169] "imc_categorias" "bial01"
## [171] "bial02" "bial03"
## [173] "bial03_99" "bial04"
## [175] "bial04_99" "bial05"
## [177] "bial05_99" "bial06"
## [179] "bial06_99" "bial07"
## [181] "bial08" "bial09"
## [183] "bial10" "promedio_fv_diario"
## [185] "consumo_fv" "barreras_fyv"
## [187] "tipo_dieta_razones" "bico01"
## [189] "bico02" "bico03"
## [191] "bico04" "bico05_01"
## [193] "bico05_02" "control_colesterol"
## [195] "prevalencia_colesterol" "bica01"
## [197] "bica02" "bica03_01"
## [199] "bica03_02" "bica03_99"
## [201] "bica04_01_b" "bica04_01_c"
## [203] "bica04_02_b" "bica04_02_c"
## [205] "bica04_03_b" "bica04_03_c"
## [207] "bica04_04" "bica05_01_b"
## [209] "bica05_01_c" "bica05_02_b"
## [211] "bica05_02_c" "bica05_03_b"
## [213] "bica05_03_c" "bica05_04"
## [215] "bica06" "bica07"
## [217] "consumo_regular_riesgo" "consumo_episodico_excesivo"
## [219] "bidi01" "bidi02"
## [221] "bidi03" "bidi04_01"
## [223] "bidi04_02" "bidi04_03"
## [225] "bidi05" "bidi06_01"
## [227] "bidi06_02" "bidi07"
## [229] "bidi08" "bidi09"
## [231] "bidi10" "bidi11"
## [233] "bidi12" "bidi13"
## [235] "bidi14" "control_diabetes"
## [237] "prevalencia_diabetes" "bile01"
## [239] "bile02" "bile03"
## [241] "bipp01" "bipp02"
## [243] "bipp03" "bipp04"
## [245] "control_mamografia" "control_pap"
## [247] "bicc01_01" "bicc01_02"
## [249] "bicc01_03" "bicc02"
## [251] "bicc03" "control_colon"
## [253] "bima01" "bima02"
## [255] "bima03" "bima04_01_a"
## [257] "bima04_01_b" "bima04_02_a"
## [259] "bima04_02_b" "bima04_03_a"
## [261] "bima04_03_b" "promedio_sistolica"
## [263] "promedio_diastolica" "ta_elevada"
## [265] "prevalencia_hipertension_combina" "bima06"
## [267] "bima07" "bima09"
## [269] "bima10" "imc_bima"
## [271] "imc_categorias_bima" "bima12"
## [273] "bima13" "bima14"
## [275] "bimq01" "bimq05"
## [277] "bimq05_01" "glucemia_elevada"
## [279] "prevalencia_glucemia_elevada_com" "findrisc"
## [281] "bimq06" "bimq06_01"
## [283] "colesterol_elevado" "prevalencia_colesterol_combinada"
## [285] "wf1p" "wf2p"
## [287] "wf3p"
Selecciono las columnas con las que voy a trabajar y realizo un filto para trabajar sólo con los datos corresponondientes a mujeres de 50 o más años de edad, eliminando además del dataset las observaciones con datos faltantes de la variable control_mamografía (código 99 se usa como dato faltante en las ENFR):
data_select <- data %>%
filter(bhch03 == 2, bhch04 >= 50, control_mamografia != "99") %>%
select(control_mamografia,
bhch04, bhch05, nivel_instruccion,
cobertura_salud, condicion_actividad,
bisg01, bisg02, bisg06,
control_pap,
region, localidades_150, cant_componentes, tipo_hogar, quintil_uc,
nivel_actividad_fisica, bial08) %>%
drop_na(control_mamografia)%>%
droplevels()
Me aseguro que la variable de estratificación (control_mamografia) y demás variables categóricas (salvo la edad), estén definidas como un factor:
data_select <- data_select %>%
mutate(
across(-bhch04, as.factor)
)
Con la libreria rsample vamos a dividir los datos en 80/20, utilizando control_mamografia como variable de estratificación;
set.seed(123) # Para reproducibilidad
split <- initial_split(data_select, prop = 0.8, strata = control_mamografia)
train_data <- training(split)
test_data <- testing(split)
Controlo los datos de entrenamiento:
train_data
## # A tibble: 5,597 × 17
## control_mamografia bhch04 bhch05 nivel_instruccion cobertura_salud
## <fct> <dbl> <fct> <fct> <fct>
## 1 1 65 2 7 1
## 2 1 84 4 7 1
## 3 1 63 6 6 1
## 4 1 59 2 5 1
## 5 1 60 4 6 1
## 6 1 58 4 7 1
## 7 1 50 6 7 1
## 8 1 70 2 7 1
## 9 1 77 2 7 1
## 10 1 56 3 5 1
## # ℹ 5,587 more rows
## # ℹ 12 more variables: condicion_actividad <fct>, bisg01 <fct>, bisg02 <fct>,
## # bisg06 <fct>, control_pap <fct>, region <fct>, localidades_150 <fct>,
## # cant_componentes <fct>, tipo_hogar <fct>, quintil_uc <fct>,
## # nivel_actividad_fisica <fct>, bial08 <fct>
Controlo los datos de test:
test_data
## # A tibble: 1,400 × 17
## control_mamografia bhch04 bhch05 nivel_instruccion cobertura_salud
## <fct> <dbl> <fct> <fct> <fct>
## 1 1 56 4 7 1
## 2 1 55 5 7 1
## 3 2 80 5 3 1
## 4 2 63 2 5 2
## 5 1 58 3 7 1
## 6 1 60 2 7 1
## 7 2 84 5 3 1
## 8 2 70 5 7 1
## 9 1 67 2 7 1
## 10 1 54 3 5 1
## # ℹ 1,390 more rows
## # ℹ 12 more variables: condicion_actividad <fct>, bisg01 <fct>, bisg02 <fct>,
## # bisg06 <fct>, control_pap <fct>, region <fct>, localidades_150 <fct>,
## # cant_componentes <fct>, tipo_hogar <fct>, quintil_uc <fct>,
## # nivel_actividad_fisica <fct>, bial08 <fct>
Configuro la validación cruzada de 5 folds sobre el conjunto de entrenamiento, manteniendo la proporción de clases (de contol de mamografía) en cada fold:
set.seed(123)
folds <- vfold_cv(train_data, v = 5, strata = control_mamografia)
Controlo:
folds$splits
## [[1]]
## <Analysis/Assess/Total>
## <4476/1121/5597>
##
## [[2]]
## <Analysis/Assess/Total>
## <4478/1119/5597>
##
## [[3]]
## <Analysis/Assess/Total>
## <4478/1119/5597>
##
## [[4]]
## <Analysis/Assess/Total>
## <4478/1119/5597>
##
## [[5]]
## <Analysis/Assess/Total>
## <4478/1119/5597>
Controlo si hay desbalance en mis datos de entrenamiento:
table(train_data$control_mamografia)
##
## 1 2
## 3146 2451
prop.table(table(train_data$control_mamografia))
##
## 1 2
## 0.5620868 0.4379132
Se observa que no existe un desbalance importante en mis datos (56/44), en relación a la variable de estratificación.
Realizo la receta de preprocesamiento, definiendo el set de variables predictoras:
rf_recipe <- recipe(control_mamografia ~
bhch04 + bhch05 + nivel_instruccion + cobertura_salud + # variables demográficas
bisg01 + bisg02 + bisg06 + # salud y cuidados
control_pap +
region + localidades_150 + # contexto de residencia
cant_componentes + tipo_hogar + # rcaracterísticas del hogar/familia
quintil_uc + condicion_actividad + # ingresos y empleo
nivel_actividad_fisica + bial08, data = train_data) # var relacionadas a estilos de vida
El modelo para el entrenamiento elegido en este caso fue un árbol de decisión. La librería de R y tipo de tarea seleccionada para ello fue rpart y clasificación, respectivamente.
tree_spec <- decision_tree() %>%
set_engine("rpart") %>%
set_mode("classification")
Controlo:
tree_spec
## Decision Tree Model Specification (classification)
##
## Computational engine: rpart
A continuación armo mi workflow: * Agrego la receta mediante add_recipe(), y luego, * agrego el modelo.
rf_wf <- workflow() %>%
add_recipe(rf_recipe) %>%
add_model(tree_spec)
Luego imprimo el workflow:
rf_wf
## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: decision_tree()
##
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 0 Recipe Steps
##
## ── Model ───────────────────────────────────────────────────────────────────────
## Decision Tree Model Specification (classification)
##
## Computational engine: rpart
En esta etapa, con la función fit(), realizo el entrenamiento del modelo.
set.seed(123)
tree_fit <- tree_spec %>%
fit(control_mamografia ~ . , data= train_data )
#imprimo el modelo
tree_fit
## parsnip model object
##
## n= 5597
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 5597 2451 1 (0.5620868 0.4379132)
## 2) control_pap=1 2915 302 1 (0.8963979 0.1036021) *
## 3) control_pap=2,99 2682 533 2 (0.1987323 0.8012677) *
Vemos que el modelo genera un árbol simple con dos nodos terminales, lo que sugiere una fuerte asociación entre la prácticas preventivas de control mamográfico y realización de PAP.
Agrego la columna de la predicción a los datos de TRAIN, mediante la función bind_cols, para ver la diferencia entre la predicción y el valor verdadero.
results <- train_data %>%
bind_cols(predict(tree_fit, train_data)%>%
rename(.pred_tree = .pred_class))
Controlo:
head(results)
## # A tibble: 6 × 18
## control_mamografia bhch04 bhch05 nivel_instruccion cobertura_salud
## <fct> <dbl> <fct> <fct> <fct>
## 1 1 65 2 7 1
## 2 1 84 4 7 1
## 3 1 63 6 6 1
## 4 1 59 2 5 1
## 5 1 60 4 6 1
## 6 1 58 4 7 1
## # ℹ 13 more variables: condicion_actividad <fct>, bisg01 <fct>, bisg02 <fct>,
## # bisg06 <fct>, control_pap <fct>, region <fct>, localidades_150 <fct>,
## # cant_componentes <fct>, tipo_hogar <fct>, quintil_uc <fct>,
## # nivel_actividad_fisica <fct>, bial08 <fct>, .pred_tree <fct>
Comparando la columna de la variable de estratificación (control_mamografias) con la de predicciones (.pred_tree) se observa que el modelo clasificó correctamente las primeras observaciones.
metrics(results, truth = control_mamografia, estimate = .pred_tree)
## # A tibble: 2 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.851
## 2 kap binary 0.700
La métrica accuracy obtenida con el modelo ajustado en train es bastante buena (>0,80).
Esta etapa es la predicción del modelo de ML:
results_test <- predict(tree_fit, new_data = test_data, type = "prob")%>%
bind_cols(predict(tree_fit, new_data = test_data, type = "class")%>%
rename(.pred_tree = .pred_class))%>%
bind_cols(test_data)
head(results_test)
## # A tibble: 6 × 20
## .pred_1 .pred_2 .pred_tree control_mamografia bhch04 bhch05 nivel_instruccion
## <dbl> <dbl> <fct> <fct> <dbl> <fct> <fct>
## 1 0.896 0.104 1 1 56 4 7
## 2 0.896 0.104 1 1 55 5 7
## 3 0.199 0.801 2 2 80 5 3
## 4 0.199 0.801 2 2 63 2 5
## 5 0.896 0.104 1 1 58 3 7
## 6 0.896 0.104 1 1 60 2 7
## # ℹ 13 more variables: cobertura_salud <fct>, condicion_actividad <fct>,
## # bisg01 <fct>, bisg02 <fct>, bisg06 <fct>, control_pap <fct>, region <fct>,
## # localidades_150 <fct>, cant_componentes <fct>, tipo_hogar <fct>,
## # quintil_uc <fct>, nivel_actividad_fisica <fct>, bial08 <fct>
.pred_1 y .pred_2 me muestra las probabilidades de clasificar cada dato como con y sin control de mamografía (respectivamente). Comparando las predicciones resultantes (.pred_tree) con la clase de la columna control_mamografia (clase observada), vemos que los datos predichos (6 primeras observaciones) se ajustan correctamente a lo observado.
Selecciono las métricas del modelo que me interesan:
selected_metrics <- metric_set(accuracy, sensitivity, specificity, precision, f_meas)
selected_metrics(results_test, truth = control_mamografia, estimate = .pred_tree)
## # A tibble: 5 × 3
## .metric .estimator .estimate
## <chr> <chr> <dbl>
## 1 accuracy binary 0.844
## 2 sensitivity binary 0.823
## 3 specificity binary 0.871
## 4 precision binary 0.891
## 5 f_meas binary 0.856
Vemos que el modelo obtuvo un buen valor de accuracy (exactitud), el cual indica que clasifica bien el 84% de los casos. Además, el modelo identifica correctamente 82% de los casos con control mamográfico (sensibilidad), y el 87% de los casos sin control mamográfico (especificidad). Asimismo, la alta precisión (89%) y el F1-score (0,86) indican un adecuado balance entre detección de casos y control de falsos positivos.
Calculo la matriz de confusión y la imprimo:
library(yardstick)
conf_matrix <- conf_mat(data = results_test,
truth = control_mamografia,
estimate = .pred_tree)
print(conf_matrix)
## Truth
## Prediction 1 2
## 1 648 79
## 2 139 534
Segun lo anterior, vemos que el modelo clasifica bastante bien.
Grafico la curva ROC:
rf_fit <- last_fit(rf_wf, split)
preds <- collect_predictions(rf_fit)
names(preds)
## [1] ".pred_class" ".pred_1" ".pred_2"
## [4] "id" "control_mamografia" ".row"
## [7] ".config"
preds %>%
roc_curve(control_mamografia, .pred_1) %>% #aca va las predicciones "positivas" (exito) (pred_1)
ggplot(aes(x = 1 - specificity, y = sensitivity)) +
geom_line(size = 1.5, color = "midnightblue") +
geom_abline(lty = 2, alpha = 0.5, color = "grey50", size = 1.2) +
coord_equal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
La gráfica es consistente con las métricas reportadas, mostrando un buen ajuste del modelo.
Ploteo las variables más importantes del MODELO, mediante la librería vip.
library(vip)
# Fit the model to the entire training data to get the final model
final_fit <- rf_wf %>%
fit(data = train_data)
# Plot variable importance
vip(pull_workflow_fit(final_fit)$fit, geom = "point")
La gráfica muestra que la variable más discriminante en la clasificación de acuerdo al control mamográfico es la práctica de PAP. Otras variables predictoras que le siguen en orden de importancia (decreciente) fueron: edad, nivel de instrucción, situación conyugal, condición de actividad y capacidad de movilidad (caminar).