1. Exploración Analítica de Datos

El análisis exploratorio de datos (EDA) constituye la primera etapa del proceso analítico y tiene como propósito conocer la estructura, calidad y comportamiento general de la información antes de realizar análisis estadísticos de mayor complejidad. En este informe se utiliza una base de datos obtenida por encuesta de vigilancia de factores de riesgo en población adulta de EE.UU, que contiene aproximadamente 401.958 registros y 279 variables en la versión oficial de 2020, y pertenece al Sistema de Vigilancia de Factores de Riesgo Conductuales (BRFSS), de donde se conservan solo un subconjunto de 19 variables como propósito de análisis del informe.

1.1 Pregunta central y objetivo

Pregunta central: ¿Cómo se acumulan los factores de riesgo modificables (tabaquismo, sedentarismo, obesidad y sueño insuficiente) y las barreras económicas de acceso según el ingreso y la edad, y cómo se relacionan con la mala salud autopercibida y la enfermedad coronaria en adultos de EE. UU.?

Objetivo del EDA: Describir la estructura y la calidad de la base, identificar los grupos con mayor carga de riesgo y traducir los hallazgos en prioridades de gestión del riesgo y prevención. Este capítulo 1 se enfoca en dónde se concentra el riesgo


1.2 Preparación del entorno y de los datos

1.2.1 Paquetes y paleta de colores

Se instalan solo los paquetes que falten y se cargan.

paquetes <- c("tidyverse", "skimr", "janitor", "plotly", "knitr")
faltantes_pkg <- paquetes[!paquetes %in% installed.packages()[, "Package"]]
if (length(faltantes_pkg) > 0) install.packages(faltantes_pkg)

library(tidyverse)
library(skimr)
library(janitor)
library(plotly)
library(knitr)
library(ggpubr)

# Paleta institucional del informe
paleta <- c(azul = "#1F4E79", naranja = "#E67E22", verde = "#2E8B57", gris = "#7F8C8D")

1.2.2 Importación y selección de variables

# Importación de la base de datos
ruta_archivo <- "brfss2020.csv"
base_completa <- read_csv(ruta_archivo, col_types = cols(.default = col_double()))

# Datos generales de la base completa 
n_registros  <- nrow(base_completa)
n_variables  <- ncol(base_completa)
tipos_var    <- table(sapply(base_completa, class))
n_estados    <- n_distinct(base_completa$`_STATE`)
faltantes_por_variable <- base_completa |>
  summarise(across(everything(), ~ mean(is.na(.x)) * 100)) |>
  pivot_longer(everything(), names_to = "variable", values_to = "pct_faltante")

# Visualización inicial
head(base_completa)
# Nombres originales -> nombres en español
nombres_es <- c(
  estado_eeuu            = "_STATE",   sexo                   = "SEXVAR",
  salud_general          = "GENHLTH",  dias_mala_salud_fisica = "PHYSHLTH",
  dias_mala_salud_mental = "MENTHLTH", tiene_seguro           = "HLTHPLN1",
  barrera_costo          = "MEDCOST",  ultimo_chequeo         = "CHECKUP1",
  actividad_fisica       = "EXERANY2", horas_sueno            = "SLEPTIM1",
  tabaquismo             = "_SMOKER3", enf_coronaria          = "_MICHD",
  acv                    = "CVDSTRK3", diabetes               = "DIABETE4",
  imc                    = "_BMI5",    edad                   = "_AGE80",
  ingreso                = "_INCOMG",  educacion              = "_EDUCAG",
  peso_muestral          = "_LLCPWT"
)

base_sel <- base_completa |>
  select(all_of(nombres_es)) |>   
  clean_names()

rm(base_completa); invisible(gc())  

1.2.3 Normalización al español (recodificación)

Se convierten en NA (dato faltante) o en su valor real, cada categoría se traduce al español y se crean los indicadores de riesgo (1 = presente, 0 = ausente)

dias_validos <- function(x) case_when(x == 88 ~ 0, x >= 1 & x <= 30 ~ x, TRUE ~ NA_real_)

base_eda <- base_sel |>
  mutate(
    sexo             = factor(sexo, levels = 1:2, labels = c("Hombre", "Mujer")),
    salud_general    = factor(salud_general, levels = 1:5,
                              labels = c("Excelente", "Muy buena", "Buena", "Regular", "Mala")),
    dias_mala_salud_fisica = dias_validos(dias_mala_salud_fisica),
    dias_mala_salud_mental = dias_validos(dias_mala_salud_mental),
    tiene_seguro     = factor(tiene_seguro,  levels = 1:2, labels = c("Sí", "No")),
    barrera_costo    = factor(barrera_costo, levels = 1:2, labels = c("Sí", "No")),
    ultimo_chequeo   = factor(ultimo_chequeo, levels = c(1, 2, 3, 4, 8),
                              labels = c("< 1 año", "1-2 años", "2-5 años", "5+ años", "Nunca")),
    actividad_fisica = factor(actividad_fisica, levels = 1:2, labels = c("Sí", "No")),
    horas_sueno      = if_else(horas_sueno <= 24, horas_sueno, NA_real_),
    tabaquismo       = factor(tabaquismo, levels = 1:4,
                              labels = c("Fuma diario", "Fuma algunos días", "Exfumador", "Nunca fumó")),
    enf_coronaria    = factor(enf_coronaria, levels = 1:2, labels = c("Sí", "No")),
    acv              = factor(acv, levels = 1:2, labels = c("Sí", "No")),
    diabetes         = factor(diabetes, levels = 1:4,
                              labels = c("Sí", "Solo en embarazo", "No", "Prediabetes")),
    imc              = imc / 100,   # la base trae el IMC con 2 decimales implícitos
    ingreso          = factor(ingreso, levels = 1:5, ordered = TRUE,
                              labels = c("< 15 mil", "15-25 mil", "25-35 mil", "35-50 mil", "50 mil o más")),
    educacion        = factor(educacion, levels = 1:4,
                              labels = c("Sin bachillerato", "Bachiller", "Superior incompleta", "Universitario")),
    grupo_edad       = cut(edad, breaks = c(17, 34, 49, 64, 80),
                           labels = c("18-34", "35-49", "50-64", "65+")),
    # Indicadores de riesgo (1 = presente)
    fuma_actual      = case_when(tabaquismo %in% c("Fuma diario", "Fuma algunos días") ~ 1,
                                 tabaquismo %in% c("Exfumador", "Nunca fumó") ~ 0),
    sedentarismo     = case_when(actividad_fisica == "No" ~ 1, actividad_fisica == "Sí" ~ 0),
    obesidad         = if_else(imc >= 30, 1, 0),
    sueno_corto      = if_else(horas_sueno < 7, 1, 0),
    n_factores_riesgo = fuma_actual + sedentarismo + obesidad + sueno_corto,
    mala_salud       = case_when(salud_general %in% c("Regular", "Mala") ~ 1,
                                 !is.na(salud_general) ~ 0)
  )

1.2.4 Diccionario de variables

diccionario <- tibble(
  variable_original = unname(nombres_es),
  variable_espanol  = names(nombres_es),
  descripcion = c(
    "Estado o territorio de residencia", "Sexo del encuestado",
    "Salud general autopercibida", "Días con mala salud física (últimos 30)",
    "Días con mala salud mental (últimos 30)", "Tiene algún seguro de salud",
    "No consultó al médico por costo (último año)", "Tiempo desde el último chequeo",
    "Actividad física fuera del trabajo (último mes)", "Horas de sueño en 24 h",
    "Condición de fumador", "Enfermedad coronaria o infarto (calculada)",
    "Antecedente de ACV", "Diagnóstico de diabetes", "Índice de masa corporal (kg/m²)",
    "Edad en años (tope 80)", "Ingreso anual del hogar (USD)", "Nivel educativo",
    "Peso muestral (factor de expansión)"),
  tipo = c("Categórica", "Categórica", "Ordinal", "Numérica", "Numérica", "Categórica",
           "Categórica", "Ordinal", "Categórica", "Numérica", "Categórica", "Categórica",
           "Categórica", "Categórica", "Numérica", "Numérica", "Ordinal", "Ordinal", "Numérica"),
  dimension = c("Identificación", "Demografía", "Resultado en salud", "Resultado en salud",
                "Resultado en salud", "Acceso / financiero", "Acceso / financiero", "Uso de servicios",
                "Riesgo modificable", "Riesgo modificable", "Riesgo modificable", "Enfermedad crónica",
                "Enfermedad crónica", "Enfermedad crónica", "Riesgo modificable", "Demografía",
                "Socioeconómica / financiero", "Socioeconómica", "Diseño muestral")
)
diccionario

1.3 Comprensión de la base

1.3.1 Registros, variables y unidad de análisis

