Caracterización de perfiles de alto riesgo asociados al desempleo en la República de Panamá, mediante Regresión Logística: años 2019, 2021 y 2024

Librerias de R utilizadas para el procesamiento, modelado y visualización:

pacman

foreign

ggplot2

dplyr

pROC

caret

---

A continuación, el código:

# Limpieza de entorno
rm(list=ls())

# 1. Cargar librerías
pacman::p_load(foreign, readr, dplyr, ggplot2, car, pROC, caret)

# Cargar bases .dbf
base2024 <- read.dbf("C:/Users/Oliver/OneDrive/Universidad/unknown/Data bruta/Extraida/Data tesis/EML_OCTUBRE_2024/persona.dbf", as.is = TRUE)

# Agregar columna para ID de cada individuo
base2024 <- base2024 %>% mutate(ID_individuo = paste0("ENC_", sprintf("%05d", row_number())))

# Filtrar población entre 15 y 80 años
base2024_fil <- base2024 %>% filter(!is.na(P3) & P3 >= 15 & P3 <= 80)

# Eliminar los individuos "NEA"
base2024_fil2 <- base2024_fil %>% filter(PEA_NEA != "NEA")

# 2. Preparación y Limpieza
datos2024 <- base2024_fil2 %>%
  mutate(ID_individuo = paste0("ENC_2024_", row_number())) %>%
  select(ID_individuo, OCU_DES, PROVINCIA, P1, P2, P3, P3_RECO, P4D_INDIGE, 
         P4F_AFROD, P5_CONYUGA, P7, GRADO_AP, AREARECO) %>% # Nota: P4 (Seguro) quitada del select para evitar tentación
  mutate(
    # Variable dependiente (1=Desempleado)
    desempleado = as.factor(ifelse(OCU_DES == "Desocupados", 1, 0)),
    
    # Región
    REGION = factor(case_when(
      PROVINCIA %in% c("08","13","03") ~ "Metropolitana",
      PROVINCIA %in% c("04","01","12") ~ "Occidental",
      PROVINCIA %in% c("02","06","07","09") ~ "Central",
      PROVINCIA %in% c("05","10","11") ~ "Oriental"
    )),
    
    # Variables sociodemográficas
    jefe_hogar = factor(ifelse(P1 == 1, 1, 0)),
    sexo = factor(ifelse(P2 == 1, 0, 1)), # 0=Hombre, 1=Mujer
    identidad_indigena = factor(ifelse(P4D_INDIGE == "11", 0, 1)),
    identidad_afro = factor(ifelse(P4F_AFROD == 8, 0, 1)),
    asiste_escuela = factor(ifelse(P7 == 1, 1, 0)),
    area = factor(ifelse(AREARECO %in% c("R","Rural"), 1, 0)),
    
    #Nivel educativo
    nivel_edu = factor(case_when(
      GRADO_AP == "NG" ~ "Ningún Grado",
      GRADO_AP %in% c("P1A3","P4A6","S1A3") ~ "Educación Básica General",
      GRADO_AP %in% c("S4A6","VOC") ~ "Educación Media",
      GRADO_AP %in% c("No U","UNIV") ~ "Educación Superior"
    )),
    
    # Estado Conyugal
    estado_conyugal = factor(case_when(
      P5_CONYUGA == 1 ~ "Unido", P5_CONYUGA == 4 ~ "Casado",
      P5_CONYUGA %in% c(2,3,5) ~ "SeparadoDivorciado",
      P5_CONYUGA == 6 ~ "Viudo", P5_CONYUGA == 7 ~ "Soltero"
    ))
  ) %>%
  na.omit()

# ==============================================================================
# 2.1. DEFINICIÓN EXPLÍCITA DE CATEGORÍAS DE REFERENCIA (BASES)
# ==============================================================================
# Esto asegura que la interpretación sea siempre "Riesgo respecto a la base"

# Educación: Base = Universitaria (Para ver el riesgo de tener menos estudios)
datos2024$nivel_edu <- relevel(datos2024$nivel_edu, ref = "Educación Superior")

