1 Síntesis

Este estudio tiene como objetivo fortalecer el tratamiento de la no respuesta en encuestas económicas mediante la evaluación comparativa de diversos métodos de imputación. Se analizan métodos aplicados a dos variables clave: personal ocupado e ingresos.

Se comparan los métodos institucionales actuales —vecino más cercano ajustado por dominio–estrato— con alternativas multivariadas de mayor complejidad como MICE y k-NN. Asimismo, se optimizan hiperparámetros y se evalúa el desempeño de cada método mediante simulaciones controladas.

2 Introducción

El tratamiento adecuado de la no respuesta es fundamental para mantener la calidad de las estadísticas derivadas de encuestas económicas. En la EMEC, los métodos institucionales han sido diseñados para garantizar la coherencia histórica de las series; sin embargo, dichos métodos pueden presentar limitaciones ante patrones complejos de no respuesta.

Este documento analiza alternativas modernas de imputación multivariada y las contrasta con los enfoques tradicionales, utilizando para ello datos sintéticos con estructura estadística similar a la EMEC. El objetivo es identificar configuraciones que minimicen el error de imputación sin comprometer la coherencia temporal de las series económicas.

3 Paquetes y opciones

# --- cargamos paquetes necesarios ---
library(tidyverse)    # Manipulación de datos y gráficos (ggplot2)
library(openxlsx)     # Lectura de archivos Excel
library(naniar)       # Visualización de datos faltantes
library(finalfit)     # Prueba de Little (MCAR)
library(mice)         # Imputación Múltiple (MICE)
library(VIM)          # Imputación k-NN
library(missForest)   # Imputación Random Forest (Machine Learning)
library(Metrics)      # Cálculo de métricas de error (RMSE)
library(zoo)          # Funciones de series de tiempo (rollmean)
library(lubridate)    # Manejo de fechas

4 Carga y limpieza de datos

Dada la necesidad de evaluar métodos basados en series de tiempo —que requieren historial suficiente— se generó un conjunto sintético controlado que replica las características estadísticas de la EMEC. Se simulan tendencias, ruido y estructura estratificada tanto para ingresos como para personal ocupado.

Este enfoque permite evaluar el desempeño de los distintos métodos de imputación en un entorno donde se conoce la “verdad” de los datos, estableciendo así un marco justo para la comparación.

# --- generamos datos sintéticos ---
# creamos series de ingreso y personal para cada empresa con tendencias y ruido
set.seed(2025)
n_emp <- 50
# 36 meses para tener historia suficiente para variación anual
n_months <- 36 
start_period <- as.Date("2021-01-01")
periods <- seq.Date(from = start_period, by = "month", length.out = n_months)

rows <- list()

for(i in 1:n_emp){
  dominio <- sample(c(43,46), size = 1)
  estrato <- sample(c(1,2,3), size = 1)
  
  if(estrato == 1){ factor <- 1 } else { factor <- round(runif(1, 1.5, 25), 3) }

  # INGRESOS
  if(estrato == 1){
    base_start <- runif(1, 200000, 600000)
    trend <- rnorm(1, mean = 800, sd = 150)
    noise <- rnorm(n_months, mean = 0, sd = base_start * 0.10) 
  } else if(estrato == 2){
    base_start <- runif(1, 80000, 200000)
    trend <- rnorm(1, mean = 300, sd = 80)
    noise <- rnorm(n_months, mean = 0, sd = base_start * 0.06)
  } else {
    base_start <- runif(1, 30000, 100000)
    trend <- rnorm(1, mean = 150, sd = 60)
    noise <- rnorm(n_months, mean = 0, sd = base_start * 0.04)
  }
  
  ingresos_series <- base_start + cumsum(trend + noise)
  ingresos_series[ingresos_series < 1000] <- 1000
  
  # PERSONAL
  if(estrato == 1){
    emp_base <- rpois(1, lambda = 20) + 5
    emp_changes <- rnorm(n_months, mean = 0.5, sd = 1.5)
  } else if(estrato == 2){
    emp_base <- rpois(1, lambda = 8) + 2
    emp_changes <- rnorm(n_months, mean = 0.2, sd = 1.0)
  } else {
    emp_base <- rpois(1, lambda = 4) + 1
    emp_changes <- rnorm(n_months, mean = 0.1, sd = 0.8)
  }
  
  emp_series <- pmax(1, round(emp_base + cumsum(emp_changes)))

  for(j in seq_along(periods)){
    rows[[length(rows) + 1]] <- data.frame(
      clee = as.character(i), 
      fecha = periods[j],
      personal = as.integer(emp_series[j]),
      ingreso = ingresos_series[j],
      dominio = dominio,
      estrato = estrato,
      imputado = 0, 
      stringsAsFactors = FALSE
    )
  }
}

datos <- bind_rows(rows)
datos_reales <- datos 

cat("Datos sintéticos generados. Total:", nrow(datos), "\n")
## Datos sintéticos generados. Total: 1800
# --- función: promedio móvil ---
# calculamos la imputación usando rollapply con ventana k
imputar_promedio_movil <- function(data, variable, k) {
  data %>%
    group_by(clee) %>%
    arrange(fecha) %>%
    mutate(
      imputado = rollapply(get(variable), width = k, FUN = mean, 
                           align = "right", fill = NA, partial = FALSE, na.rm = TRUE)
    ) %>%
    ungroup() %>%
    pull(imputado)
}

# --- función: variación anual ---
# imputamos usando tasa anual por dominio-estrato
imputar_variacion_anual <- function(data, variable) {
  tasas <- data %>%
    group_by(dominio, estrato, fecha) %>%
    summarise(tasa = median(get(variable)/lag(get(variable), 12) - 1, na.rm=TRUE), .groups="drop")
  
  data %>%
    left_join(tasas, by=c("dominio","estrato","fecha")) %>%
    group_by(clee) %>%
    arrange(fecha) %>%
    mutate(
      imputado = lag(get(variable), 12) * (1 + replace_na(tasa, 0))
    ) %>%
    ungroup() %>%
    pull(imputado)
}

# --- función: knn estratificado ---
# imputamos usando vecino más cercano dentro de dominio y estrato
imputar_knn_estratificado <- function(data, variable) {
  res <- VIM::kNN(data, variable = variable, 
                  dist_var = c("dominio", "estrato"), 
                  k = 1, imp_var = FALSE)
  return(res[[variable]])
}

# --- función maestra: cascada institucional ---
# aplicamos variación anual, promedios móviles y al final knn si aún faltan valores
imputar_institucional_integrado <- function(data, variable) {
  
  tasas <- data %>%
    group_by(dominio, estrato, fecha) %>%
    summarise(g_sector = median(get(variable)/lag(get(variable), 12) - 1, na.rm=TRUE), .groups="drop")
  
  data_proc <- data %>%
    left_join(tasas, by=c("dominio","estrato","fecha")) %>%
    group_by(clee) %>%
    arrange(fecha) %>%
    mutate(
      c_var = lag(get(variable), 12) * (1 + replace_na(g_sector, 0)),
      c_ma6 = rollapply(get(variable), 6, mean, align="right", fill=NA, partial=FALSE, na.rm=TRUE),
      c_ma3 = rollapply(get(variable), 3, mean, align="right", fill=NA, partial=FALSE, na.rm=TRUE),
      c_ma2 = rollapply(get(variable), 2, mean, align="right", fill=NA, partial=FALSE, na.rm=TRUE)
    ) %>%
    ungroup()
  
  data_proc$final <- coalesce(data_proc[[variable]], data_proc$c_var, data_proc$c_ma6, data_proc$c_ma3, data_proc$c_ma2)
  
  faltan <- which(is.na(data_proc$final))
  if(length(faltan) > 0) {
    knn_imp <- VIM::kNN(data_proc, variable=variable, dist_var=c("dominio","estrato"), k=1, imp_var=FALSE)
    data_proc$final[faltan] <- knn_imp[[variable]][faltan]
  }
  
  return(data_proc$final)
}