tibble(
  Indicador = c("Número de registros", "Número de variables", "Estados/territorios",
                "Variables seleccionadas para el EDA"),
  Valor = c(format(n_registros, big.mark = "."), n_variables, n_estados, ncol(base_sel))
)
  • Unidad de análisis: cada fila es un adulto (18 años o más) encuestado por teléfono en 2020. No son pacientes ni atenciones: es población general, lo que la hace útil para prevención y no para medir producción de servicios.
  • La base tiene 401.958 registros, 279 variables y cubre 53 estados y territorios (50 estados, Distrito de Columbia, Guam y Puerto Rico).
  • Datos financieros: la base no contiene costos ni facturación. Se usan tres variables como aproximación (proxy, variable indirecta) de la dimensión financiera: ingreso del hogar, haber dejado de consultar por costo y tener seguro de salud.
  • Ponderación: el BRFSS es una encuesta compleja; para cifras oficiales debe usarse peso_muestral. Este EDA reporta valores sin ponderar; la diferencia es pequeña (mala salud: 15,4 % sin ponderar vs. 14,7 % ponderado; obesidad: 32,0 % vs. 31,9 %).

1.3.2 Tipos de variables

tipos_var                       # tipo en el archivo original
## 
## numeric 
##     279
table(sapply(base_eda, function(x) class(x)[1]))   # tipo después de normalizar
## 
##  factor numeric ordered 
##      12      13       1

1.3.3 Datos faltantes

faltantes_por_variable |>
  mutate(rango = cut(pct_faltante, breaks = c(-Inf, 0, 10, 50, Inf),
                     labels = c("Sin faltantes", "Hasta 10 %", "10 % a 50 %", "Más de 50 %"))) |>
  count(rango, name = "n_variables")
# Faltantes en las variables del EDA, después de convertir códigos 7/9/77/99 a NA
base_eda |>
  summarise(across(c(salud_general, dias_mala_salud_fisica, dias_mala_salud_mental, barrera_costo,
                     horas_sueno, tabaquismo, enf_coronaria, imc, ingreso, n_factores_riesgo),
                   ~ round(mean(is.na(.x)) * 100, 1))) |>
  pivot_longer(everything(), names_to = "variable", values_to = "pct_faltante") |>
  arrange(desc(pct_faltante))

Interpretación:

  • En la base completa, el 45,7 % de las celdas está vacío y 142 de las 279 variables tienen más de 50 % de faltantes; ningún registro está completo. Esto no es un error: el BRFSS usa módulos opcionales (cada estado decide cuáles aplica) y preguntas condicionadas (por ejemplo, la tamización de próstata solo se hace a hombres).
  • Implicación: no se debe eliminar filas con cualquier faltante, porque se perdería toda la base.
  • Alerta: el ingreso tiene cerca de 20 % sin dato; es la variable financiera clave y su faltante no es aleatorio (las personas suelen negarse a declarar ingresos). El IMC falta en ~10 % y el número de factores de riesgo en ~13 %.

1.4 Calidad y distribución

1.4.1 Estructura de los datos

str(base_eda)
## tibble [401,958 × 26] (S3: tbl_df/tbl/data.frame)
##  $ estado_eeuu           : num [1:401958] 1 1 1 1 1 1 1 1 1 1 ...
##  $ sexo                  : Factor w/ 2 levels "Hombre","Mujer": 2 2 2 2 2 1 2 2 2 2 ...
##  $ salud_general         : Factor w/ 5 levels "Excelente","Muy buena",..: 2 3 3 1 2 4 3 4 2 4 ...
##  $ dias_mala_salud_fisica: num [1:401958] 3 0 0 0 0 20 0 15 28 6 ...
##  $ dias_mala_salud_mental: num [1:401958] 30 0 0 0 0 30 0 10 0 0 ...
##  $ tiene_seguro          : Factor w/ 2 levels "Sí","No": 2 1 1 1 1 1 1 1 1 1 ...
##  $ barrera_costo         : Factor w/ 2 levels "Sí","No": 1 1 2 2 2 2 2 2 2 2 ...
##  $ ultimo_chequeo        : Factor w/ 5 levels "< 1 año","1-2 años",..: 4 1 1 2 1 2 1 1 1 1 ...
##  $ actividad_fisica      : Factor w/ 2 levels "Sí","No": 1 1 1 2 1 1 2 1 1 2 ...
##  $ horas_sueno           : num [1:401958] 5 7 7 6 7 8 6 6 8 12 ...
##  $ tabaquismo            : Factor w/ 4 levels "Fuma diario",..: 1 NA 4 4 4 3 4 1 4 3 ...
##  $ enf_coronaria         : Factor w/ 2 levels "Sí","No": 2 2 2 2 2 2 2 2 2 1 ...
##  $ acv                   : Factor w/ 2 levels "Sí","No": 2 2 2 2 1 2 2 2 2 2 ...
##  $ diabetes              : Factor w/ 4 levels "Sí","Solo en embarazo",..: 1 3 3 3 3 1 3 1 3 3 ...
##  $ imc                   : num [1:401958] 16.6 29.2 NA NA 20.3 ...
##  $ edad                  : num [1:401958] 56 65 65 80 80 66 75 69 41 78 ...
##  $ ingreso               : Ord.factor w/ 5 levels "< 15 mil"<"15-25 mil"<..: 1 NA 5 NA NA 3 4 3 4 NA ...
##  $ educacion             : Factor w/ 4 levels "Sin bachillerato",..: 4 4 3 2 4 2 2 2 4 1 ...
##  $ peso_muestral         : num [1:401958] 284 171 1334 1297 455 ...
##  $ grupo_edad            : Factor w/ 4 levels "18-34","35-49",..: 3 4 4 4 4 4 4 4 2 4 ...
##  $ fuma_actual           : num [1:401958] 1 NA 0 0 0 0 0 1 0 0 ...
##  $ sedentarismo          : num [1:401958] 0 0 0 1 0 0 1 0 0 1 ...
##  $ obesidad              : num [1:401958] 0 0 NA NA 0 0 0 0 0 0 ...
##  $ sueno_corto           : num [1:401958] 1 0 0 1 0 0 1 1 0 0 ...
##  $ n_factores_riesgo     : num [1:401958] 2 NA NA NA 0 0 2 2 0 1 ...
##  $ mala_salud            : num [1:401958] 0 0 0 0 0 1 0 1 0 1 ...
summary(select(base_eda, edad, imc, horas_sueno, dias_mala_salud_fisica, dias_mala_salud_mental))
##       edad            imc         horas_sueno     dias_mala_salud_fisica
##  Min.   :18.00   Min.   :12.02   Min.   : 1.000   Min.   : 0.000        
##  1st Qu.:40.00   1st Qu.:23.99   1st Qu.: 6.000   1st Qu.: 0.000        
##  Median :57.00   Median :27.32   Median : 7.000   Median : 0.000        
##  Mean   :54.43   Mean   :28.31   Mean   : 7.097   Mean   : 3.458        
##  3rd Qu.:69.00   3rd Qu.:31.38   3rd Qu.: 8.000   3rd Qu.: 2.000        
##  Max.   :80.00   Max.   :98.43   Max.   :24.000   Max.   :30.000        
##                  NA's   :41357   NA's   :4694     NA's   :8691          
##  dias_mala_salud_mental
##  Min.   : 0.000        
##  1st Qu.: 0.000        
##  Median : 0.000        
##  Mean   : 3.916        
##  3rd Qu.: 3.000        
##  Max.   :30.000        
##  NA's   :7929

1.4.2 Variables numéricas

Resume en una tabla: faltantes, media, desviación estándar, percentiles y un mini histograma.

base_eda |>
  select(edad, imc, horas_sueno, dias_mala_salud_fisica, dias_mala_salud_mental, n_factores_riesgo) |>
  skim()
Data summary
Name select(…)
Number of rows 401958
Number of columns 6
_______________________
Column type frequency:
numeric 6
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
edad 0 1.00 54.43 17.67 18.00 40.00 57.00 69.00 80.00 ▃▅▆▇▇
imc 41357 0.90 28.31 6.38 12.02 23.99 27.32 31.38 98.43 ▇▅▁▁▁
horas_sueno 4694 0.99 7.10 1.47 1.00 6.00 7.00 8.00 24.00 ▁▇▁▁▁
dias_mala_salud_fisica 8691 0.98 3.46 8.08 0.00 0.00 0.00 2.00 30.00 ▇▁▁▁▁
dias_mala_salud_mental 7929 0.98 3.92 8.05 0.00 0.00 0.00 3.00 30.00 ▇▁▁▁▁
n_factores_riesgo 53108 0.87 1.00 0.95 0.00 0.00 1.00 2.00 4.00 ▇▇▅▂▁

Interpretación:

  • Edad: mediana 57 años; población adulta mayor que el promedio, consistente con encuestas telefónicas.
  • IMC: mediana 27,3 kg/m² (rango de sobrepeso); la mitad central de la población está entre 24,0 y 31,4.
  • Horas de sueño: mediana 7 h, pero el 30,6 % duerme menos de 7 h (sueño insuficiente).
  • Días de mala salud física y mental: la mediana es 0, pero el promedio es mayor. Son variables asimétricas (la mayoría reporta 0 días y un grupo pequeño reporta muchos), por eso la mediana describe mejor al grupo típico que la media.

1.4.3 Variables categóricas

variables_cat <- c("sexo", "salud_general", "tiene_seguro", "barrera_costo", "ultimo_chequeo",
                   "actividad_fisica", "tabaquismo", "diabetes", "enf_coronaria", "ingreso")