# Región: Base = Metropolitana (Centro económico)
datos2024$REGION <- relevel(datos2024$REGION, ref = "Metropolitana")

# Sexo: Base = 0 (Hombre). Resultado mostrará efecto de ser Mujer.
datos2024$sexo <- relevel(datos2024$sexo, ref = "0")

# Jefe de Hogar: Base = 1 (Jefe). Resultado mostrará efecto de NO ser jefe (Hijo/Dependiente).
datos2024$jefe_hogar <- relevel(datos2024$jefe_hogar, ref = "1")

# Etnia: Base = 0 (No perteneciente). Resultado mostrará efecto de pertenecer.
datos2024$identidad_indigena <- relevel(datos2024$identidad_indigena, ref = "0")
datos2024$identidad_afro <- relevel(datos2024$identidad_afro, ref = "0")

# Área: Base = 0 (Urbana). Resultado mostrará efecto Rural.
datos2024$area <- relevel(datos2024$area, ref = "0")

# Estado Civil: Base = Casado (Estabilidad). Resultado mostrará efecto Soltero/Unido.
# Nota: Asegúrate de que "Casado" esté escrito exactamente igual que en tus datos.
datos2024$estado_conyugal <- relevel(datos2024$estado_conyugal, ref = "Casado")

# Asistencia Escolar: Base = 0 (No asiste). Resultado mostrará efecto de Asistir.
datos2024$asiste_escuela <- relevel(datos2024$asiste_escuela, ref = "0")

# 3. Modelo de Regresión Logística
# Crear variable centrada
media_edad_2024 <- mean(datos2024$P3, na.rm = TRUE) # Guardamos la media para usarla luego en el gráfico
datos2024 <- datos2024 %>% mutate(P3c = P3 - media_edad_2024)

# Modelo FINAL
modelo_2024 <- glm(desempleado ~ REGION + jefe_hogar + sexo + P3c + I(P3c^2) + 
                     identidad_indigena + identidad_afro + estado_conyugal + 
                     asiste_escuela + nivel_edu + area, data = datos2024, family = binomial())
summary(modelo_2024)

# ==============================================================================
# TRADUCCIÓN Y ESTÉTICA DE VARIABLES (DICCIONARIO)
# ==============================================================================
# Creamos un vector con los nombres
# Izquierda: Nombre que sale en R  ->  Derecha: Nombre para tu Tesis
nombres_bonitos <- c(
  "(Intercept)"                       = "(Intercepto)",
  "REGIONCentral"                     = "Región Central",
  "REGIONOccidental"                  = "Región Occidental",
  "REGIONOriental"                    = "Región Oriental",
  "jefe_hogar0"                       = "No Jefe de Hogar",
  "sexo1"                             = "Sexo Femenino (Vs. Masculino)",
  "P3c"                               = "Edad Centrada",
  "I(P3c^2)"                          = "Edad Centrada^2 (Efecto Curva)",
  "identidad_indigena1"               = "Identidad Indígena",
  "identidad_afro1"                   = "Identidad Afrodescendiente",
  "estado_conyugalSeparadoDivorciado" = "Estado Conyugal: Separado/Divorciado",
  "estado_conyugalSoltero"            = "Estado Conyugal: Soltero",
  "estado_conyugalUnido"              = "Estado Conyugal: Unido",
  "estado_conyugalViudo"              = "Estado Conyugal: Viudo",
  "asiste_escuela1"                   = "Asiste a Escuela",
  "nivel_eduNingún Grado"              = "Educación: Ningún Grado",
  "nivel_eduEducación Superior"          = "Educación: Superior",
  "nivel_eduEducación Básica General"                 = "Educación: Básica General",
  "nivel_eduEducación Media"               = "Educación: Media",
  "area1"                             = "Área Rural"
)