5 Exploración inicial

El resumen descriptivo muestra valores coherentes con la estructura simulada y confirma correlación positiva entre personal ocupado e ingresos.

En la gráfica de dispersión se observa una relación creciente entre ambas variables, lo cual es consistente con la dinámica económica simulada.

A continuación se presentan hallazgos derivados de la optimización del algoritmo MICE.

# --- exploración inicial de datos ---
# resumen y correlaciones
datos %>% select(personal, ingreso) %>% summary()
##     personal        ingreso      
##  Min.   : 1.00   Min.   : 30665  
##  1st Qu.: 9.00   1st Qu.: 77803  
##  Median :17.00   Median :162205  
##  Mean   :20.43   Mean   :216795  
##  3rd Qu.:31.00   3rd Qu.:325511  
##  Max.   :64.00   Max.   :867190
# correlación entre personal y ingreso
correlacion <- cor(datos$personal, datos$ingreso, use = "complete.obs")
correlacion
## [1] 0.6684274
# scatter
ggplot(datos, aes(x = personal, y = ingreso)) +
  geom_point(alpha = 0.4, color = "darkblue") +
  labs(title = paste0("Correlación: PO vs ING (r = ", round(correlacion, 2), ")"),
       subtitle = "Relación entre Personal Ocupado e Ingresos",
       x = "PO (Personal Ocupado)",
       y = "ING (Ingresos)") +
  theme_minimal()

# --- distribución por estrato ---
# vemos la heterogeneidad entre grupos para justificar métodos estratificados

ggplot(datos, aes(x = factor(estrato), y = ingreso, fill = factor(estrato))) +
  geom_boxplot(alpha = 0.6) +
  
  # --- escala log ---
  # colocamos escala log para visualizar mejor diferencias grandes
  scale_y_log10(labels = scales::comma) + 
  
  labs(title = "Distribución de Ingresos por Estrato",
       subtitle = "Se observa una clara diferenciación: Estrato 1 concentra mayores ingresos",
       x = "Estrato (1=Grande, 3=Pequeño)", 
       y = "Ingresos (Escala Logarítmica)",
       fill = "Estrato") +
  theme_minimal()

El gráfico de cajas y bigotes compara la variable Ingresos (eje Y, en escala logarítmica) desagregada por Estrato (eje X).

6 Diseño del experimento (simulaciones)

Se define un escenario MCAR (Missing Completely at Random) borrando el 20% de los datos de forma aleatoria para evaluar la capacidad de recuperación de cada método.

# --- preparamos ground truth ---
# tomamos registros completos para evaluar los métodos
datos_reales <- datos %>% 
  select(clee, fecha, personal, ingreso, dominio, estrato) %>% 
  drop_na(personal, ingreso)

n <- nrow(datos_reales)
cat("n registros reales (completos):", n, "\n")
## n registros reales (completos): 1800

6.1 Simulación MCAR

# --- simulamos mecanismo mar ---
# la probabilidad depende de estrato, tamaño y dominio
set.seed(2025)
prop_faltantes <- 0.20
idx_personal_mcar <- sample(1:n, size = floor(prop_faltantes * n))
idx_ingreso_mcar  <- sample(1:n, size = floor(prop_faltantes * n)) # Cambio de nombre

datos_mcar <- datos_reales
datos_mcar$personal[idx_personal_mcar] <- NA
datos_mcar$ingreso[idx_ingreso_mcar] <- NA

6.2 Simulación MAR

# --- simulación mar ---
# la probabilidad de falta depende de estrato, tamaño y dominio
set.seed(2025)

datos_mar <- datos_reales

n <- nrow(datos_mar)


# --- Probabilidad según estrato ---
# Estrato 1 suele ser el que más no responde
p_estrato <- case_when(
  datos_mar$estrato == 1 ~ 0.30,   # pequeñas
  datos_mar$estrato == 2 ~ 0.20,
  datos_mar$estrato == 3 ~ 0.10,
  TRUE ~ 0.05
)


# --- Probabilidad según tamaño ---
cuartil_personal <- ntile(datos_mar$personal, 4)

p_personal <- case_when(
  cuartil_personal == 1 ~ 0.30,   # empresas chicas responden menos
  cuartil_personal == 2 ~ 0.20,
  cuartil_personal == 3 ~ 0.10,
  cuartil_personal == 4 ~ 0.05,   # grandes responden más
  TRUE ~ 0.15
)

# --- Probabilidad según dominio ---
p_dominio <- ifelse(datos_mar$dominio == 46, 0.25, 0.10)

# --- Combinado de las probabilidades ---
# Tomamos un promedio ponderado (puedes ajustar los pesos)
prob_falta_total <- 0.5 * p_estrato + 
                    0.3 * p_personal + 
                    0.2 * p_dominio

# Para asegurar que max no exceda 1
prob_falta_total <- pmin(prob_falta_total, 0.90)

# --- Aplicar NA en ingreso y personal ---
random_uniform <- runif(n)

datos_mar$ingreso[ random_uniform < prob_falta_total ] <- NA
datos_mar$personal[ runif(n) < prob_falta_total ] <- NA
# 1. Unificar los datasets simulados para compararlos
df_comparativo <- bind_rows(
  datos_mcar %>% mutate(Mecanismo = "MCAR"),
  datos_mar  %>% mutate(Mecanismo = "MAR")
) %>%
  mutate(
    # Recalculamos cuartiles sobre el dato de personal
    cuartil_personal = ntile(personal, 4)
  )

# 2. Calcular la proporción de NAs en 'ingreso' por grupo
heatmap_data <- df_comparativo %>%
  group_by(Mecanismo, estrato, cuartil_personal) %>%
  summarise(
    # AQUÍ ESTÁ LA CLAVE: Usamos 'ingreso' que es como se llama en tus datos sintéticos
    prop_na = mean(is.na(ingreso)), 
    .groups = "drop"
  )

# 3. Generar el Mapa de Calor (Heatmap)
ggplot(heatmap_data, aes(x = factor(cuartil_personal), y = factor(estrato), fill = prop_na)) +
  geom_tile(color = "white") + # Bordes blancos
  facet_wrap(~Mecanismo) +     # Dividir paneles
  
  # Escala de colores (Rojo para NA)
  scale_fill_gradient(low = "#fee5d9", high = "#de2d26", 
                      limits = c(0, 1), name = "Prop. NA") +
  
  # Etiquetas
  labs(title = "Comparación MCAR vs MAR",
       subtitle = "Proporción de NA en 'ingresos'",
       x = "Cuartil de Personal",
       y = "Estrato") +
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    axis.text = element_text(color = "gray30")
  )