base_eda |>
  select(all_of(variables_cat)) |>
  mutate(across(everything(), as.character)) |>
  pivot_longer(everything(), names_to = "variable", values_to = "categoria") |>
  mutate(categoria = replace_na(categoria, "Sin dato")) |>
  count(variable, categoria) |>
  group_by(variable) |>
  mutate(porcentaje = round(100 * n / sum(n), 1)) |>
  ungroup()

Se genera Tabla cruzada de la dimensión financiera: porcentaje que no consultó por costo según ingreso (con janitor::tabyl).

base_eda |>
  tabyl(ingreso, barrera_costo, show_na = FALSE) |>
  adorn_percentages("row") |>
  adorn_pct_formatting(digits = 1) |>
  adorn_ns()

Interpretación:

  • Salud general: 15,4 % se califica en salud regular o mala; este es el indicador resultado más completo de la base.
  • Acceso: 8,5 % no tiene seguro y 8,6 % dejó de consultar por costo; 77 % tuvo chequeo en el último año.
  • Riesgo: 14 % fuma actualmente, 24 % es sedentario y 13 % reporta diabetes.
  • Gradiente financiero: en hogares con ingreso < 15 mil USD, 18,2 % no consultó por costo, frente a 4,8 % en hogares con 50 mil o más (casi 4 veces más). El costo es una barrera concentrada, no general.

1.4.4 Valores extremos (boxplot)

Un boxplot (diagrama de caja) muestra la mediana (línea central), el 50 % central de los datos (la caja) y los valores atípicos (puntos fuera de 1,5 veces el rango intercuartílico, es decir, la distancia entre el percentil 25 y el 75).

datos_box <- base_eda |>
  select(`IMC (kg/m²)` = imc, `Horas de sueño` = horas_sueno,
         `Días mala salud física` = dias_mala_salud_fisica,
         `Días mala salud mental` = dias_mala_salud_mental) |>
  pivot_longer(everything(), names_to = "variable", values_to = "valor") |>
  drop_na()