# ==============================================================================
# 4. VALIDACIÓN Y JUSTIFICACIÓN DEL UMBRAL (0.5 vs YOUDEN)
# ==============================================================================
# A. Generar predicciones (Probabilidades)
prob_2024 <- predict(modelo_2024, type = "response")

# --- ESCENARIO 1: UMBRAL POR DEFECTO (0.5) ---
# Este es el estándar, pero suele fallar en datos desbalanceados como el desempleo
pred_default <- factor(ifelse(prob_2024 > 0.5, 1, 0), levels = c(0, 1))
cm_default <- confusionMatrix(pred_default, datos2024$desempleado, positive = "1")

# --- ESCENARIO 2: UMBRAL OPTIMIZADO (ÍNDICE DE YOUDEN) ---
# Busca el punto de corte que maximiza la Sensibilidad + Especificidad
roc_2024 <- roc(datos2024$desempleado, prob_2024, quiet = TRUE)
umbral_youden <- coords(roc_2024, "best", ret = "threshold", best.method = "youden")[[1]]

pred_youden <- factor(ifelse(prob_2024 > umbral_youden, 1, 0), levels = c(0, 1))
cm_youden <- confusionMatrix(pred_youden, datos2024$desempleado, positive = "1")

# --- COMPARACIÓN Y JUSTIFICACIÓN ---
# Creamos una tabla para ver la mejora claramente
comparacion_umbrales <- data.frame(
  Metrica = c("Umbral Utilizado", "Accuracy (Global)", "Sensibilidad (Detectar Desempleo)", "Especificidad (Detectar Empleo)", "Balanced Accuracy", "Índice Kappa"),
  Modelo_Estandar_0.5 = c(
    0.5,
    round(cm_default$overall["Accuracy"], 3),
    round(cm_default$byClass["Sensitivity"], 3),
    round(cm_default$byClass["Specificity"], 3),
    round(cm_default$byClass["Balanced Accuracy"], 3),
    round(cm_default$overall["Kappa"], 3)
  ),
  Modelo_Ajustado_Youden = c(
    round(umbral_youden, 3),
    round(cm_youden$overall["Accuracy"], 3),
    round(cm_youden$byClass["Sensitivity"], 3),
    round(cm_youden$byClass["Specificity"], 3),
    round(cm_youden$byClass["Balanced Accuracy"], 3),
    round(cm_youden$overall["Kappa"], 3)
  )
)

print("--- JUSTIFICACIÓN TÉCNICA DEL AJUSTE DE UMBRAL ---")
print(comparacion_umbrales)

# Matriz de confusion con el umbral estandard
print("--- Matriz de Confusión con Umbral Estándar (0.5) ---")
print(cm_default)

# --- VISUALIZACIÓN DE RESULTADOS ---
# A. Tabla de Odds Ratios
tabla_OR_2024 <- exp(cbind(OR = coef(modelo_2024), confint(modelo_2024)))
print("--- Odds Ratios 2024 ---")
print(round(tabla_OR_2024, 3))

# B. Matriz de Confusión
prediccion_clasificada <- factor(ifelse(prob_2024 > umbral_youden, 1, 0), levels = c(0, 1))
conf_matrix <- confusionMatrix(prediccion_clasificada, datos2024$desempleado, positive = "1")
print(conf_matrix)

# ==============================================================================
# C. FOREST PLOT ESTÉTICO
# ==============================================================================

# 1. Extraemos los datos y aplicamos los nombres bonitos
tabla_OR_plot <- as.data.frame(exp(cbind(OR = coef(modelo_2024), confint(modelo_2024)))) %>%
  mutate(
    Variable_Cruda = rownames(.),
    Variable_Bonita = nombres_bonitos[Variable_Cruda],
    inf = `2.5 %`,
    sup = `97.5 %`
  ) %>%
  # Si alguna variable no estaba en el diccionario, usa la original para no romper nada
  mutate(Variable_Bonita = ifelse(is.na(Variable_Bonita), Variable_Cruda, Variable_Bonita)) %>%
  filter(!Variable_Cruda %in% c("(Intercept)", "P3c", "I(P3c^2)"))