El mapa de calor compara la estructura de los datos faltantes en la variable Ingresos bajo dos mecanismos simulados: MCAR y MAR.

  • Ejes del gráfico:

    • Eje X: Tamaño de la empresa medido mediante Cuartiles de Personal.
    • Eje Y: Estrato económico de la empresa.

  • Panel MCAR (derecha): Los colores aparecen dispersos sin un patrón definido. Las intensidades son suaves e irregulares, lo que es consistente con un mecanismo completamente aleatorio donde los datos se pierden sin relación con las características de la empresa.

  • Panel MAR (izquierda): Se aprecia una estructura claramente definida. Destaca un bloque de color más intenso en el Estrato 2, Cuartil 4, indicando una concentración sistemática de no respuesta en un grupo específico.

  • Hallazgo: El gráfico valida la correcta implementación del diseño experimental.

    • En MCAR, los faltantes provienen de variabilidad aleatoria.
    • En MAR, la falta de respuesta depende de características observables (tamaño y estrato), reproduciendo condiciones realistas en encuestas económicas.

  • Implicaciones para la imputación: El patrón estructurado del escenario MAR constituye el principal desafío metodológico. Cualquier método que no sea capaz de manejar este sesgo (ej. promedios globales) producirá imputaciones distorsionadas.

  • Justificación metodológica: La presencia de un bloque sistemático de no respuesta en MAR confirma la necesidad de utilizar métodos avanzados, como Random Forest, que pueden identificar patrones complejos y generar imputaciones coherentes según el perfil de cada empresa. Esto es especialmente importante para instituciones como el INEGI, donde ignorar dicha estructura generaría estimaciones sesgadas.

7 Métodos de imputación (implementación)

Se implementan los métodos: - Institucionales: Promedios Móviles, Variación Anual y Estrategia en Cascada. - Alternativos: Media, MICE, k-NN y Random Forest.

# --- función: rmse ---
# calculamos rmse solo sobre posiciones donde hubo na
calcular_rmse <- function(real, imputado, mascara) {
  # rmse sobre las posiciones donde mascara==TRUE
  if(sum(mascara) == 0) return(NA)
  return(rmse(real[mascara], imputado[mascara]))
}

7.1 Optimización de hiperparámetros

7.1.1 MICE: búsqueda de maxit y m (iteraciones e imputaciones)

# --- optimizamos parámetros de mice ---
# probamos combinaciones de maxit y m sobre datos mcar
df_opt <- datos_mcar
mascara_na <- as.data.frame(is.na(df_opt))

iteraciones <- c(5, 10, 20)
m_values <- c(5, 10)
res_mice <- tibble()

for(it in iteraciones) {
  for(mv in m_values) {
    imp <- try(mice(df_opt, method = 'pmm', m = mv, maxit = it, printFlag = FALSE, seed = 123), silent = TRUE)
    if(inherits(imp, "try-error")) next
    comp <- mice::complete(imp)
    rmse_p <- calcular_rmse(datos_reales$personal, comp$personal, mascara_na$personal)
    rmse_m <- calcular_rmse(datos_reales$ingreso, comp$ingreso, mascara_na$ingreso)
    res_mice <- bind_rows(res_mice, tibble(maxit = it, m = mv, rmse_personal = rmse_p, rmse_ingreso = rmse_m, rmse_promedio = mean(c(rmse_p, rmse_m), na.rm = TRUE)))
  }
}

# --- GRÁFICA CON DESTAQUE (CÍRCULO ROJO) ---
mejor_config <- res_mice %>% filter(rmse_promedio == min(rmse_promedio))

ggplot(res_mice, aes(x = factor(maxit), y = rmse_promedio, group = factor(m), color = factor(m))) +
  geom_line(size = 0.8) + 
  geom_point(size = 2) +
  # Círculo rojo en el mínimo
  geom_point(data = mejor_config, aes(x = factor(maxit), y = rmse_promedio), 
             color = "red", size = 5, shape = 1, stroke = 1.5) +
  # Etiqueta
  geom_label(data = mejor_config, 
             aes(label = paste0("Min RMSE:\n", format(round(rmse_promedio, 0), big.mark=","))),
             vjust = -0.5, color = "black", fontface = "bold", fill="white", alpha=0.9) +
  labs(title = "Optimización MICE: Selección de Hiperparámetros", 
       subtitle = paste0("Configuración óptima detectada: m=", mejor_config$m, " | maxit=", mejor_config$maxit),
       x = "Iteraciones (maxit)", y = "RMSE Promedio", color = "m") +
  theme_minimal() + theme(legend.position = "bottom")

Los valores m = 5 y m = 10 se seleccionaron siguiendo las recomendaciones clásicas de Rubin (1987) y van Buuren (2018), quienes sugieren trabajar con entre 5 y 10 imputaciones para obtener estimadores estables sin incrementar innecesariamente el costo computacional.

Se evaluaron maxit = 5, 10 y 20 porque MICE usualmente converge entre 5 y 10 iteraciones (van Buuren, 2018). Incluir 20 permite examinar la estabilidad de la cadena y verificar si un número mayor de iteraciones reduce adicionalmente el RMSE.

El gráfico de líneas evalúa el desempeño del algoritmo MICE midiendo el RMSE promedio (eje Y) bajo distintas combinaciones de dos parámetros:
- maxit (eje X): número de iteraciones del algoritmo.
- m (color de las líneas): cantidad de imputaciones múltiples.

  • Línea azul (\(m = 10\)):
    • Presenta el mejor desempeño global con el error más bajo cuando maxit = 5.
    • Sin embargo, muestra inestabilidad: el RMSE aumenta abruptamente al pasar a 10 iteraciones y luego disminuye moderadamente al llegar a 20.
  • Línea roja (\(m = 5\)):
    • Exhibe el patrón contrario: inicia con un RMSE alto en 5 iteraciones.
    • Mejora de manera sustancial en maxit = 10, donde alcanza su punto óptimo y se estabiliza.
  • Hallazgos:
    • Punto óptimo: Se logró reducir el error al mínimo de todo el experimento, alcanzando un RMSE promedio de 69,283. Este valor se obtuvo con la configuración \(m = 10\) y \(maxit = 5\), lo que indica que el algoritmo converge rápidamente sin requerir una carga computacional excesiva.
    • Rendimiento decreciente: Incrementar el número de iteraciones no garantiza mejoras. En el caso de \(m = 10\), aumentar el procesamiento (pasar de 5 a 10 iteraciones) introdujo ruido que elevó el error, alejándose del óptimo encontrado.

7.1.2 k-NN: búsqueda de k

# --- optimizamos knn ---
# probamos valores de k y evaluamos rmse
k_values <- c(3, 5, 10, 20)
res_knn <- tibble()

for(k in k_values) {
  imp_knn <- kNN(df_opt, variable = c("personal", "ingreso"), k = k, imp_var = FALSE)
  rmse_p <- calcular_rmse(datos_reales$personal, imp_knn$personal, mascara_na$personal)
  rmse_m <- calcular_rmse(datos_reales$ingreso, imp_knn$ingreso, mascara_na$ingreso)
  res_knn <- bind_rows(res_knn, tibble(k = k, rmse_personal = rmse_p, rmse_ingreso = rmse_m, rmse_promedio = mean(c(rmse_p, rmse_m), na.rm = TRUE)))
}
##   ingreso   dominio   estrato   ingreso   dominio   estrato 
##  30664.53     43.00      1.00 867190.28     46.00      3.00 
## personal  dominio  estrato personal  dominio  estrato 
##        1       43        1       64       46        3 
##   ingreso   dominio   estrato   ingreso   dominio   estrato 
##  30664.53     43.00      1.00 867190.28     46.00      3.00 
## personal  dominio  estrato personal  dominio  estrato 
##        1       43        1       64       46        3 
##   ingreso   dominio   estrato   ingreso   dominio   estrato 
##  30664.53     43.00      1.00 867190.28     46.00      3.00
## personal  dominio  estrato personal  dominio  estrato 
##        1       43        1       64       46        3 
##   ingreso   dominio   estrato   ingreso   dominio   estrato 
##  30664.53     43.00      1.00 867190.28     46.00      3.00
## personal  dominio  estrato personal  dominio  estrato 
##        1       43        1       64       46        3
# --- GRÁFICA CON DESTAQUE ---
mejor_knn <- res_knn %>% filter(rmse_promedio == min(rmse_promedio))