ggplot(datos_box, aes(x = variable, y = valor, fill = variable)) +
  geom_boxplot(outlier.alpha = 0.05, outlier.colour = paleta["gris"], width = 0.5) +
  facet_wrap(~ variable, scales = "free") +
  scale_fill_manual(values = unname(paleta)) +
  labs(title = "Valores extremos en variables numéricas clave", x = NULL, y = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none", axis.text.x = element_blank())

# Cuantificación de atípicos
datos_box |>
  group_by(variable) |>
  summarise(mediana = median(valor), q1 = quantile(valor, 0.25), q3 = quantile(valor, 0.75),
            limite_superior = q3 + 1.5 * (q3 - q1),
            pct_atipicos = round(100 * mean(valor < q1 - 1.5 * (q3 - q1) | valor > limite_superior), 1),
            maximo = max(valor))

Interpretación:

  • IMC: 3,3 % de atípicos con máximo de 98 kg/m²; valores > 70 son biológicamente posibles pero raros y pueden ser errores de digitación de peso o talla. Decisión: se conservan para el EDA y se marcan para revisión antes de cualquier modelamiento.
  • Sueño: 1,5 % de atípicos (personas que reportan 1-3 h o más de 12 h). Pueden reflejar trabajo nocturno o enfermedad.
  • Días de mala salud física (15,1 %) y mental (16,2 %): aquí los “atípicos” no son errores: son las personas con mayor carga de enfermedad. Desde la gestión del riesgo, estos extremos son la población objetivo, no ruido que se deba eliminar.

1.4.5 Distribución del IMC

ggplot(filter(base_eda, imc <= 60), aes(x = imc)) +
  geom_histogram(binwidth = 1, fill = paleta["azul"], colour = "white") +
  geom_vline(xintercept = c(25, 30), linetype = "dashed", colour = paleta["naranja"], linewidth = 0.8) +
  annotate("text", x = c(25, 30), y = Inf, vjust = 1.5, hjust = -0.1,
           label = c("Sobrepeso", "Obesidad"), colour = paleta["naranja"]) +
  labs(title = "Distribución del índice de masa corporal",
       subtitle = "Se muestran valores hasta 60 kg/m² para facilitar la lectura",
       x = "IMC (kg/m²)", y = "Número de personas") +
  theme_minimal(base_size = 12)

Interpretación: la distribución tiene sesgo a la derecha (cola larga hacia valores altos). El pico está en sobrepeso y el 32 % de los adultos con dato tiene obesidad (IMC ≥ 30). La obesidad es el factor de riesgo modificable más frecuente de los cuatro analizados.


1.5 Relaciones entre variables

Los siguientes gráficos son interactivos: al pasar el cursor se ve el valor exacto

1.5.1 Gradiente socioeconómico del riesgo y del acceso

  • Variables resultado: mala salud autopercibida, no consultar por costo, no tener seguro, sedentarismo (%).
  • Variable de comparación: ingreso anual del hogar.
gradiente_ingreso <- base_eda |>
  filter(!is.na(ingreso)) |>
  group_by(ingreso) |>
  summarise(`Mala salud autopercibida` = mean(mala_salud, na.rm = TRUE) * 100,
            `No consultó por costo`    = mean(barrera_costo == "Sí", na.rm = TRUE) * 100,
            `Sin seguro de salud`      = mean(tiene_seguro == "No", na.rm = TRUE) * 100,
            `Sedentarismo`             = mean(sedentarismo, na.rm = TRUE) * 100) |>
  pivot_longer(-ingreso, names_to = "indicador", values_to = "porcentaje")

plot_ly(gradiente_ingreso, x = ~ingreso, y = ~porcentaje, color = ~indicador,
        colors = unname(paleta), type = "bar",
        text = ~paste0(round(porcentaje, 1), "%"), textposition = "outside",
        hovertemplate = "%{x}<br>%{y:.1f}%<extra></extra>") |>
  layout(title = "A menor ingreso, mayor riesgo y menor acceso",
         barmode = "group",
         xaxis = list(title = "Ingreso anual del hogar (USD)"),
         yaxis = list(title = "% de adultos", range = c(0, 50)),
         legend = list(orientation = "h", y = -0.2))

Interpretación: los cuatro indicadores descienden de forma escalonada a medida que sube el ingreso (gradiente). En el grupo de menor ingreso, 38,9 % reporta mala salud, frente a 7,1 % en el de mayor ingreso (5,5 veces más). El sedentarismo pasa de 41,7 % a 14,5 %, la falta de seguro de 16,2 % a 4,1 % y la barrera de costo de 18,2 % a 4,8 %.

Para la gestión: el ingreso funciona como un criterio de focalización: un mismo recurso preventivo tiene más impacto potencial en los quintiles bajos, donde coinciden más riesgo y menos acceso.

1.5.2 Superficie 3D: edad × factores de riesgo × enfermedad coronaria

  • Variable resultado: prevalencia (porcentaje de personas que tienen la condición) de enfermedad coronaria o infarto.
  • Variables de comparación: grupo de edad y número de factores de riesgo modificables (0 a 4).
matriz_coronaria <- base_eda |>
  filter(!is.na(grupo_edad), !is.na(n_factores_riesgo), !is.na(enf_coronaria)) |>
  group_by(grupo_edad, n_factores_riesgo) |>
  summarise(prevalencia = mean(enf_coronaria == "Sí") * 100, .groups = "drop") |>
  pivot_wider(names_from = n_factores_riesgo, values_from = prevalencia)

z_prevalencia <- as.matrix(matriz_coronaria[, -1])   # filas = edad, columnas = n.º factores

plot_ly(x = 0:4, y = 1:4, z = z_prevalencia, type = "surface",
        colorscale = list(c(0, paleta[["verde"]]), c(0.5, paleta[["gris"]]), c(1, paleta[["naranja"]])),
        colorbar = list(title = "% coronaria"),
        hovertemplate = "Factores de riesgo: %{x}<br>Grupo de edad: %{y}<br>Prevalencia: %{z:.1f}%<extra></extra>") |>
  layout(title = "Riesgo acumulado: la edad no modificable y los factores modificables se suman",
         scene = list(
           xaxis = list(title = "N.º factores de riesgo"),
           yaxis = list(title = "Edad", tickvals = 1:4, ticktext = levels(base_eda$grupo_edad)),
           zaxis = list(title = "% enfermedad coronaria", color = paleta[["azul"]]),
           camera = list(eye = list(x = 1.6, y = -1.6, z = 0.9))))

Interpretación: la superficie sube en las dos direcciones. Con la edad (factor no modificable) y con el número de factores de riesgo (factores modificables).

  • En 50-64 años, la prevalencia pasa de 4,1 % sin factores de riesgo a 21,7 % con los cuatro (más de 5 veces).
  • En 65+, va de 13,3 % a 29,8 %.
  • Una persona de 50-64 años con 3 factores (15,5 %) tiene una prevalencia mayor que una de 65+ sin factores (13,3 %).
  • El 27,9 % de los adultos con dato completo acumula 2 o más factores.

Para la gestión: el número de factores de riesgo es una regla simple de estratificación (clasificar a la población por nivel de riesgo) que se puede aplicar en el primer nivel de atención. La ventana de mayor rendimiento preventivo está en 35-64 años, donde los factores modificables todavía cambian mucho el resultado.

Nota: es una asociación transversal (medida en un solo momento), no causal; las celdas de 4 factores en jóvenes tienen pocos casos.

1.5.3 Burbujas 3D: ingreso × edad × mala salud, con barrera de costo

  • Variable resultado: % de mala salud autopercibida (altura).
  • Variables de comparación: ingreso, grupo de edad y % que no consultó por costo (color). El tamaño de la burbuja indica cuántas personas hay en el grupo.
burbujas <- base_eda |>
  filter(!is.na(ingreso), !is.na(grupo_edad)) |>
  group_by(ingreso, grupo_edad) |>
  summarise(pct_mala_salud = mean(mala_salud, na.rm = TRUE) * 100,
            pct_costo      = mean(barrera_costo == "Sí", na.rm = TRUE) * 100,
            n = n(), .groups = "drop") |>
  mutate(tamano = scales::rescale(sqrt(n), to = c(8, 30)))

plot_ly(burbujas, x = ~as.numeric(ingreso), y = ~as.numeric(grupo_edad), z = ~pct_mala_salud,
        type = "scatter3d", mode = "markers",
        marker = list(size = ~tamano, color = ~pct_costo, opacity = 0.85,
                      colorscale = list(c(0, paleta[["azul"]]), c(0.5, paleta[["gris"]]), c(1, paleta[["naranja"]])),
                      colorbar = list(title = "% no consultó<br>por costo"), line = list(width = 0)),
        text = ~paste0("Ingreso: ", ingreso, "<br>Edad: ", grupo_edad,
                       "<br>Mala salud: ", round(pct_mala_salud, 1), "%",
                       "<br>Barrera de costo: ", round(pct_costo, 1), "%",
                       "<br>Personas: ", format(n, big.mark = ".")),
        hoverinfo = "text") |>
  layout(title = "Dónde se concentran la mala salud y la barrera económica",
         scene = list(
           xaxis = list(title = "Ingreso (USD)", tickvals = 1:5, ticktext = levels(base_eda$ingreso)),
           yaxis = list(title = "Edad", tickvals = 1:4, ticktext = levels(base_eda$grupo_edad)),
           zaxis = list(title = "% mala salud")))

Interpretación:

  • El punto más alto está en 50-64 años con ingreso < 15 mil USD: casi la mitad (49,7 %) reporta mala salud, frente a 7,1 % en la misma edad con ingreso alto (7 veces más).
  • La barrera de costo (color naranja) se concentra en menores de 65 años de bajo ingreso (21-24 %). En 65+ cae a 2-10 %, lo que es coherente con la cobertura pública universal desde los 65 años en EE. UU. (Medicare).
  • Las burbujas grandes (más población) están en ingreso alto, mientras que la mayor carga está en grupos pequeños. Por eso los promedios generales esconden el problema.

Para la gestión: el segmento adulto 35-64 de bajo ingreso combina la peor salud y la mayor barrera económica. Es el grupo prioritario para programas de prevención con subsidio o con atención sin copago.

1.5.4 Correlación preliminar entre variables

Se usa la correlación de Spearman (mide si dos variables aumentan o disminuyen juntas, basada en el orden de los valores; va de -1 a +1). Es adecuada porque varias variables son asimétricas u ordinales.

matriz_cor <- base_eda |>
  transmute(Edad = edad, IMC = imc, `Horas sueño` = horas_sueno,
            `Días salud física` = dias_mala_salud_fisica, `Días salud mental` = dias_mala_salud_mental,
            Ingreso = as.numeric(ingreso), `Salud general (1=excelente)` = as.numeric(salud_general),
            `N.º factores riesgo` = n_factores_riesgo) |>
  cor(method = "spearman", use = "pairwise.complete.obs")

plot_ly(x = colnames(matriz_cor), y = rownames(matriz_cor), z = matriz_cor, type = "heatmap",
        zmin = -1, zmax = 1,
        colorscale = list(c(0, paleta[["naranja"]]), c(0.5, "#F2F2F2"), c(1, paleta[["azul"]])),
        text = round(matriz_cor, 2), texttemplate = "%{text}",
        hovertemplate = "%{x} vs %{y}: %{z:.2f}<extra></extra>") |>
  layout(title = "Matriz de correlación (Spearman)",
         xaxis = list(tickangle = -35), yaxis = list(autorange = "reversed"))

Interpretación (en salud general, valores altos significan peor salud):

  • La correlación más fuerte con la mala salud general es con los días de mala salud física (0,41), seguida del número de factores de riesgo (0,34) y del ingreso (-0,30): a más ingreso, mejor salud.
  • Salud física y mental se mueven juntas (0,29), lo que apoya un abordaje integral y no por separado.
  • La edad se asocia negativamente con los días de mala salud mental (-0,22): los adultos jóvenes reportan más carga mental.
  • Las correlaciones del número de factores con IMC (0,49) y sueño (-0,46) son altas por construcción (son parte del indicador), no son hallazgos.
  • En general las correlaciones son débiles a moderadas (< 0,5): ninguna variable explica por sí sola la salud. Esto justifica pasar, en una fase posterior, a un modelo que combine varias variables.

1.6 Hallazgos principales y su influencia en la toma de decisiones

# Hallazgo Implicación para la gestión del riesgo
1 Base grande (401.958 adultos), pero con 45,7 % de celdas vacías por diseño modular; ingreso sin dato en ~20 %. Analizar por módulo y documentar faltantes; no eliminar registros en bloque. Mejorar la captura de variables socioeconómicas en los propios sistemas de información.
2 Gradiente socioeconómico claro: la mala salud es 5,5 veces mayor en ingreso bajo que en alto; la barrera de costo, casi 4 veces mayor. Focalizar recursos de prevención por nivel socioeconómico; el costo de bolsillo es una barrera de acceso medible y modificable.
3 Los factores de riesgo modificables se acumulan (28 % tiene 2 o más) y la prevalencia coronaria crece con cada factor adicional en todos los grupos de edad. Usar el conteo de factores como herramienta sencilla de estratificación del riesgo en atención primaria.
4 El segmento de 35-64 años de bajo ingreso concentra la peor salud y la mayor barrera económica, mientras que en 65+ la barrera cae. Priorizar este segmento en programas de detección temprana y control de riesgo; la cobertura financiera reduce la barrera de acceso.
5 Los valores extremos en días de mala salud (~15 %) representan población de alta carga, no errores. Identificar y gestionar a los usuarios de alta carga (gestión de casos); revisar solo los extremos de IMC biológicamente improbables.
6 Las correlaciones son débiles a moderadas y multifactoriales. Justifica la siguiente fase: modelamiento con varias variables y uso del peso muestral.

2. Probabilidad y Estadística Aplicada (Descriptiva e Inferencial)

2.1 Propósito y hoja de ruta

Este capítulo responde cuánto, con qué certeza y para qué alcanzan los recursos. Para eso se aplican tres herramientas:

Herramienta Pregunta de gestión que responde Funciones de R
Estadística descriptiva ¿Cómo se comportan los indicadores clave y cuánto varían? mean(), median(), sd(), var(), range(), summary()
Probabilidad y simulación ¿Qué tan probable es un evento y cuántos recursos debo prever? sample(), table(), runif(), rbinom(), dpois(), rnorm(), rexp(), dnorm(), pnorm(), dexp(), pexp()
Estadística inferencial ¿Las diferencias entre grupos son reales o pueden deberse al azar? t.test(), aov(), chisq.test(), prop.test()

2.1.1 Preparación

if (!"ggpubr" %in% installed.packages()[, "Package"]) install.packages("ggpubr")
library(ggpubr)

# semilla: hace que las simulaciones den el mismo resultado cada vez que se genera el informe
set.seed(2026)   

# Función de apoyo para mostrar porcentajes en el texto
pct <- function(x, dec = 1) paste0(format(round(100 * x, dec), decimal.mark = ","), " %")

2.2 Estadística descriptiva de los indicadores clave

2.2.1 Tendencia central y dispersión

  • Media: Promedio.
  • Mediana: Valor del centro; la mitad de las personas está por debajo.
  • Desviación estándar: Cuánto se alejan en promedio los datos de la media.
  • Varianza: Desviación estándar al cuadrado.
  • Rango: Valores mínimo y máximo.
  • Coeficiente de variación (CV): Desviación estándar dividida entre la media, en %. Permite comparar la variabilidad de indicadores con unidades distintas.
indicadores <- base_eda |>
  select(Edad = edad, IMC = imc, `Horas de sueño` = horas_sueno,
         `Días mala salud física` = dias_mala_salud_fisica,
         `Días mala salud mental` = dias_mala_salud_mental)

tabla_descriptiva <- map_dfr(indicadores, function(x) {
  x <- na.omit(x)
  tibble(n = length(x),
         Media = mean(x), Mediana = median(x),
         `Desv. estándar` = sd(x), Varianza = var(x),
         Mínimo = range(x)[1], Máximo = range(x)[2],
         `CV (%)` = sd(x) / mean(x) * 100)
}, .id = "Indicador") |>
  mutate(across(where(is.numeric), ~ round(.x, 2)))

tabla_descriptiva
summary(indicadores)   # resumen completo: mínimo, cuartiles, media, máximo y faltantes (NA's)
##       Edad            IMC        Horas de sueño   Días mala salud física
##  Min.   :18.00   Min.   :12.02   Min.   : 1.000   Min.   : 0.000        
##  1st Qu.:40.00   1st Qu.:23.99   1st Qu.: 6.000   1st Qu.: 0.000        
##  Median :57.00   Median :27.32   Median : 7.000   Median : 0.000        
##  Mean   :54.43   Mean   :28.31   Mean   : 7.097   Mean   : 3.458        
##  3rd Qu.:69.00   3rd Qu.:31.38   3rd Qu.: 8.000   3rd Qu.: 2.000        
##  Max.   :80.00   Max.   :98.43   Max.   :24.000   Max.   :30.000        
##                  NA's   :41357   NA's   :4694     NA's   :8691          
##  Días mala salud mental
##  Min.   : 0.000        
##  1st Qu.: 0.000        
##  Median : 0.000        
##  Mean   : 3.916        
##  3rd Qu.: 3.000        
##  Max.   :30.000        
##  NA's   :7929

Interpretación:

  • IMC: media 28,3 y mediana 27,3 kg/m². Que la media quede por encima de la mediana indica asimetría (cola larga hacia valores altos): unos pocos IMC muy altos suben el promedio.
  • Horas de sueño: media 7,1 h con CV de 21 %; es el indicador más estable.
  • Días de mala salud física (media 3,5; mediana 0) y mental (media 3,9; mediana 0): tienen un CV mayor a 200 %, el más alto de todos. La mayoría reporta 0 días y una minoría reporta muchos.
  • Para la gestión: en indicadores tan dispersos, el promedio por sí solo engaña. Conviene reportar mediana y percentil 75 y gestionar por separado al grupo de alta carga, como en un modelo de gestión de casos.

2.2.2 Resumen por nivel de ingreso

tapply(base_eda$dias_mala_salud_fisica, base_eda$ingreso, summary)
## $`< 15 mil`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   0.000   0.000   7.906  15.000  30.000    1050 
## 
## $`15-25 mil`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   0.000   0.000   5.498   5.000  30.000    1357 
## 
## $`25-35 mil`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   0.000   0.000   4.155   3.000  30.000     760 
## 
## $`35-50 mil`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   0.000   0.000   3.355   2.000  30.000     708 
## 
## $`50 mil o más`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   0.000   0.000   2.115   0.000  30.000    1713

Interpretación: la media de días con mala salud física baja de forma escalonada: 7,9 días en ingreso < 15 mil USD, luego 5,5, 4,2, 3,4 y 2,1 días en ingreso de 50 mil o más. El percentil 75 (3er cuartil) pasa de 15 días a 0: en el grupo de menor ingreso, uno de cada cuatro adultos pierde la mitad del mes por mala salud física. Este resumen anticipa la prueba ANOVA de la sección 2.6.2.


2.3 Probabilidad empírica con datos reales

La probabilidad empírica es la proporción de veces que un evento ocurre en los datos observados. Por ejemplo, si 32 de cada 100 adultos tienen obesidad, la probabilidad de obesidad es 0,32.

2.3.1 Espacio muestral y selección aleatoria: brigada de tamización

Caso práctico: una aseguradora organiza brigadas de tamización (detección temprana) cardiovascular. En cada brigada se cita al azar a 50 adultos.

  • El espacio muestral (todos los resultados posibles) para cada persona es su número de factores de riesgo: {0, 1, 2, 3, 4}.
  • El evento A es que una persona tenga 2 o más factores de riesgo y requiera consulta de seguimiento.

¿Cuántas consultas de seguimiento debe reservar cada brigada?

poblacion_tamizable <- base_eda |> filter(!is.na(n_factores_riesgo)) |> pull(n_factores_riesgo)

espacio_muestral <- sort(unique(poblacion_tamizable))
espacio_muestral
## [1] 0 1 2 3 4
# Una brigada: selección aleatoria de 50 personas sin reemplazo
una_brigada <- sample(poblacion_tamizable, size = 50)
table(una_brigada)
## una_brigada
##  0  1  2 
## 18 16 16
sum(una_brigada >= 2)   # ocurrencias del evento A en esta brigada
## [1] 16
# 1.000 brigadas simuladas
brigadas <- replicate(1000, sum(sample(poblacion_tamizable, size = 50) >= 2))

tibble(`Promedio de casos` = mean(brigadas),
       `Mínimo` = min(brigadas), `Máximo` = max(brigadas),
       `Percentil 5` = quantile(brigadas, 0.05),
       `Percentil 95` = quantile(brigadas, 0.95))
p95_brigadas <- quantile(brigadas, 0.95)

g_brigadas <- tibble(casos = brigadas) |>
  count(casos) |>
  mutate(situacion = if_else(casos > p95_brigadas, "Supera la reserva (P95)", "Cubierto por la reserva"),
         porcentaje = round(100 * n / sum(n), 1)) |>
  ggplot(aes(x = casos, y = porcentaje, fill = situacion)) +
  geom_col(colour = "white") +
  geom_vline(xintercept = mean(brigadas), linetype = "dashed", colour = paleta[["gris"]]) +
  scale_fill_manual(values = c(paleta[["azul"]], paleta[["naranja"]])) +
  labs(title = "1.000 brigadas simuladas de 50 adultos",
       x = "Personas con 2 o más factores de riesgo", y = "% de brigadas", fill = NULL) +
  theme_minimal(base_size = 12)

ggplotly(g_brigadas) |> layout(legend = list(orientation = "h", y = -0.25))

Interpretación:

  • En promedio, cada brigada encuentra 13.7 personas con 2 o más factores de riesgo (27,9 % de la población). Ese número varía de una brigada a otra solo por azar, entre 3 y 26.
  • Si se reservan consultas solo para el promedio, cerca de la mitad de las brigadas quedan cortas.
  • Decisión: reservar 19 cupos de seguimiento por brigada (percentil 95) cubre 95 de cada 100 brigadas. Así, la reserva queda sustentada en datos en lugar de un número fijo arbitrario.

2.3.2 Frecuencias y probabilidades condicionales con table()

La probabilidad condicional, P(A | B), es la probabilidad de A sabiendo que ya ocurrió B. Por ejemplo: la probabilidad de mala salud sabiendo que el ingreso es bajo.

# Frecuencias absolutas y probabilidades empíricas
tabla_salud <- table(base_eda$salud_general)
tabla_salud
## 
## Excelente Muy buena     Buena   Regular      Mala 
##     81660    138139    119502     46239     15457
round(prop.table(tabla_salud) * 100, 1)
## 
## Excelente Muy buena     Buena   Regular      Mala 
##      20.4      34.4      29.8      11.5       3.9
# Tabla de doble entrada: ingreso x mala salud (1 = regular/mala)
tabla_ingreso_salud <- table(Ingreso = base_eda$ingreso, `Mala salud` = base_eda$mala_salud)
tabla_ingreso_salud
##               Mala salud
## Ingreso             0      1
##   < 15 mil      16181  10301
##   15-25 mil     35279  13360
##   25-35 mil     25170   6178
##   35-50 mil     37450   6329
##   50 mil o más 158963  12154
round(prop.table(tabla_ingreso_salud, margin = 1) * 100, 1)   # % por fila = P(mala salud | ingreso)
##               Mala salud
## Ingreso           0    1
##   < 15 mil     61.1 38.9
##   15-25 mil    72.5 27.5
##   25-35 mil    80.3 19.7
##   35-50 mil    85.5 14.5
##   50 mil o más 92.9  7.1
p_mala        <- mean(base_eda$mala_salud, na.rm = TRUE)
p_mala_bajo   <- prop.table(tabla_ingreso_salud, 1)["< 15 mil", "1"]
p_mala_alto   <- prop.table(tabla_ingreso_salud, 1)["50 mil o más", "1"]

tibble(`P(mala salud)` = pct(p_mala),
       `P(mala salud | ingreso < 15 mil)` = pct(p_mala_bajo),
       `P(mala salud | ingreso ≥ 50 mil)` = pct(p_mala_alto),
       `Veces más en ingreso bajo vs. alto` = round(p_mala_bajo / p_mala_alto, 1))

Interpretación:

  • Si se elige un adulto al azar, la probabilidad de que reporte salud regular o mala es 15,4 %. Si se sabe que su ingreso es menor de 15 mil USD, esa probabilidad sube a 38,9 % (2,5 veces el promedio). Si su ingreso es de 50 mil o más, baja a 7,1 %.
  • Decisión: conocer el ingreso cambia la probabilidad de mala salud. Por eso es un criterio válido para priorizar la búsqueda activa de pacientes y la asignación de citas preventivas.

2.4 Distribuciones discretas: planear recursos ante eventos que se cuentan

Una distribución discreta describe variables que se cuentan (0, 1, 2, …): pacientes, ingresos, camas.

2.4.1 Binomial: jornada de valoración preanestésica

Caso práctico: una jornada de valoración preanestésica atiende a 30 pacientes. La probabilidad de obesidad se toma de la base (32 %). Cada paciente con obesidad exige recursos adicionales: manguito de presión grande, camilla con mayor capacidad, evaluación de vía aérea difícil y más tiempo de consulta.

La distribución binomial cuenta cuántos “éxitos” (aquí, pacientes con obesidad) hay en un número fijo de intentos (30 pacientes) con la misma probabilidad cada uno. Se simula de dos formas:

  • runif() genera un número aleatorio entre 0 y 1 para cada paciente. Si es menor que 0,32, el paciente “tiene obesidad”.
  • rbinom() hace lo mismo de forma directa.
n_pacientes <- 30
p_obesidad  <- mean(base_eda$obesidad, na.rm = TRUE)
n_jornadas  <- 1000

jornadas_runif  <- replicate(n_jornadas, sum(runif(n_pacientes) < p_obesidad))
jornadas_rbinom <- rbinom(n_jornadas, size = n_pacientes, prob = p_obesidad)

tibble(Método = c("Simulación con runif()", "Simulación con rbinom()", "Valor teórico"),
       `Promedio de pacientes con obesidad` = round(c(mean(jornadas_runif), mean(jornadas_rbinom),
                                                      n_pacientes * p_obesidad), 1),
       `P(15 o más)` = pct(c(mean(jornadas_runif >= 15), mean(jornadas_rbinom >= 15),
                            1 - pbinom(14, n_pacientes, p_obesidad))))
cupos_95 <- qbinom(0.95, n_pacientes, p_obesidad)   # cantidad que cubre el 95 % de las jornadas
cupos_95
## [1] 14
comparacion_binomial <- bind_rows(
  tibble(metodo = "Simulación con runif()",  obesos = jornadas_runif),
  tibble(metodo = "Simulación con rbinom()", obesos = jornadas_rbinom)) |>
  count(metodo, obesos) |>
  mutate(probabilidad = n / n_jornadas)

teorico_binomial <- tibble(obesos = 0:22, probabilidad = dbinom(0:22, n_pacientes, p_obesidad))

g_binomial <- ggplot(comparacion_binomial, aes(x = obesos, y = probabilidad, fill = metodo)) +
  geom_col(position = "dodge", alpha = 0.85) +
  geom_line(data = teorico_binomial, aes(x = obesos, y = probabilidad), inherit.aes = FALSE,
            colour = paleta[["naranja"]], linewidth = 1) +
  geom_vline(xintercept = cupos_95 + 0.5, linetype = "dashed", colour = paleta[["gris"]]) +
  scale_fill_manual(values = c(paleta[["verde"]], paleta[["azul"]])) +
  labs(title = "Pacientes con obesidad en una jornada de 30 valoraciones",
       subtitle = "Barras: 1.000 jornadas simuladas | Línea naranja: valor teórico",
       x = "Pacientes con obesidad", y = "Probabilidad", fill = NULL) +
  theme_minimal(base_size = 12)

ggplotly(g_binomial) |> layout(legend = list(orientation = "h", y = -0.25))

Interpretación:

  • Las dos simulaciones y la curva teórica coinciden: en promedio hay 9,6 pacientes con obesidad por jornada (casi 1 de cada 3). Una jornada puede tener desde 3-4 hasta 16-17 según el azar.
  • La probabilidad de recibir 15 o más en un día es de 3 %.
  • Decisión: preparar equipo e insumos para 14 pacientes con obesidad por jornada cubre el 95 % de los días. Además, se puede reservar más tiempo de consulta para ese volumen de pacientes desde la programación.

2.4.2 Poisson: escenarios de demanda de camas en la unidad coronaria

Caso práctico: un hospital de referencia necesita definir cuántas camas de unidad coronaria abrir. Para eso analiza el número de ingresos por síndrome coronario agudo por día en tres escenarios. Los promedios de 2, 4 y 6 ingresos diarios son supuestos ilustrativos; cada institución usa su histórico. La capacidad actual es de 5 camas.

La distribución de Poisson describe cuántos eventos ocurren en un intervalo de tiempo cuando llegan de forma independiente y con un promedio conocido (λ, lambda). dpois() da la probabilidad exacta de cada valor.

escenarios <- c("Baja demanda (media 2)" = 2, "Demanda habitual (media 4)" = 4, "Alta demanda (media 6)" = 6)
camas_disponibles <- 5

distribucion_poisson <- map_dfr(names(escenarios), function(nombre) {
  tibble(escenario = nombre, ingresos = 0:14,
         probabilidad = dpois(0:14, lambda = escenarios[[nombre]]))
}) |>
  mutate(escenario = factor(escenario, levels = names(escenarios)))

distribucion_poisson |>
  group_by(escenario) |>
  summarise(`P(ningún ingreso)` = pct(probabilidad[ingresos == 0]),
            `P(se superan las 5 camas)` = pct(1 - sum(probabilidad[ingresos <= camas_disponibles])),
            `Camas para cubrir el 95 % de los días` = min(ingresos[cumsum(probabilidad) >= 0.95]))
g_poisson <- ggplot(distribucion_poisson, aes(x = ingresos, y = probabilidad, fill = escenario)) +
  geom_col(position = "dodge") +
  geom_vline(xintercept = camas_disponibles + 0.5, linetype = "dashed", colour = paleta[["gris"]]) +
  annotate("text", x = camas_disponibles + 0.7, y = 0.26, label = "Capacidad: 5 camas",
           hjust = 0, colour = paleta[["gris"]]) +
  scale_fill_manual(values = c(paleta[["verde"]], paleta[["azul"]], paleta[["naranja"]])) +
  scale_x_continuous(breaks = 0:14) +
  labs(title = "Ingresos diarios a unidad coronaria en tres escenarios",
       x = "Ingresos por día", y = "Probabilidad", fill = NULL) +
  theme_minimal(base_size = 12)

ggplotly(g_poisson) |> layout(legend = list(orientation = "h", y = -0.25))

Interpretación:

  • Demanda baja: las 5 camas se superan solo el 1,7 % de los días (unos 6 días al año). La capacidad es suficiente.
  • Demanda habitual: se superan el 21,5 % de los días (1 de cada 5). Para cubrir el 95 % de los días se necesitan 8 camas.
  • Demanda alta: se superan más de la mitad de los días (55,4 %); se necesitarían 10 camas.
  • Decisión: la capacidad no se dimensiona con el promedio, sino con la probabilidad de superarla. Con estos cálculos se puede justificar ante la junta directiva:
    • camas flexibles o convenios de referencia para los días que superan la capacidad;
    • el costo de oportunidad de las camas vacías en los días de baja demanda: con media 2, el 13,5 % de los días no hay ningún ingreso.

2.5 Distribuciones continuas: mediciones y tiempos

Una distribución continua describe variables que pueden tomar cualquier valor dentro de un intervalo (kg/m², minutos). Para estas distribuciones se usan tres tipos de función:

  • Densidad (dnorm, dexp, dunif): la forma de la curva.
  • Distribución acumulada (pnorm, pexp): la probabilidad de quedar por debajo de un valor.
  • Generación aleatoria (rnorm, rexp, runif): crea datos simulados.

2.5.1 Normal: ¿sirve para estimar cuántos pacientes con obesidad severa llegarán?

Caso práctico: compras quiere estimar cuántos pacientes con obesidad grado III (IMC ≥ 40) se atenderán, para adquirir camillas y mesas quirúrgicas bariátricas. Un analista propone suponer que el IMC sigue una distribución normal (curva en forma de campana, simétrica). ¿Es correcto?

media_imc <- mean(base_eda$imc, na.rm = TRUE)
de_imc    <- sd(base_eda$imc, na.rm = TRUE)

imc_simulado <- rnorm(10000, mean = media_imc, sd = de_imc)   # 10.000 IMC simulados con la normal

tibble(Umbral = c("Bajo peso (< 18,5)", "Obesidad (≥ 30)", "Obesidad grado III (≥ 40)"),
       `Real (datos BRFSS)` = pct(c(mean(base_eda$imc < 18.5, na.rm = TRUE),
                                    mean(base_eda$imc >= 30, na.rm = TRUE),
                                    mean(base_eda$imc >= 40, na.rm = TRUE))),
       `Teórico normal (pnorm)` = pct(c(pnorm(18.5, media_imc, de_imc),
                                        1 - pnorm(30, media_imc, de_imc),
                                        1 - pnorm(40, media_imc, de_imc))),
       `Simulado (rnorm)` = pct(c(mean(imc_simulado < 18.5), mean(imc_simulado >= 30),
                                  mean(imc_simulado >= 40))))
curva_normal <- tibble(imc = seq(10, 60, by = 0.5),
                       densidad = dnorm(imc, mean = media_imc, sd = de_imc))

g_normal <- ggplot(filter(base_eda, imc <= 60), aes(x = imc)) +
  geom_histogram(aes(y = after_stat(density)), binwidth = 1,
                 fill = paleta[["azul"]], colour = "white", alpha = 0.8) +
  geom_line(data = curva_normal, aes(x = imc, y = densidad),
            colour = paleta[["naranja"]], linewidth = 1.2) +
  geom_vline(xintercept = c(30, 40), linetype = "dashed", colour = paleta[["gris"]]) +
  labs(title = "IMC real (barras azules) frente a la curva normal teórica (línea naranja)",
       x = "IMC (kg/m²)", y = "Densidad") +
  theme_minimal(base_size = 12)

ggplotly(g_normal)

Interpretación:

  • La curva normal no se ajusta bien. El IMC real tiene la cola derecha más larga que la campana.
  • La normal sobreestima la obesidad (39,5 % frente al 32,0 % real) y el bajo peso (6,2 % frente a 1,7 %). En cambio, subestima la obesidad grado III (3,3 % frente al 5,1 % real). Por cada 1.000 pacientes, se planearían 33 camillas bariátricas cuando la necesidad real es de 51: un déficit de casi 40 %.
  • Decisión: antes de usar la distribución normal para planear recursos, hay que compararla con los datos reales. Para extremos clínicamente críticos se debe usar la proporción observada.

2.5.2 Exponencial y uniforme: tiempos de espera y de recambio de quirófano

Caso práctico A, tiempo de espera en consulta prioritaria (exponencial): el tablero de indicadores reporta una espera promedio de 30 minutos (supuesto ilustrativo). La distribución exponencial describe tiempos de espera hasta que ocurre un evento: muchas esperas cortas y pocas muy largas.

Caso práctico B, recambio de quirófano (uniforme): el tiempo entre la salida de un paciente y la entrada del siguiente (limpieza y alistamiento) toma entre 20 y 40 minutos, con cualquier valor igualmente probable. Eso es una distribución uniforme.

espera_media <- 30
esperas   <- rexp(1000, rate = 1 / espera_media)    # 1.000 pacientes simulados
recambios <- runif(1000, min = 20, max = 40)        # 1.000 recambios simulados

tibble(Indicador = c("P(espera ≤ 30 min)", "P(espera > 60 min)", "Mediana de espera (min)",
                     "P(recambio > 35 min)", "Promedio de recambio (min)"),
       Teórico  = round(c(pexp(30, 1 / espera_media), 1 - pexp(60, 1 / espera_media),
                          qexp(0.5, 1 / espera_media), (40 - 35) / (40 - 20), 30), 3),
       Simulado = round(c(mean(esperas <= 30), mean(esperas > 60), median(esperas),
                          mean(recambios > 35), mean(recambios)), 3))
# Sala con 8 cirugías = 7 recambios: tiempo total perdido en recambios por día
recambio_diario <- replicate(1000, sum(runif(7, min = 20, max = 40)))
datos_tiempos <- bind_rows(
  tibble(proceso = "A. Espera en consulta prioritaria (exponencial)", minutos = esperas),
  tibble(proceso = "B. Recambio de quirófano (uniforme)", minutos = recambios))

curvas_tiempos <- bind_rows(
  tibble(proceso = "A. Espera en consulta prioritaria (exponencial)",
         minutos = seq(0, 200, 1), densidad = dexp(seq(0, 200, 1), rate = 1 / espera_media)),
  tibble(proceso = "B. Recambio de quirófano (uniforme)",
         minutos = seq(15, 45, 0.5), densidad = dunif(seq(15, 45, 0.5), min = 20, max = 40)))

g_tiempos <- ggplot(datos_tiempos, aes(x = minutos)) +
  geom_histogram(aes(y = after_stat(density), fill = proceso), bins = 30, colour = "white", alpha = 0.8) +
  geom_line(data = curvas_tiempos, aes(x = minutos, y = densidad), colour = paleta[["naranja"]], linewidth = 1) +
  facet_wrap(~ proceso, scales = "free") +
  scale_fill_manual(values = c(paleta[["azul"]], paleta[["verde"]])) +
  labs(title = "Tiempos operativos simulados (barras) y curva teórica (línea naranja)",
       x = "Minutos", y = "Densidad") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "none")