# 2. Definir colores
c_riesgo <- "#b03a2e" # Rojo
c_protec <- "#2b4c7e" # Azul
c_neutro <- "#bdc3c7" # Gris

# 3. Clasificar para colorear
tabla_OR_plot <- tabla_OR_plot %>%
  mutate(
    Tipo = case_when(
      OR > 1 & inf > 1 ~ "Riesgo Aumentado",    # Significativo y mayor a 1
      OR < 1 & sup < 1 ~ "Factor Protector",    # Significativo y menor a 1
      TRUE ~ "No Significativo"                 # Cruza el 1
    )
  )

# 4. Graficar
ggplot(tabla_OR_plot, aes(x = reorder(Variable_Bonita, OR), y = OR, color = Tipo)) +
  geom_hline(yintercept = 1, linetype = "dashed", color = "black") +
  geom_errorbar(aes(ymin = inf, ymax = sup), width = 0.2, size = 0.8) +
  geom_point(size = 3.5) +
  coord_flip() +
  scale_y_log10() +
  scale_color_manual(values = c("Riesgo Aumentado" = c_riesgo, 
                                "Factor Protector" = c_protec, 
                                "No Significativo" = c_neutro)) +
  labs(title = "Determinantes del Desempleo en Panamá (2024)",
       x = "", 
       y = "Odds Ratio (Escala Logarítmica)",
       caption = "Fuente: Elaboración propia con datos EML 2024.") +
  theme_minimal() +
  theme(
    text = element_text(family = "sans", size = 12),
    plot.title = element_text(face = "bold", size = 14),
    axis.text.y = element_text(size = 11, color = "#1f1f20"),
    legend.position = "bottom"
  )

# ==============================================================================
# C. FOREST PLOT CON ENCABEZADO PARA LA COLUMNA DE DATOS
# ==============================================================================
# 1. Preparación de datos (Igual que antes, pero asegurando el texto)
tabla_OR_plot <- tabla_OR_plot %>%
  mutate(
    Etiqueta_Texto = sprintf("%.2f (%.2f - %.2f)", OR, inf, sup),
    Tipo = case_when(
      OR > 1 & inf > 1 ~ "Riesgo Aumentado",
      OR < 1 & sup < 1 ~ "Factor Protector",
      TRUE ~ "No Significativo"
    )
  )

# Determinamos la posición X (vertical en el plot) para el encabezado
# Será el número de variables + 1
posicion_encabezado <- nrow(tabla_OR_plot) + 1
max_or_eje <- max(tabla_OR_plot$sup) * 2 # Espacio para que quepa el texto largo

# 2. Graficar
ggplot(tabla_OR_plot, aes(x = reorder(Variable_Bonita, OR), y = OR)) +
  # Línea de referencia
  # geom_hline(yintercept = 1, linetype = "dashed", color = "gray60") +
  geom_hline(yintercept = 1, linetype = "dashed", color = "red", size = 1)+

# Barras y puntos
geom_errorbar(aes(ymin = inf, ymax = sup, color = Tipo), width = 0.3, size = 0.8) +
  geom_point(aes(color = Tipo), size = 3.5) +
  
  # VALORES A LA DERECHA
  geom_text(aes(label = Etiqueta_Texto, y = max_or_eje),
            hjust = 1, size = 3.5, color = "black") +
  
  # --- EL ENCABEZADO DE LA COLUMNA ---
  annotate("text", x = posicion_encabezado, y = max_or_eje,
           label = "OR (IC 95%)", fontface = "bold", size = 4, hjust = 1) +
  # ----------------------------------