ggplot(res_knn, aes(x = k, y = rmse_promedio)) + 
  geom_line() + geom_point() +
  # Círculo rojo
  geom_point(data = mejor_knn, aes(x = k, y = rmse_promedio), 
             color = "red", size = 5, shape = 1, stroke = 1.5) +
  geom_label(data = mejor_knn, 
             aes(label = paste0("Min RMSE:\n", format(round(rmse_promedio, 0), big.mark=","))), 
             vjust = -0.5, fontface = "bold") +
  labs(title = "Optimización k-NN: selección de k", 
       subtitle = paste0("El error se minimiza en k=", mejor_knn$k),
       x = "Vecinos (k)", y = "RMSE Promedio") + 
  theme_minimal()

El gráfico de línea muestra la evolución del RMSE promedio (eje Y) conforme aumenta el número de vecinos (\(k\), eje X) utilizados por el algoritmo k-NN.

  • Descenso inicial:
    • El error disminuye abruptamente al pasar de \(k = 3\) (RMSE > 27,000) a \(k = 5\), evidenciando una rápida mejora al permitir que el modelo considere más información vecinal.
  • Punto mínimo (“codo”):
    • El RMSE alcanza su valor más bajo en \(k = 10\), con un error aproximado de 24,200.
    • Este punto representa el equilibrio óptimo entre captar estructura y evitar ruido.
  • Ascenso posterior:
    • Al incrementar hasta \(k = 20\), el error vuelve a aumentar, lo que indica un deterioro del rendimiento del modelo.
  • Hallazgo analítico:
    • El comportamiento del gráfico refleja de manera clara el compromiso entre sesgo y varianza (bias–variance trade-off):
      • \(k\) bajos (3–5): El modelo sufre sobreajuste (overfitting). Es demasiado sensible al ruido local y a valores atípicos.
      • \(k\) altos (20): Aparece subajuste (underfitting). La predicción se vuelve excesivamente suave al promediar demasiados vecinos.
  • Conclusión operativa:
    • El valor óptimo para la imputación es k = 10, logrando minimizar el error a un RMSE de 24,172. Esta configuración captura eficientemente la estructura local de los datos económicos sin verse afectada por la variabilidad extrema de las empresas individuales.

8 Ajuste final de modelos (usar mejores hiperparámetros)

# --- preparamos dataset base ---
# añadimos variables auxiliares para los modelos

df_base <- datos_mcar 
df_base$clee <- datos_reales$clee
df_base$fecha <- datos_reales$fecha
df_base$dominio <- datos_reales$dominio
df_base$estrato <- datos_reales$estrato

mascara_na <- as.data.frame(is.na(df_base))

# --- imputación media ---
# reemplazamos por medias por dominio-estrato
imp_media <- df_base %>%
  group_by(dominio, estrato) %>%
  mutate(
    personal = ifelse(is.na(personal), mean(personal, na.rm=TRUE), personal),
    ingreso  = ifelse(is.na(ingreso),  mean(ingreso, na.rm=TRUE),  ingreso)
  ) %>% ungroup()
# Relleno de seguridad por si algún estrato es todo NA
imp_media$personal[is.na(imp_media$personal)] <- mean(df_base$personal, na.rm=TRUE)
imp_media$ingreso[is.na(imp_media$ingreso)] <- mean(df_base$ingreso, na.rm=TRUE)

# --- imputación mice ---
# usamos parámetros optimizados: m = 1, maxit = 1
cols_mice <- c("personal", "ingreso", "dominio", "estrato")
imp_mice_mod <- mice(df_base[, cols_mice], method='pmm', m=10, maxit=5, printFlag=FALSE, seed=123)
imp_mice_data <- mice::complete(imp_mice_mod)
imp_mice <- df_base
imp_mice$personal <- imp_mice_data$personal
imp_mice$ingreso <- imp_mice_data$ingreso

# --- imputación knn ---
# usamos knn con k=5 y distancia dominio-estrato
imp_knn <- kNN(df_base, variable=c("personal","ingreso"), dist_var=c("dominio","estrato"), k=10, imp_var=FALSE)
## dominio estrato dominio estrato 
##      43       1      46       3
## dominio estrato dominio estrato 
##      43       1      46       3
# --- Random Forest ---
# Aseguramos factores para que RF entienda grupos
df_rf_in <- df_base[, c("personal","ingreso","dominio","estrato")]
df_rf_in$dominio <- as.factor(df_rf_in$dominio)
df_rf_in$estrato <- as.factor(df_rf_in$estrato)
imp_rf_obj <- missForest(df_rf_in, verbose=FALSE)
imp_rf <- df_base
imp_rf$personal <- imp_rf_obj$ximp$personal
imp_rf$ingreso <- imp_rf_obj$ximp$ingreso

# --- MÉTODOS INSTITUCIONALES (INDIVIDUALES) ---
# Generamos los vectores necesarios para la tabla detallada
df_inst <- df_base

# Personal
per_pm2 <- imputar_promedio_movil(df_inst, "personal", 2)
per_pm3 <- imputar_promedio_movil(df_inst, "personal", 3)
per_pm6 <- imputar_promedio_movil(df_inst, "personal", 6)
per_anual <- imputar_variacion_anual(df_inst, "personal")
per_knn_est <- imputar_knn_estratificado(df_inst, "personal")
## dominio estrato dominio estrato 
##      43       1      46       3
# Ingreso
ing_pm2 <- imputar_promedio_movil(df_inst, "ingreso", 2)
ing_pm3 <- imputar_promedio_movil(df_inst, "ingreso", 3)
ing_pm6 <- imputar_promedio_movil(df_inst, "ingreso", 6)
ing_anual <- imputar_variacion_anual(df_inst, "ingreso")
ing_knn_est <- imputar_knn_estratificado(df_inst, "ingreso")
## dominio estrato dominio estrato 
##      43       1      46       3
# --- ESTRATEGIA INSTITUCIONAL INTEGRADA (CASCADA) ---
# Solución completa
ingreso_cascada <- imputar_institucional_integrado(df_inst, "ingreso")
## dominio estrato dominio estrato 
##      43       1      46       3
# --- CONSOLIDACIÓN ---
calc_rmse_safe <- function(real, estimado, mask) {
  indices <- which(mask & !is.na(estimado))
  if(length(indices) == 0) return(NA)
  return(rmse(real[indices], estimado[indices]))
}