ggplotly(g_tiempos)

Interpretación:

  • Espera exponencial: con un promedio de 30 minutos, la mitad de los pacientes espera menos de 21 minutos. Sin embargo, el 13,5 % espera más de 1 hora (pexp). El promedio del tablero oculta a ese grupo, que es el que genera quejas y riesgo clínico.
    • Decisión: el indicador de oportunidad debe reportar el porcentaje de pacientes que superan el estándar, no solo el promedio.
  • Recambio uniforme: el 25 % de los recambios supera 35 minutos. En una sala con 8 cirugías, los 7 recambios suman en promedio 210 minutos al día (unas 3,5 horas). Llegan a superar 235 minutos en el 5 % de los días.
    • Decisión: estandarizar el alistamiento (listas de chequeo, equipo dedicado) para bajar el recambio promedio de 30 a 22 minutos libera cerca de una hora de quirófano por sala al día (7 recambios × 8 min), equivalente a una cirugía corta adicional.

2.6 Estadística inferencial: ¿las diferencias son reales?

La estadística inferencial evalúa si una diferencia observada puede deberse al azar. Conceptos clave:

  • Hipótesis nula (H0): supone que no hay diferencia entre los grupos.
  • Valor p: probabilidad de observar una diferencia como la encontrada (o mayor) si H0 fuera cierta. Si p < 0,05, se rechaza H0: la diferencia es estadísticamente significativa.
  • Intervalo de confianza del 95 % (IC 95 %): rango de valores donde, con 95 % de confianza, está el valor real.
  • Advertencia: con 400.000 registros casi cualquier diferencia resulta significativa. Por eso, en cada prueba se evalúa también si la diferencia es grande o relevante para la gestión, no solo el valor p.