coord_flip(clip = "off") + # 'clip=off' evita que el encabezado se corte arriba
  scale_y_log10(limits = c(min(tabla_OR_plot$inf)*0.7, max_or_eje)) +
  scale_color_manual(values = c("Riesgo Aumentado" = "#b03a2e",
                                "Factor Protector" = "#2b4c7e",
                                "No Significativo" = "#bdc3c7")) +
  labs(
    title = "Determinantes del Desempleo en Panamá (2024)",
    subtitle = "Análisis de Odds Ratios mediante Regresión Logística",
    x = "",
    y = "Probabilidad Relativa (Escala Logarítmica)",
    caption = "Fuente: Elaboración propia con datos EML 2024."
  ) +
  theme_minimal() +
  theme(
    text = element_text(size = 12),
    plot.title = element_text(face = "bold", size = 14),
    axis.text.y = element_text(size = 10, face = "bold", color = "black"),
    plot.margin = margin(t = 20, r = 10, b = 10, l = 10),
    legend.position = "bottom",
    legend.title = element_blank(),
    panel.grid.minor = element_blank(), # ya lo tienes
    # 👇 Personalización de la cuadrícula
    panel.grid.major = element_line(color = "#AFE0D9", size = 0.5),
    # Fondo transparente
    panel.background = element_rect(fill = NA, color = NA),
    plot.background  = element_rect(fill = NA, color = NA)
  )
################################################################################
# D. Gráfico de la Edad en 'U' (CORREGIDO LÓGICA DE CENTRADO)
################################################################################
# 1. Creamos un rango de edades REALES (15 a 80)
rango_edad_real <- seq(15, 80, 1)

# 2. Creamos el dataframe de predicción
perfil_ejemplo <- data.frame(
  # Restamos la media guardada arriba para que el modelo entienda el valor
  P3c = rango_edad_real - media_edad_2024, 
  
  # Variables fijas para el perfil "promedio" o de interés
  jefe_hogar = "1", 
  REGION = "Metropolitana", 
  sexo = "0", # Hombre
  identidad_indigena = "0", 
  identidad_afro = "0", 
  estado_conyugal = "Soltero", 
  asiste_escuela = "0", 
  nivel_edu = "Educación Superior", 
  area = "0"
)

# 3. Predecimos
perfil_ejemplo$prob <- predict(modelo_2024, newdata = perfil_ejemplo, type = "response")

# 4. Agregamos la edad REAL al dataframe para el eje X del gráfico
perfil_ejemplo$Edad_Real <- rango_edad_real

# 5. Graficamos usando Edad_Real en X
ggplot(perfil_ejemplo, aes(x = Edad_Real, y = prob)) +
  geom_line(color = "firebrick", size = 1.2) +
  labs(title = "Probabilidad de Desempleo según Edad (2024)", 
       subtitle = "Perfil ajustado: Hombre, Edu. Superior, Metro, Soltero",
       x = "Edad (Años)", y = "Probabilidad Estimada") +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold"))

# E. Diagnóstico VIF
print("--- Verificación de Multicolinealidad (VIF) ---")
print(vif(modelo_2024))

# ==============================================================================
# F. IDENTIFICACIÓN DEL TOP 5 DE PERFILES DE MAYOR RIESGO (INDIVIDUAL)
# ==============================================================================
# Se identifican los perfiles reales existentes en la muestra con la mayor
# probabilidad estimada, sin agrupar por rangos de edad.

# Seleccionamos las variables relevantes para caracterizar el perfil
# y ordenamos de mayor a menor probabilidad.
# Usamos distinct() para que si hay dos personas idénticas (mismo perfil), 
# no nos repita la fila, sino que muestre los 5 perfiles distintos más altos.

# 5. Creación de Perfiles de Vulnerabilidad
datos2024_final <- datos2024 %>%
  mutate(
    probabilidad_riesgo = prob_2024,
    perfil_vulnerable = ifelse(prob_2024 > umbral_youden, "SÍ", "NO")
  )