res_final <- tibble(
  metodo = c("Media (Estrato)", "MICE", "k-NN", "Random Forest", 
             "Prom. Movil 2", "Prom. Movil 3", "Prom. Movil 6", 
             "Var. Anual", "Vecino (Estrato)", "Estrategia Cascada (INEGI)"),
  
  rmse_personal = c(
    calc_rmse_safe(datos_reales$personal, imp_media$personal, mascara_na$personal),
    calc_rmse_safe(datos_reales$personal, imp_mice$personal, mascara_na$personal),
    calc_rmse_safe(datos_reales$personal, imp_knn$personal, mascara_na$personal),
    calc_rmse_safe(datos_reales$personal, imp_rf$personal, mascara_na$personal),
    calc_rmse_safe(datos_reales$personal, per_pm2, mascara_na$personal),
    calc_rmse_safe(datos_reales$personal, per_pm3, mascara_na$personal),
    calc_rmse_safe(datos_reales$personal, per_pm6, mascara_na$personal),
    calc_rmse_safe(datos_reales$personal, per_anual, mascara_na$personal),
    calc_rmse_safe(datos_reales$personal, per_knn_est, mascara_na$personal),
    NA # No calculamos cascada para personal en este ejemplo
  ),
  
  rmse_ingreso = c(
    calc_rmse_safe(datos_reales$ingreso, imp_media$ingreso, mascara_na$ingreso),
    calc_rmse_safe(datos_reales$ingreso, imp_mice$ingreso, mascara_na$ingreso),
    calc_rmse_safe(datos_reales$ingreso, imp_knn$ingreso, mascara_na$ingreso),
    calc_rmse_safe(datos_reales$ingreso, imp_rf$ingreso, mascara_na$ingreso),
    calc_rmse_safe(datos_reales$ingreso, ing_pm2, mascara_na$ingreso),
    calc_rmse_safe(datos_reales$ingreso, ing_pm3, mascara_na$ingreso),
    calc_rmse_safe(datos_reales$ingreso, ing_pm6, mascara_na$ingreso),
    calc_rmse_safe(datos_reales$ingreso, ing_anual, mascara_na$ingreso),
    calc_rmse_safe(datos_reales$ingreso, ing_knn_est, mascara_na$ingreso),
    calc_rmse_safe(datos_reales$ingreso, ingreso_cascada, mascara_na$ingreso)
  )
)
# --- análisis de cobertura ---
# ver qué porcentaje de los huecos pudo llenar cada método

calcular_cobertura <- function(imp_vec, mask_vec) {
  total_huecos <- sum(mask_vec)
  llenos <- sum(mask_vec & !is.na(imp_vec))
  return((llenos / total_huecos) * 100)
}

res_final <- res_final %>%
  mutate(
    cob_personal = c(
      calcular_cobertura(imp_media$personal, mascara_na$personal),
      calcular_cobertura(imp_mice$personal, mascara_na$personal),
      calcular_cobertura(imp_knn$personal, mascara_na$personal),
      calcular_cobertura(imp_rf$personal, mascara_na$personal),
      calcular_cobertura(per_pm2, mascara_na$personal), 
      calcular_cobertura(per_pm3, mascara_na$personal), 
      calcular_cobertura(per_pm6, mascara_na$personal), 
      calcular_cobertura(per_anual, mascara_na$personal), 
      calcular_cobertura(per_knn_est, mascara_na$personal),
      NA # <--- 10. ESTRATEGIA CASCADA (No la aplicamos a Personal en este ejemplo)
    ),
    cob_ingreso = c(
      calcular_cobertura(imp_media$ingreso, mascara_na$ingreso),
      calcular_cobertura(imp_mice$ingreso, mascara_na$ingreso),
      calcular_cobertura(imp_knn$ingreso, mascara_na$ingreso),
      calcular_cobertura(imp_rf$ingreso, mascara_na$ingreso),
      calcular_cobertura(ing_pm2, mascara_na$ingreso),
      calcular_cobertura(ing_pm3, mascara_na$ingreso),
      calcular_cobertura(ing_pm6, mascara_na$ingreso),
      calcular_cobertura(ing_anual, mascara_na$ingreso),
      calcular_cobertura(ing_knn_est, mascara_na$ingreso),
      calcular_cobertura(ingreso_cascada, mascara_na$ingreso) # <--- 10. ESTRATEGIA CASCADA (Aquí sí la calculamos)
    )
  )

knitr::kable(res_final, digits = 2, caption = "Tabla 1: Resumen de Desempeño (RMSE y Cobertura)")
Tabla 1: Resumen de Desempeño (RMSE y Cobertura)
metodo rmse_personal rmse_ingreso cob_personal cob_ingreso
Media (Estrato) 6.80 113006.6 100.00 100.00
MICE 9.71 158663.5 100.00 100.00
k-NN 7.58 140013.3 100.00 100.00
Random Forest 6.48 112142.1 100.00 100.00
Prom. Movil 2 18.79 239838.9 93.06 92.22
Prom. Movil 3 18.65 241396.6 93.89 95.28
Prom. Movil 6 18.21 237068.6 85.56 86.94
Var. Anual 17.25 202261.4 53.89 50.56
Vecino (Estrato) 8.17 118473.2 100.00 100.00
Estrategia Cascada (INEGI) NA 237806.1 NA 100.00

Análisis de desempeño (RMSE) y viabilidad operativa (Cobertura)

# --- gráfica comparativa de rmse ---
# visualización de barras para comparar desempeño entre métodos

# Convertimos los datos a formato largo (ya lo tenías bien)
res_final_long <- res_final %>% 
  pivot_longer(cols = starts_with("rmse_"), 
               names_to = "variable", 
               values_to = "rmse") %>%
  # Opcional: Quitamos el prefijo "rmse_" para que los títulos se vean más limpios
  mutate(variable = str_to_title(str_remove(variable, "rmse_")))

# Creamos el gráfico con facetas
ggplot(res_final_long, aes(x = metodo, y = rmse, fill = variable)) +
  # Usamos geom_col simple, ya no necesitamos position="dodge"
  geom_col() +
  facet_wrap(~variable, scales = "free_y", ncol = 1) +
  
  # Estética y etiquetas ---
  scale_fill_manual(values = c("Ingreso" = "#F8766D", "Personal" = "#00BFC4")) +
  
  labs(title = "Comparativa Final de RMSE por Método",
       subtitle = "Nota: Las escalas del eje Y son diferentes para cada variable.",
       y = "RMSE (Valor del error)",
       x = NULL) + # Quitamos etiqueta X para no saturar
  
  theme_minimal() +
  theme(
    # Rotamos las etiquetas del eje X para que se lean bien los nombres largos
    axis.text.x = element_text(angle = 45, hjust = 1, size = 9),
    # Ocultamos la leyenda porque el título de cada panel ya dice qué es
    legend.position = "none",
    # Negritas para los títulos de los paneles
    strip.text = element_text(face = "bold", size = 12)
  )

Comparativa Final de RMSE por Método (Desagregado)

Lo que significa (Hallazgo analítico): Este gráfico visualiza la superioridad de los métodos transversales (ML) frente a los longitudinales en escenarios de información limitada:

8.1 Comparacion MAR MCAR

# --- PROCESAMIENTO RÁPIDO DE MAR ---
# Preparamos datos MAR
df_mar_proc <- datos_mar
df_mar_proc$clee    <- datos_reales$clee
df_mar_proc$fecha   <- datos_reales$fecha
df_mar_proc$dominio <- datos_reales$dominio
df_mar_proc$estrato <- datos_reales$estrato
mask_mar <- as.data.frame(is.na(df_mar_proc %>% select(personal, ingreso)))

# Ejecución de modelos en MAR
# A. Media
imp_media_mar <- df_mar_proc %>% group_by(dominio, estrato) %>% 
  mutate(ingreso = ifelse(is.na(ingreso), mean(ingreso, na.rm=TRUE), ingreso),
         personal = ifelse(is.na(personal), mean(personal, na.rm=TRUE), personal)) %>% ungroup()
imp_media_mar$ingreso[is.na(imp_media_mar$ingreso)] <- mean(datos_reales$ingreso, na.rm=TRUE)

# B. Random Forest
vars_rf_mar <- df_mar_proc %>% select(personal, ingreso, dominio, estrato) %>%
  mutate(dominio = as.factor(dominio), estrato = as.factor(estrato))
set.seed(2025)
imp_rf_mar <- missForest(vars_rf_mar, maxiter=5, ntree=50, verbose=FALSE)$ximp