2.6.1 Prueba t: ¿la barrera de costo se asocia con más días de mala salud?

  • Variable resultado: días de mala salud física en el último mes.
  • Grupos: personas que dejaron de consultar por costo vs. personas que no.
  • La prueba t compara las medias de dos grupos.
prueba_t <- t.test(dias_mala_salud_fisica ~ barrera_costo, data = base_eda)
prueba_t
## 
##  Welch Two Sample t-test
## 
## data:  dias_mala_salud_fisica by barrera_costo
## t = 59.114, df = 37218, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Sí and group No is not equal to 0
## 95 percent confidence interval:
##  3.332561 3.561132
## sample estimates:
## mean in group Sí mean in group No 
##         6.602538         3.155691
# Intervalo de confianza de una media: ¿los adultos duermen en promedio las 7 h recomendadas?
prueba_sueno <- t.test(base_eda$horas_sueno, mu = 7)
prueba_sueno
## 
##  One Sample t-test
## 
## data:  base_eda$horas_sueno
## t = 41.688, df = 397263, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 7
## 95 percent confidence interval:
##  7.092534 7.101664
## sample estimates:
## mean of x 
##  7.097099
base_eda |>
  filter(!is.na(barrera_costo), !is.na(dias_mala_salud_fisica)) |>
  ggbarplot(x = "barrera_costo", y = "dias_mala_salud_fisica", add = "mean_ci",
            fill = "barrera_costo", palette = c(paleta[["naranja"]], paleta[["azul"]]),
            xlab = "¿Dejó de consultar por costo?", ylab = "Días con mala salud física (media e IC 95 %)",
            title = "Barrera de costo y carga de enfermedad", legend = "none") +
  stat_compare_means(method = "t.test", label.x = 1.3, label.y = 7.5)