top10_riesgo <- datos2024_final %>%
  arrange(desc(probabilidad_riesgo)) %>% # Ordenar: el más alto primero
  select(
    ID_individuo,
    sexo, 
    P3,           # Edad exacta (sin agrupar)
    nivel_edu, 
    REGION, 
    jefe_hogar,
    estado_conyugal,
    identidad_afro, 
    identidad_indigena, 
    area,  
    asiste_escuela,
    probabilidad_riesgo
  ) %>%
  distinct(probabilidad_riesgo, .keep_all = TRUE) %>% 
  arrange(desc(probabilidad_riesgo)) %>% 
  head(20)     

print("--- TOP 10 PERFILES DE RIESGO (INDIVIDUOS REALES EN LA MUESTRA) ---")
print(top10_riesgo)

# ==============================================================================
# G. VISUALIZACIÓN CREATIVA (INSPIRADA EN REFERENCIA BIBLIOGRÁFICA)
# ==============================================================================
# 1. ANÁLISIS DESCRIPTIVO (LO OBSERVADO, NO LO PREDICHO)
# ------------------------------------------------------------------------------
# A. Tasa de Desempleo Real por Nivel Educativo
tasa_edu <- datos2024 %>%
  group_by(nivel_edu) %>%
  summarise(Tasa_Desempleo = mean(as.numeric(as.character(desempleado))) * 100)

g1 <- ggplot(tasa_edu, aes(x = reorder(nivel_edu, Tasa_Desempleo), y = Tasa_Desempleo)) +
  geom_col(fill = "#2c3e50", width = 0.7) +
  geom_text(aes(label = round(Tasa_Desempleo, 1)), vjust = -0.5, fontface = "bold") +
  labs(title = "Tasa de Desempleo Observada por Nivel Educativo (2024)",
       subtitle = "Datos directos de la muestra (Sin modelar)",
       x = "", y = "% Desempleo") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

print(g1)


# 2. CURVAS DE RESPUESTA (LO PREDICHO POR EL MODELO)
# ------------------------------------------------------------------------------
# Creamos datos para simular la vida de una persona de 15 a 80 años
curvas_data <- expand.grid(
  P3 = seq(15, 80, 1),
  nivel_edu = c("Educación Superior", "Educación Media", "Educación Básica General", "Ningún Grado"),
  # nivel_edu = c("Universitaria", "Vocacional", "Secundaria", "Primaria", "Ningun Grado"),
  sexo = "1",              # Fijamos en MUJER (grupo más vulnerable)
  REGION = "Metropolitana",# Fijamos región base
  jefe_hogar = "0",        # No jefe
  identidad_indigena = "0",
  identidad_afro = "0",
  estado_conyugal = "Soltero",
  asiste_escuela = "0",
  area = "0"
)

# IMPORTANTE: Calcular la edad centrada P3c para la predicción
curvas_data$P3c <- curvas_data$P3 - media_edad_2024

# Predecir
curvas_data$Probabilidad <- predict(modelo_2024, newdata = curvas_data, type = "response")
# Gráfico de Curvas Suaves (Estilo TFG pág. 34)
g2 <- ggplot(curvas_data, aes(x = P3, y = Probabilidad, color = nivel_edu, group = nivel_edu)) +
  geom_line(size = 1.5) +
  #scale_y_continuous(labels = scales::percent, limits = c(0, 0.55)) + # Ajusta límite si es necesario
  scale_x_continuous(limits = c(15, 80), breaks = seq(15, 80, by = 5)) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 0.5), breaks = seq(0, 0.5, by = 0.10)) +
  scale_color_manual(values = c(
    "Educación Superior" = "#51A3CC",
    "Educación Media"    = "#CC5151",
    "Educación Básica General"    = "#A3CC51",
    "Ningún Grado"      = "#6551CC"
  )) + 
  labs(title = "Curvas de Respuesta: Ciclo de Vida del Riesgo de Desempleo 2024",
    x = "Edad (Años)", y = "Probabilidad Estimada",
    color = "Nivel Educativo") +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "#AFE0D9", size = 0.5),
    panel.background = element_rect(fill = NA, color = NA),
    plot.background  = element_rect(fill = NA, color = NA),
        legend.position = c(0.85, 0.8),
        legend.background = element_rect(fill = "white"))