# C. Cascada (Solo Ingreso)
cascada_mar <- imputar_institucional_integrado(df_mar_proc, "ingreso")
## dominio estrato dominio estrato 
##      43       1      46       3
# D. k-NN
imp_knn_mar <- VIM::kNN(df_mar_proc, variable=c("personal","ingreso"), dist_var=c("dominio","estrato"), k=10, imp_var=FALSE)
## dominio estrato dominio estrato 
##      43       1      46       3
## dominio estrato dominio estrato 
##      43       1      46       3
# --- CONSTRUCCIÓN DE LA TABLA MAR ---
rmse_mar_calc <- function(real, est, mask) {
  idx <- which(mask & !is.na(est))
  if(length(idx)==0) return(NA)
  return(rmse(real[idx], est[idx]))
}

res_mar_table <- tibble(
  metodo = c("Media (Estrato)", "Random Forest", "k-NN", "Estrategia Cascada (INEGI)"),
  rmse_personal = c(
    rmse_mar_calc(datos_reales$personal, imp_media_mar$personal, mask_mar$personal),
    rmse_mar_calc(datos_reales$personal, imp_rf_mar$personal, mask_mar$personal),
    rmse_mar_calc(datos_reales$personal, imp_knn_mar$personal, mask_mar$personal),
    NA
  ),
  rmse_ingreso = c(
    rmse_mar_calc(datos_reales$ingreso, imp_media_mar$ingreso, mask_mar$ingreso),
    rmse_mar_calc(datos_reales$ingreso, imp_rf_mar$ingreso, mask_mar$ingreso),
    rmse_mar_calc(datos_reales$ingreso, imp_knn_mar$ingreso, mask_mar$ingreso),
    rmse_mar_calc(datos_reales$ingreso, cascada_mar, mask_mar$ingreso)
  )
) %>%
  mutate(
    rmse_promedio = rowMeans(select(., rmse_personal, rmse_ingreso), na.rm = TRUE),
    Mecanismo = "MAR"
  )

# --- UNIÓN CON TU TABLA 'res_final' (MCAR) ---
res_mcar_table <- res_final %>%
  filter(metodo %in% c("Media (Estrato)", "Random Forest", "k-NN", "Estrategia Cascada (INEGI)")) %>%
  mutate(
    rmse_promedio = rowMeans(select(., rmse_personal, rmse_ingreso), na.rm = TRUE),
    Mecanismo = "MCAR"
  ) %>%
  select(Mecanismo, metodo, rmse_promedio, rmse_personal, rmse_ingreso)

tabla_consolidada <- bind_rows(res_mcar_table, res_mar_table) %>%
  arrange(Mecanismo, rmse_promedio)

# Imprimimos tabla
knitr::kable(tabla_consolidada %>% 
               mutate(across(where(is.numeric), ~ format(round(., 2), big.mark=",", scientific=FALSE))), 
             caption = "Comparación Estructural: MCAR vs MAR")
Comparación Estructural: MCAR vs MAR
Mecanismo metodo rmse_promedio rmse_personal rmse_ingreso
MAR Random Forest 57,617.47 6.42 115,228.5
MAR Media (Estrato) 59,152.50 6.60 118,298.4
MAR k-NN 65,117.81 7.23 130,228.4
MAR Estrategia Cascada (INEGI) 238,408.37 NA 238,408.4
MCAR Random Forest 56,074.30 6.48 112,142.1
MCAR Media (Estrato) 56,506.72 6.80 113,006.6
MCAR k-NN 70,010.42 7.58 140,013.3
MCAR Estrategia Cascada (INEGI) 237,806.07 NA 237,806.1 
# --- GRÁFICA ---
ggplot(tabla_consolidada, aes(x = reorder(metodo, rmse_ingreso), y = rmse_ingreso, fill = Mecanismo)) +
  geom_col(position = "dodge") +
  coord_flip() +
  facet_wrap(~Mecanismo) +
  scale_fill_manual(values = c("MCAR" = "#F8766D", "MAR" = "#CD5C5C")) + 
  labs(title = "Impacto del Mecanismo de No Respuesta (Solo Ingresos)",
       subtitle = "Note el incremento del error en el escenario MAR (derecha)",
       y = "RMSE (Ingresos)", x = NULL) +
  theme_minimal() +
  theme(legend.position = "none", strip.text = element_text(face="bold", size=12))

Se observa una vulnerabilidad Crítica de la Estrategia en Cascada: La Estrategia en Cascada (barra superior) presenta el error más alto en ambos escenarios, con un RMSE que supera los 230,000 . El hecho de que su error sea masivo y casi idéntico en ambos paneles revela que su debilidad principal no es el tipo de no respuesta, sino la falta de historial (“arranque en frío”). Al no contar con los 12 meses previos de datos, el método falla estructuralmente, sin importar si el dato falta por azar (MCAR) o por características de la empresa (MAR).

Incremento del Error en Métodos Simples: Aunque sutil visualmente, los métodos como k-NN y la Media por Estrato muestran variaciones en su desempeño al pasar de un escenario a otro. Sin embargo, no logran igualar la minimización del error que ofrece Random Forest, el cual se confirma como la opción más segura ante la incertidumbre del mecanismo de pérdida de datos.

# --- GRÁFICA DE DENSIDAD AVANZADA (MCAR vs MAR) ---

# 1. Preparamos datos de Random Forest (el mejor modelo) para ambos escenarios
df_densidad_final <- bind_rows(
  # Datos MCAR (que ya tenías en 'imp_rf')
  tibble(Valor = datos_reales$ingreso, Tipo = "Real", Variable = "Ingresos", Mecanismo = "MCAR"),
  tibble(Valor = imp_rf$ingreso, Tipo = "Imputado", Variable = "Ingresos", Mecanismo = "MCAR"),
  tibble(Valor = datos_reales$personal, Tipo = "Real", Variable = "Personal Ocupado", Mecanismo = "MCAR"),
  tibble(Valor = imp_rf$personal, Tipo = "Imputado", Variable = "Personal Ocupado", Mecanismo = "MCAR"),
  
  # Datos MAR (Calculados en el paso anterior 'imp_rf_mar')
  tibble(Valor = datos_reales$ingreso, Tipo = "Real", Variable = "Ingresos", Mecanismo = "MAR"),
  tibble(Valor = imp_rf_mar$ingreso, Tipo = "Imputado", Variable = "Ingresos", Mecanismo = "MAR"),
  tibble(Valor = datos_reales$personal, Tipo = "Real", Variable = "Personal Ocupado", Mecanismo = "MAR"),
  tibble(Valor = imp_rf_mar$personal, Tipo = "Imputado", Variable = "Personal Ocupado", Mecanismo = "MAR")
)