Interpretación:

  • Quienes dejaron de consultar por costo reportan 6,6 días de mala salud física al mes, frente a 3,2 días de quienes no tuvieron esa barrera. La diferencia es de 3,4 días (IC 95 %: 3,3 a 3,6), con p < 0,001: la diferencia no se debe al azar y es relevante, porque es el doble de días.
  • Prueba de sueño: la media es 7,10 h (IC 95 %: 7,09 a 7,10). Es estadísticamente distinta de 7 h (p < 0,001), pero la diferencia es de 6 minutos: no tiene relevancia práctica. Es el ejemplo de por qué no basta con el valor p.
  • Decisión: reducir copagos y cuotas moderadoras en población de riesgo no es solo un tema financiero. Esa población tiene el doble de carga de enfermedad, que luego se traduce en consultas, urgencias y ausentismo.

2.6.2 ANOVA: ¿los días de mala salud difieren según el nivel de ingreso?

  • Variable resultado: días de mala salud física.
  • Grupos: cinco niveles de ingreso.
  • El ANOVA (análisis de varianza) compara las medias de tres o más grupos a la vez.
  • La prueba de Tukey indica después cuáles pares de grupos difieren.
datos_anova <- base_eda |>
  filter(!is.na(ingreso), !is.na(dias_mala_salud_fisica)) |>
  mutate(ingreso = factor(ingreso, ordered = FALSE))