print(g2)

# 3. IMPACTO DEL GÉNERO (Curva Edad-Sexo)
# ------------------------------------------------------------------------------
sexo_data <- expand.grid(
  P3 = seq(15, 80, 1),
  sexo = c("0", "1"), # Hombre y Mujer
  nivel_edu = "Educación Superior", # Fijamos un nivel educativo común
  REGION = "Metropolitana",
  jefe_hogar = "0",
  identidad_indigena = "0",
  identidad_afro = "0",
  estado_conyugal = "Soltero",
  asiste_escuela = "0",
  area = "0"
)

sexo_data$P3c <- sexo_data$P3 - media_edad_2024
sexo_data$Probabilidad <- predict(modelo_2024, newdata = sexo_data, type = "response")
sexo_data$Sexo <- ifelse(sexo_data$sexo == "0", "Hombre", "Mujer")

g3 <- ggplot(sexo_data, aes(x = P3, y = Probabilidad, color = Sexo)) +
  geom_line(size = 1.2) +
  scale_color_manual(values = c("Mujer" = "#E74C3C", "Hombre" = "#3498DB")) +
  labs(
    title = "Brecha de Género Estructural (2024)",
       subtitle = "Diferencia de riesgo Masculino Vs. Femenino a lo largo de la vida",
    x = "Edad", y = "Probabilidad de Desempleo") +
  scale_x_continuous(limits = c(15, 80), breaks = seq(15, 80, by = 5)) +
  # ajuste de escala en eje Y
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0, 0.4), breaks = seq(0, 0.4 , by = 0.05)) +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(color = "#AFE0D9", size = 0.5),
    panel.background = element_rect(fill = NA, color = NA),
    plot.background  = element_rect(fill = NA, color = NA),
        legend.position = c(0.90, 0.9),
        legend.title = element_text(size = 12, face = "bold"),
        legend.text = element_text(size = 11, face = "bold"),
        legend.background = element_rect(fill = "white"))
print(g3)

################################################################################
# Curvas AUC
################################################################################
# 1. Calcular y redondear los valores de AUC
val_auc_2024 <- round(auc(roc_2024), 3)
val_auc_2021 <- round(auc(roc_2021), 3)
val_auc_2019 <- round(auc(roc_2019), 3)

# 2. Generar el gráfico base con el año 2024
plot(roc_2024, col = "#E74C3C", lwd = 3, 
     main = "Rendimiento del Modelo: Curvas ROC Comparativas",
     xlab = "Especificidad (1 - Tasa Falsos Positivos)",
     ylab = "Sensibilidad (Tasa Verdaderos Positivos)")

# 2. Generar el gráfico base con el año 2019
plot(roc_2019, col = "green4", lwd = 3, 
     main = "Rendimiento del Modelo: Curvas ROC Comparativas",
     xlab = "Especificidad (1 - Tasa Falsos Positivos)",
     ylab = "Sensibilidad (Tasa Verdaderos Positivos)")

# 3. Añadir las líneas de los otros años
lines(roc_2021, col = "#2b4c7e", lwd = 3)
lines(roc_2024, col = "#E74C3C", lwd = 3)

# 4. Leyenda profesional con los valores de AUC integrados
legend("bottomright", 
       inset = 0.05, 
       legend = c(paste("Año 2019 (AUC =", val_auc_2019, ")"),
                  paste("Año 2021 (AUC =", val_auc_2021, ")"), 
                  paste("Año 2024 (AUC =", val_auc_2024, ")")), 
       col = c("green4", "#2b4c7e", "#E74C3C"), 
       lwd = 4, 
       cex = 0.9,      # Tamaño de la letra
       bty = "n")      # Quita el borde de la caja para que se vea más moderno
################################################################################

    

```

---

Análisis desarrollado en RStudio Desktop

Versión estable utilizada: 2024.12.0+467 "Cranberry Hibiscus"
(Última versión disponible a la fecha de este estudio)

República de Panamá | 2026