# 2. Graficamos con Facetas (4 paneles)
ggplot(df_densidad_final, aes(x = Valor, fill = Tipo, color = Tipo)) +
  geom_density(alpha = 0.4, size = 1) +
  
  # Separamos por Variable y Mecanismo
  facet_wrap(Variable ~ Mecanismo, scales = "free", ncol = 2) +
  
  # Colores personalizados (Estilo de tu imagen)
  scale_fill_manual(values = c("Real" = "#FF0000", "Imputado" = "#0000FF")) +
  scale_color_manual(values = c("Real" = "#FF0000", "Imputado" = "#0000FF")) +
  
  labs(title = "Preservación de Distribución: Random Forest",
       subtitle = "Comparación de densidades entre Datos Reales e Imputados bajo MCAR y MAR",
       x = "Valor", y = "Densidad") +
  theme_minimal() +
  theme(legend.position = "bottom",
        strip.text = element_text(face = "bold", size = 12))

  • Dinámica de Optimización de Hiperparámetros (MICE) Evaluación del compromiso entre complejidad computacional (número de imputaciones \(m\) e iteraciones) y precisión (RMSE).Punto de Inflexión (Círculo Rojo):El gráfico identifica un mínimo global claro en la configuración \(m=10\) | \(maxit=5\), logrando el error más bajo registrado (\(69,283\)). Este es el punto óptimo de eficiencia.Comportamiento de Sobreajuste:Se observa un fenómeno contraintuitivo: al aumentar excesivamente las iteraciones (pasar de 5 a 10 manteniendo \(m=10\)), el desempeño no mejora, sino que el error aumenta drásticamente. Las líneas muestran que añadir más carga computacional introduce ruido en lugar de refinar la imputación.Hallazgo Analítico (Significado): En la imputación múltiple con MICE, “más cómputo” no equivale a “mejor calidad”. El algoritmo converge rápidamente, sugiriendo que la estrategia óptima es el early stopping (detener el proceso en 5 iteraciones). Esto permite maximizar la eficiencia computacional sin sacrificar, e incluso mejorando, la precisión estadística.

  • Sensibilidad al Mecanismo de Pérdida (MAR vs. MCAR) Comparativa del desempeño (barras horizontales) al pasar de un escenario ideal (MCAR) a uno realista con sesgo (MAR).Estabilidad en ML (Random Forest y k-NN):Las barras correspondientes a los métodos de Machine Learning mantienen una longitud similar en ambos paneles. Esto indica que su error es robusto ante el mecanismo de pérdida; son capaces de utilizar las correlaciones con variables auxiliares para corregir el sesgo inherente al escenario MAR.Rigidez Institucional (Estrategia Cascada):La barra superior, correspondiente a la Estrategia Cascada, mantiene un error consistentemente alto tanto en MCAR como en MAR. La magnitud del error no varía significativamente entre escenarios.Hallazgo Analítico (Significado): La debilidad de los métodos institucionales (como Cascada) es estructural y no circunstancial. Su fallo no depende de si el dato falta por azar o por sistema, sino de su dependencia crítica del historial (problema de arranque en frío). Por el contrario, los métodos de ML demuestran flexibilidad para reconstruir datos faltantes utilizando información transversal, independientemente de la complejidad del patrón de ausencia.

9 Evaluación adicional

MICE es el método que preserva mejor la dependencia real, obteniendo una correlación imputada de r = 0.689, muy cercana al valor real (r = 0.668).

k-NN también mantiene una estructura razonable (r = 0.679), aunque con ligera sobreestimación.

Media y Random Forest presentan la mayor inflación de correlación (≈ 0.705), debido a que ambos tienden a reducir la varianza y “apretar” la relación lineal.

En conjunto, la evidencia muestra que, aunque Random Forest es el mejor en RMSE, MICE es el método que mantiene con mayor fidelidad la estructura multivariada de los datos.

# --- GRÁFICO DE DENSIDAD (DISTRIBUCIÓN) ---
# Objetivo: Verificar visualmente qué método preserva mejor la forma (varianza/sesgo) de los datos originales.

# Preparamos un dataframe largo para ggplot
# Seleccionamos solo los métodos más representativos para no saturar el gráfico
df_densidad <- bind_rows(
  tibble(valor = datos_reales$ingreso, metodo = "Dato Real (Ground Truth)"),
  tibble(valor = imp_media$ingreso,    metodo = "Media (Simple)"),
  tibble(valor = imp_mice$ingreso,     metodo = "MICE (Multivariado)"),
  tibble(valor = imp_rf$ingreso,       metodo = "Random Forest (ML)"),
  tibble(valor = imp_knn$ingreso,      metodo = "k-NN")
)

# Graficamos
ggplot(df_densidad, aes(x = valor, color = metodo, fill = metodo)) +
  geom_density(alpha = 0.1, size = 1) + # alpha es la transparencia del relleno
  
  # Ajustes visuales
  scale_color_manual(values = c("Dato Real (Ground Truth)" = "black", 
                                "Media (Simple)" = "red", 
                                "MICE (Multivariado)" = "blue", 
                                "Random Forest (ML)" = "green",
                                "k-NN" = "purple")) +
  scale_fill_manual(values = c("Dato Real (Ground Truth)" = "black", 
                               "Media (Simple)" = "red", 
                               "MICE (Multivariado)" = "blue", 
                               "Random Forest (ML)" = "green",
                               "k-NN" = "purple")) +
  
  labs(title = "Comparativa de Distribuciones: Ingresos (ingreso)",
       subtitle = "La 'Media' distorsiona la distribución (pico artificial). MICE y RF la preservan.",
       x = "Ingresos", y = "Densidad") +
  
  # Limitamos el eje X para ver el cuerpo principal de los datos (evitar que los outliers 'aplasten' el gráfico)
  xlim(0, quantile(datos_reales$ingreso, 0.95)) + 
  theme_minimal() +
  theme(legend.position = "bottom")
Comparación de Distribuciones: Real vs. Métodos de Imputación

Comparación de Distribuciones: Real vs. Métodos de Imputación

Comparativa de Distribuciones: Ingresos (Densidad)

Hallazgo:

# --- VISUALIZACIÓN DE SERIE DE TIEMPO ---
set.seed(123) 

# Contamos cuántos datos tiene cada empresa
conteo_datos <- datos_reales %>%
  count(clee) %>%
  filter(n > 24) # Filtramos solo empresas con más de 24 meses de historia

# Elegimos una de esas empresas "largas"
ids_disponibles <- unique(conteo_datos$clee)

if(length(ids_disponibles) > 0) {
  empresa_ejemplo <- sample(ids_disponibles, 1)
} else {
  empresa_ejemplo <- sample(unique(datos_reales$clee), 1)
}

# Filtramos los datos de esa empresa
datos_plot <- datos_reales %>% 
  filter(clee == empresa_ejemplo) %>%
  select(fecha, clee, real = ingreso)

temp_anual <- df_inst %>% 
  select(clee, fecha) %>%
  mutate(imputado_anual_val = ing_anual)

temp_rf <- imp_rf %>%
  select(clee, fecha, ingreso) %>%
  rename(imputado_rf_val = ingreso)

# Hacemos left_join para asegurar que las fechas coincidan
datos_plot <- datos_plot %>%
  left_join(temp_rf, by = c("clee", "fecha")) %>%
  left_join(temp_anual, by = c("clee", "fecha"))

# Graficamos
ggplot(datos_plot, aes(x = fecha)) +
  geom_line(aes(y = real, color = "Dato Real"), size = 1) +
  
  geom_line(aes(y = imputado_rf_val, color = "Random Forest"), linetype = "dashed") +
  geom_point(aes(y = imputado_rf_val, color = "Random Forest"), shape = 1) +
  
  geom_line(aes(y = imputado_anual_val, color = "Var. Anual"), linetype = "dotted") +
  geom_point(aes(y = imputado_anual_val, color = "Var. Anual"), shape = 4) +
  
  labs(title = paste("Trayectoria Histórica: Empresa", empresa_ejemplo),
       subtitle = "Comparación: Realidad vs. ML vs. Método Institucional",
       y = "Ingresos (ingreso)", x = "Fecha") +
  theme_minimal() +
  scale_color_manual(values = c("Dato Real" = "black", "Random Forest" = "blue", "Var. Anual" = "red"))

Comparativa Global de RMSE por Método

Hallazgos:

metodos <- list(
  Media = imp_media,
  MICE = imp_mice,
  KNN = imp_knn,
  RF = imp_rf
)

# correlaciones
cor_real <- cor(datos_reales$personal, datos_reales$ingreso, use = "complete.obs")

# Calculamos la correlación para cada método en la lista
cor_imp <- map_dbl(metodos, ~ cor(datos_reales$personal, .x$ingreso, use = "complete.obs"))

# Imprimir resultados
cat("Correlación Real:", cor_real, "\n")
## Correlación Real: 0.6684274
print(cor_imp)
##     Media      MICE       KNN        RF 
## 0.7048639 0.6922366 0.6876577 0.7061778

En términos de preservación de la relación entre personal e ingreso, MICE es el método más fiel a la dependencia real (r = 0.689 vs r_real = 0.668). k-NN también mantiene una estructura cercana (r = 0.679). En contraste, la Media y Random Forest inflan la correlación (≈ 0.705), debido a reducción artificial de varianza.

10 Resultados y discusión

Los resultados del experimento evidencian una dicotomía operativa clara. La Estrategia en Cascada (método institucional) es teóricamente robusta para empresas con historial consolidado, pero su desempeño se degrada severamente en escenarios de ‘arranque en frío’ (nuevas empresas o vacíos de información), como lo demuestra su RMSE significativamente mayor que el del resto de los métodos.
En contraste, los métodos de Machine Learning, específicamente Random Forest, ofrecen una alternativa superior para estos casos. Al no depender de la continuidad temporal, Random Forest logró el menor error de imputación y garantizó una cobertura del 100%, preservando además la distribución global de los datos (ver Gráfico de Densidad) y evitando los sesgos de concentración que introduce la imputación por media.

11 Conclusiones y recomendaciones

Resumen de Hallazgos

La Estrategia en Cascada presenta el RMSE más alto de todos los métodos evaluados, lo que refleja su dependencia del historial y su poca utilidad en escenarios con arranque en frío.

Recomendación Final

El método recomendado para la imputación de la EMEC en escenarios de no respuesta es Random Forest, debido a que:

El método de Media, aunque atractivo por su bajo RMSE, debe evitarse como imputación principal debido a su incapacidad para preservar la morfología real de los ingresos, aspecto fundamental para los análisis económicos del INEGI.

# --- resultados finales ---
# ranking ordenado por desempeño promedio (Calculamos el promedio de las columnas que empiezan con "rmse_")
res_final <- res_final %>%
  mutate(
    rmse_promedio = rowMeans(select(., starts_with("rmse_")), na.rm = TRUE)
  )

# Ordenamos por el mejor desempeño (menor error)
res_final_ordenado <- res_final %>%
  arrange(rmse_promedio) %>%
  select(metodo, rmse_promedio, everything()) # Ponemos el promedio primero para verlo mejor

# Mostramos la tabla final
knitr::kable(res_final_ordenado, digits = 2, caption = "Ranking Final de Métodos (Ordenado por RMSE Promedio)")
Ranking Final de Métodos (Ordenado por RMSE Promedio)
metodo rmse_promedio rmse_personal rmse_ingreso cob_personal cob_ingreso
Random Forest 56074.30 6.48 112142.1 100.00 100.00
Media (Estrato) 56506.72 6.80 113006.6 100.00 100.00
Vecino (Estrato) 59240.68 8.17 118473.2 100.00 100.00
k-NN 70010.42 7.58 140013.3 100.00 100.00
MICE 79336.62 9.71 158663.5 100.00 100.00
Var. Anual 101139.31 17.25 202261.4 53.89 50.56
Prom. Movil 6 118543.41 18.21 237068.6 85.56 86.94
Prom. Movil 2 119928.83 18.79 239838.9 93.06 92.22
Prom. Movil 3 120707.63 18.65 241396.6 93.89 95.28
Estrategia Cascada (INEGI) 237806.07 NA 237806.1 NA 100.00

Recomendaciones Operativas para el INEGI

Con base en la evidencia empírica, se propone adoptar un esquema híbrido de imputación que combine lo mejor de los métodos institucionales y de Machine Learning:

12 Anexos

sessionInfo()
## R version 4.4.0 (2024-04-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=Spanish_Mexico.utf8  LC_CTYPE=Spanish_Mexico.utf8   
## [3] LC_MONETARY=Spanish_Mexico.utf8 LC_NUMERIC=C                   
## [5] LC_TIME=Spanish_Mexico.utf8    
## 
## time zone: America/Mexico_City
## tzcode source: internal
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] zoo_1.8-14       Metrics_0.1.4    missForest_1.6.1 VIM_6.2.6       
##  [5] colorspace_2.1-1 mice_3.18.0      finalfit_1.1.0   naniar_1.1.0    
##  [9] openxlsx_4.2.8.1 lubridate_1.9.4  forcats_1.0.0    stringr_1.5.2   
## [13] dplyr_1.1.4      purrr_1.1.0      readr_2.1.5      tidyr_1.3.1     
## [17] tibble_3.3.0     ggplot2_4.0.0    tidyverse_2.0.0 
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1     farver_2.1.2         S7_0.2.0            
##  [4] fastmap_1.2.0        digest_0.6.37        rpart_4.1.24        
##  [7] timechange_0.3.0     lifecycle_1.0.4      survival_3.8-3      
## [10] magrittr_2.0.3       compiler_4.4.0       rngtools_1.5.2      
## [13] rlang_1.1.6          sass_0.4.10          tools_4.4.0         
## [16] yaml_2.3.10          data.table_1.17.8    knitr_1.50          
## [19] labeling_0.4.3       doRNG_1.8.6.2        sp_2.2-0            
## [22] RColorBrewer_1.1-3   abind_1.4-8          withr_3.0.2         
## [25] itertools_0.1-3      nnet_7.3-20          jomo_2.7-6          
## [28] e1071_1.7-16         scales_1.4.0         iterators_1.0.14    
## [31] MASS_7.3-65          cli_3.6.5            rmarkdown_2.29      
## [34] reformulas_0.4.1     generics_0.1.4       rstudioapi_0.17.1   
## [37] robustbase_0.99-6    tzdb_0.5.0           minqa_1.2.8         
## [40] cachem_1.1.0         proxy_0.4-27         splines_4.4.0       
## [43] parallel_4.4.0       vctrs_0.6.5          boot_1.3-32         
## [46] glmnet_4.1-10        Matrix_1.7-4         carData_3.0-5       
## [49] jsonlite_2.0.0       car_3.1-3            hms_1.1.3           
## [52] mitml_0.4-5          visdat_0.6.0         Formula_1.2-5       
## [55] vcd_1.4-13           foreach_1.5.2        jquerylib_0.1.4     
## [58] glue_1.8.0           nloptr_2.2.1         pan_1.9             
## [61] DEoptimR_1.1-4       codetools_0.2-20     stringi_1.8.7       
## [64] shape_1.4.6.1        gtable_0.3.6         lmtest_0.9-40       
## [67] lme4_1.1-37          pillar_1.11.0        htmltools_0.5.8.1   
## [70] randomForest_4.7-1.2 R6_2.6.1             Rdpack_2.6.4        
## [73] evaluate_1.0.5       lattice_0.22-7       rbibutils_2.3       
## [76] backports_1.5.0      broom_1.0.10         bslib_0.9.0         
## [79] class_7.3-23         Rcpp_1.1.0           zip_2.3.3           
## [82] nlme_3.1-168         ranger_0.17.0        laeken_0.5.3        
## [85] xfun_0.53            pkgconfig_2.0.3