anova_ingreso <- aov(dias_mala_salud_fisica ~ ingreso, data = datos_anova)
summary(anova_ingreso)
##                 Df   Sum Sq Mean Sq F value Pr(>F)    
## ingreso          4  1023788  255947    4158 <2e-16 ***
## Residuals   316308 19472465      62                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Tamaño del efecto (eta cuadrado): % de la variación de los días que se explica por el ingreso
suma_cuadrados <- summary(anova_ingreso)[[1]][["Sum Sq"]]
eta_cuadrado <- suma_cuadrados[1] / sum(suma_cuadrados)
pct(eta_cuadrado)
## [1] "5 %"
TukeyHSD(anova_ingreso)   # comparaciones por pares
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = dias_mala_salud_fisica ~ ingreso, data = datos_anova)
## 
## $ingreso
##                              diff       lwr       upr p adj
## 15-25 mil-< 15 mil     -2.4076755 -2.573761 -2.241590     0
## 25-35 mil-< 15 mil     -3.7506770 -3.931971 -3.569383     0
## 35-50 mil-< 15 mil     -4.5509295 -4.719867 -4.381992     0
## 50 mil o más-< 15 mil  -5.7901065 -5.933718 -5.646495     0
## 25-35 mil-15-25 mil    -1.3430016 -1.499867 -1.186136     0
## 35-50 mil-15-25 mil    -2.1432541 -2.285659 -2.000849     0
## 50 mil o más-15-25 mil -3.3824310 -3.493622 -3.271240     0
## 35-50 mil-25-35 mil    -0.8002525 -0.960135 -0.640370     0
## 50 mil o más-25-35 mil -2.0394295 -2.172270 -1.906589     0
## 50 mil o más-35-50 mil -1.2391770 -1.354585 -1.123769     0
ggbarplot(datos_anova, x = "ingreso", y = "dias_mala_salud_fisica", add = "mean_ci",
          fill = "ingreso",
          palette = colorRampPalette(c(paleta[["naranja"]], paleta[["gris"]], paleta[["azul"]]))(5),
          xlab = "Ingreso anual del hogar (USD)", ylab = "Días con mala salud física (media e IC 95 %)",
          title = "Gradiente de carga de enfermedad por ingreso", legend = "none") +
  stat_compare_means(method = "anova", label.x = 3, label.y = 9)

Interpretación:

  • El ANOVA es significativo (F ≈ 4.158; p < 0,001): al menos un nivel de ingreso difiere. La prueba de Tukey muestra que todos los pares difieren entre sí, en un gradiente de 7,9 a 2,1 días.
  • El eta cuadrado (~5 %) indica que el ingreso explica una parte modesta de la variación: influye, pero no es la única causa. Coherente con el capítulo 1, la carga de enfermedad es multifactorial.
  • Decisión: el ingreso sirve como criterio de priorización, pero los programas deben combinarlo con factores clínicos y de estilo de vida.

2.6.3 Chi-cuadrado: ¿los factores de riesgo se asocian con la enfermedad coronaria?

La prueba chi-cuadrado evalúa si dos variables categóricas están asociadas, es decir, si conocer una cambia la distribución de la otra.

  • Residuos: muestran qué celdas tienen más (+) o menos (−) casos de los esperados si no hubiera asociación.
  • V de Cramér: mide la fuerza de la asociación de 0 a 1. Menos de 0,1 es débil; 0,3 o más, fuerte.
datos_chi <- base_eda |>
  filter(!is.na(n_factores_riesgo), !is.na(enf_coronaria)) |>
  mutate(grupo_factores = factor(pmin(n_factores_riesgo, 3), levels = 0:3,
                                 labels = c("0", "1", "2", "3 o más")))

tabla_chi <- table(`Factores de riesgo` = datos_chi$grupo_factores,
                   `Enfermedad coronaria` = datos_chi$enf_coronaria)
tabla_chi
##                   Enfermedad coronaria
## Factores de riesgo     Sí     No
##            0         8005 117614
##            1        10736 113608
##            2         7945  62852
##            3 o más   3734  21676
prueba_chi <- chisq.test(tabla_chi)
prueba_chi
## 
##  Pearson's Chi-squared test
## 
## data:  tabla_chi
## X-squared = 2547.7, df = 3, p-value < 2.2e-16
round(prueba_chi$residuals, 1)                                   # residuos de Pearson
##                   Enfermedad coronaria
## Factores de riesgo    Sí    No
##            0       -28.9   9.0
##            1        -1.8   0.6
##            2        21.9  -6.8
##            3 o más  31.8  -9.9
v_cramer <- sqrt(prueba_chi$statistic / (sum(tabla_chi) * (min(dim(tabla_chi)) - 1)))
round(unname(v_cramer), 3)
## [1] 0.086
datos_chi |>
  group_by(grupo_factores) |>
  summarise(prevalencia = round(100 * mean(enf_coronaria == "Sí"), 1)) |>
  ggbarplot(x = "grupo_factores", y = "prevalencia", fill = "grupo_factores",
            palette = c(paleta[["verde"]], paleta[["gris"]], paleta[["azul"]], paleta[["naranja"]]),
            label = TRUE, lab.vjust = -0.4,
            xlab = "Número de factores de riesgo modificables", ylab = "% con enfermedad coronaria",
            title = "Prevalencia de enfermedad coronaria según factores de riesgo", legend = "none")

Interpretación:

  • La asociación es significativa (χ² ≈ 2.548; p < 0,001). La prevalencia coronaria sube con cada factor: 6,4 % → 8,6 % → 11,2 % → 14,7 %. Las personas con 3 o más factores tienen 2,3 veces la prevalencia de quienes no tienen ninguno.
  • Los residuos muestran menos casos de lo esperado sin factores (−28,9) y más de lo esperado con 3 o más (+31,8).
  • La V de Cramér (~0,09) indica una asociación débil. Esto es esperable, porque la edad y la genética también pesan. Aun así, es una asociación consistente y basada en factores modificables.
  • Decisión: respalda el uso del conteo de factores de riesgo para estratificar a la población en el primer nivel de atención.

2.6.4 Prueba de proporciones: ¿tener seguro aumenta el chequeo anual?

Caso práctico: se compara la proporción de adultos con chequeo médico en el último año entre asegurados y no asegurados. Es el equivalente a comparar dos clínicas o dos modelos de contratación. prop.test() compara dos proporciones y da el intervalo de confianza de su diferencia.

resumen_chequeo <- base_eda |>
  filter(!is.na(tiene_seguro), !is.na(ultimo_chequeo)) |>
  group_by(tiene_seguro) |>
  summarise(con_chequeo = sum(ultimo_chequeo == "< 1 año"), total = n()) |>
  mutate(porcentaje = round(100 * con_chequeo / total, 1))
resumen_chequeo
prueba_proporciones <- prop.test(x = resumen_chequeo$con_chequeo, n = resumen_chequeo$total)
prueba_proporciones
## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  resumen_chequeo$con_chequeo out of resumen_chequeo$total
## X-squared = 18186, df = 1, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.3126782 0.3237708
## sample estimates:
##    prop 1    prop 2 
## 0.8110629 0.4928384
ggbarplot(resumen_chequeo, x = "tiene_seguro", y = "porcentaje", fill = "tiene_seguro",
          palette = c(paleta[["azul"]], paleta[["naranja"]]),
          label = paste0(resumen_chequeo$porcentaje, " %"), lab.vjust = -0.4,
          xlab = "¿Tiene seguro de salud?", ylab = "% con chequeo en el último año",
          title = "Aseguramiento y uso de servicios preventivos", legend = "none") +
  ylim(0, 100)

Interpretación:

  • El 81,1 % de los asegurados tuvo chequeo en el último año, frente al 49,3 % de los no asegurados: una diferencia de 31,8 puntos porcentuales (IC 95 %: 31,3 a 32,4), con p < 0,001.
  • Sin seguro, la probabilidad de un chequeo anual cae casi a la mitad. Esta es la diferencia más grande del capítulo y conecta directamente lo financiero con la prevención.
  • Decisión: la afiliación y la cobertura son la puerta de entrada a la prevención. Las estrategias de búsqueda activa deben llevar los servicios preventivos a donde está la población no asegurada, en lugar de esperar a que consulte.

2.7 Síntesis de hallazgos para la toma de decisiones

# Hallazgo estadístico Decisión de gestión
1 Los días de mala salud tienen un CV > 200 %: el promedio no representa al paciente típico. Reportar mediana y percentiles; gestionar aparte a los usuarios de alta carga.
2 P(mala salud) pasa de 15,4 % en general a 38,9 % si el ingreso es bajo. Usar el ingreso como criterio de priorización en búsqueda activa.
3 Las simulaciones (brigadas, binomial, Poisson) muestran que la demanda varía por azar alrededor del promedio. Dimensionar cupos, insumos y camas con el percentil 95, no con el promedio.
4 La normal subestima la obesidad grado III (3,3 % vs. 5,1 % real); con espera promedio de 30 min, el 13,5 % espera más de 1 h. Validar los supuestos estadísticos antes de planear; medir el % de pacientes fuera del estándar.
5 Barrera de costo: el doble de días de mala salud (diferencia de 3,4 días; p < 0,001). Revisar copagos en población de riesgo.
6 El ingreso explica ~5 % de la variación (ANOVA) y la asociación factores-coronaria es débil pero consistente (V ≈ 0,09). Combinar criterios socioeconómicos y clínicos; ninguna variable sola basta.
7 Sin seguro, el chequeo anual cae de 81,1 % a 49,3 % (−31,8 puntos). Priorizar afiliación y servicios preventivos extramurales.