Fortalecimiento del tratamiento de no respuesta en la EMEC: comparación de métodos de imputación

1 Síntesis

Breve resumen del proyecto y objetivo: evaluar y comparar métodos de imputación para la EMEC (variables clave: personal y m000a), incluyendo búsqueda de hiperparámetros y pruebas de preservación de correlaciones y coherencia temporal.

2 Introducción

• contexto: la EMEC y su tratamiento institucional de imputaciones (promedios móviles, variación porcentual, vecino más cercano ajustado por dominio–estrato). Enlace al glosario: https://www.inegi.org.mx/app/glosario/default.html?p=emec2018
• objetivo del documento: proponer y evaluar alternativas (media, MICE, K-NN, Random Forest), con optimización de hiperparámetros y evaluación mediante simulaciones controladas.

3 Paquetes y opciones

# paquetes necesarios
library(tidyverse)    # manipulación y graficado
library(openxlsx)     # lectura de xlsx
library(naniar)       # visualización faltantes
library(finalfit)     # prueba de little (mcar)
library(mice)         # imputación multiple
library(VIM)          # kNN imputacion (kNN)
library(missForest)  # random forest imputacion
library(Metrics)      # rmse
library(gridExtra)    # arreglos de graficas
library(bnstruct)     # redes bayesianas (opcional)

4 Carga y limpieza de datos

#install.packages("readxl")
archivo <- "DISEÑO_2018_INFORMACION_PROCESADA_SEPTIEMBRE_2025.xlsx"
datos_raw <- openxlsx::read.xlsx(archivo)
#glimpse(datos_raw)

# normalizar nombres en minúsculas
colnames(datos_raw) <- tolower(colnames(datos_raw))

# seleccionamos variables de interes y filtramos observados reales (imputado != 4)
datos <- datos_raw %>%
  select(personal, m000a, imputado, everything()) # mantenemos otras columnas por si sirven

# convertir a numeric (limpiando caracteres extra)
datos <- datos %>%
  mutate(
    personal = as.numeric(gsub("[^0-9\\.-]", "", personal)),
    m000a = as.numeric(gsub("[^0-9\\.-]", "", m000a))
  )

cat("registros totales:", nrow(datos_raw), "\n")

## registros totales: 6134

cat("registros con personal/m000a:", sum(!is.na(datos$personal) | !is.na(datos$m000a)), "\n")

## registros con personal/m000a: 6134

5 Exploración inicial

# resumen
datos %>% select(personal, m000a) %>% summary()

##     personal          m000a         
##  Min.   :     0   Min.   :       0  
##  1st Qu.:     7   1st Qu.:     677  
##  Median :    35   Median :    9548  
##  Mean   :   363   Mean   :  138776  
##  3rd Qu.:   129   3rd Qu.:   45537  
##  Max.   :200292   Max.   :64086687

# correlación entre personal y m000a
correlacion <- cor(datos$personal, datos$m000a, use = "complete.obs")
correlacion

## [1] 0.7308078

# scatter
ggplot(datos, aes(x = personal, y = m000a)) +
  geom_point(alpha = 0.4) +
  labs(title = paste0("personal vs m000a (r = ", round(correlacion, 2), ")"),
       x = "personal ocupado", y = "m000a") +
  theme_minimal()

6 Diseño del experimento (simulaciones)

• creamos un conjunto “ground truth” con los registros no imputados por el sistema (imputado != 4) para evaluar recuperación.
• simulamos mecanismos de no respuesta: MCAR (aleatorio) y MAR (dependiente de covariable simple como tamaño/estrato si está disponible). Aquí mostramos MCAR y un ejemplo simple de MAR.

datos_reales <- datos %>% filter(is.na(imputado) | imputado != 4) %>% select(personal, m000a)

# nos quedamos con filas completas (para simplificar la validación)
datos_reales <- datos_reales %>% drop_na()

# creamos copia para experimentos
n <- nrow(datos_reales)
cat("n registros reales (completos):", n, "\n")

## n registros reales (completos): 5910

6.1 Simulación MCAR

set.seed(2025)
prop_faltantes <- 0.20
# creamos indices a borrar aleatoriamente para cada variable (MCAR)
idx_personal_mcar <- sample(1:n, size = floor(prop_faltantes * n))
idx_m000a_mcar    <- sample(1:n, size = floor(prop_faltantes * n))

datos_mcar <- datos_reales
datos_mcar$personal[idx_personal_mcar] <- NA
datos_mcar$m000a[idx_m000a_mcar] <- NA

# verificar
vis_miss(datos_mcar)

6.2 Simulación MAR (ejemplo simple)

nota: si en tu dataset hay variables que expliquen prob. de no respuesta (ej. estrato, tamaño de empresa), úsalas para diseñar MAR.

# ejemplo: crear probabilidad de faltar dependiente del cuartil de 'personal'
cut_personal <- ntile(datos_reales$personal, 4)
prob_falta <- ifelse(cut_personal == 4, 0.05, ifelse(cut_personal == 3, 0.15, 0.30))
# para simular, sacamos NA en m000a con esa prob
set.seed(2025)
datos_mar <- datos_reales
rand <- runif(n)
datos_mar$m000a[rand < prob_falta] <- NA
# añadimos faltantes en personal de forma independiente
rand2 <- runif(n)
datos_mar$personal[rand2 < 0.20] <- NA

vis_miss(datos_mar)

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

Se implementan los métodos: - media (baseline) - MICE (pmm) con optimización de maxit e m - k-NN (VIM::kNN) con selección de k - Random Forest (missForest) como método no paramétrico

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 6.1.1 MICE: búsqueda de maxit y m (iteraciones e imputaciones)

# tomamos el dataset MCAR para optimizar (es controlado)
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) {
    cat("probando maxit", it, "m=", mv, "\n")
    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$m000a, comp$m000a, mascara_na$m000a)
    res_mice <- bind_rows(res_mice, tibble(maxit = it, m = mv, rmse_personal = rmse_p, rmse_m000a = rmse_m, rmse_promedio = mean(c(rmse_p, rmse_m), na.rm = TRUE)))
  }
}

## probando maxit 5 m= 5 
## probando maxit 5 m= 10 
## probando maxit 10 m= 5 
## probando maxit 10 m= 10 
## probando maxit 20 m= 5 
## probando maxit 20 m= 10

res_mice

## # A tibble: 6 × 5
##   maxit     m rmse_personal rmse_m000a rmse_promedio
##   <dbl> <dbl>         <dbl>      <dbl>         <dbl>
## 1     5     5         5171.    341311.       173241.
## 2     5    10         5022.    419081.       212052.
## 3    10     5         3592.   1060843.       532217.
## 4    10    10         5953.    423886.       214919.
## 5    20     5         6267.    721464.       363866.
## 6    20    10         2770.    306605.       154687.

# gráfica resumen
ggplot(res_mice, aes(x = factor(maxit), y = rmse_promedio, group = factor(m), color = factor(m))) +
  geom_point() + geom_line() +
  labs(title = "optimización mice: maxit vs m", x = "maxit", y = "rmse promedio") +
  theme_minimal()

7.1.2 k-NN: búsqueda de k

k_values <- c(3, 5, 10, 20)
res_knn <- tibble()

for(k in k_values) {
  cat("probando k =", k, "\n")
  imp_knn <- kNN(df_opt, variable = c("personal", "m000a"), k = k, imp_var = FALSE)
  rmse_p <- calcular_rmse(datos_reales$personal, imp_knn$personal, mascara_na$personal)
  rmse_m <- calcular_rmse(datos_reales$m000a, imp_knn$m000a, mascara_na$m000a)
  res_knn <- bind_rows(res_knn, tibble(k = k, rmse_personal = rmse_p, rmse_m000a = rmse_m, rmse_promedio = mean(c(rmse_p, rmse_m), na.rm = TRUE)))
}

## probando k = 3 
##    m000a    m000a 
##        0 64086687 
## personal personal 
##        0   122183 
## probando k = 5 
##    m000a    m000a 
##        0 64086687 
## personal personal 
##        0   122183 
## probando k = 10 
##    m000a    m000a 
##        0 64086687 
## personal personal 
##        0   122183 
## probando k = 20 
##    m000a    m000a 
##        0 64086687 
## personal personal 
##        0   122183

res_knn

## # A tibble: 4 × 4
##       k rmse_personal rmse_m000a rmse_promedio
##   <dbl>         <dbl>      <dbl>         <dbl>
## 1     3         5765.    301030.       153397.
## 2     5         5761.    228625.       117193.
## 3    10         4657.    212803.       108730.
## 4    20         5357.    210051.       107704.

ggplot(res_knn, aes(x = k, y = rmse_promedio)) + geom_line() + geom_point() +
  labs(title = "optimización k-NN: selección de k", x = "k", y = "rmse promedio") + theme_minimal()

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

# elegimos mejores hiperparámetros según resultados
mejor_mice <- res_mice %>% arrange(rmse_promedio) %>% slice(1)
mejor_knn  <- res_knn %>% arrange(rmse_promedio) %>% slice(1)
mejor_mice; mejor_knn

## # A tibble: 1 × 5
##   maxit     m rmse_personal rmse_m000a rmse_promedio
##   <dbl> <dbl>         <dbl>      <dbl>         <dbl>
## 1    20    10         2770.    306605.       154687.

## # A tibble: 1 × 4
##       k rmse_personal rmse_m000a rmse_promedio
##   <dbl>         <dbl>      <dbl>         <dbl>
## 1    20         5357.    210051.       107704.

# aplicamos métodos al dataset MCAR (experimento principal)
df_base <- datos_mcar
mascara_na <- as.data.frame(is.na(df_base))

# 1) media simple
imp_media <- df_base
for(col in c("personal","m000a")) imp_media[[col]][is.na(imp_media[[col]])] <- mean(df_base[[col]], na.rm = TRUE)

# 2) MICE con mejores parametros
imp_mice_mod <- mice(df_base, method = 'pmm', m = as.integer(mejor_mice$m[1]), maxit = as.integer(mejor_mice$maxit[1]), printFlag = FALSE, seed = 123)
imp_mice <- mice::complete(imp_mice_mod)

# 3) kNN con mejor k
k_sel <- as.integer(mejor_knn$k[1])
imp_knn <- kNN(df_base, variable = c("personal","m000a"), k = k_sel, imp_var = FALSE)

##    m000a    m000a 
##        0 64086687 
## personal personal 
##        0   122183

# 4) Random Forest
imp_rf <- missForest(df_base, verbose = FALSE)
imp_rf <- imp_rf$ximp

# calcular rmse por metodo
metodos <- list(media = imp_media, mice = imp_mice, knn = imp_knn, rf = imp_rf)
res_final <- tibble()
for(nombre in names(metodos)){
  imp <- metodos[[nombre]]
  rmse_p <- calcular_rmse(datos_reales$personal, imp$personal, mascara_na$personal)
  rmse_m <- calcular_rmse(datos_reales$m000a, imp$m000a, mascara_na$m000a)
  res_final <- bind_rows(res_final, tibble(metodo = nombre, rmse_personal = rmse_p, rmse_m000a = rmse_m, rmse_promedio = mean(c(rmse_p, rmse_m), na.rm = TRUE)))
}
res_final

## # A tibble: 4 × 4
##   metodo rmse_personal rmse_m000a rmse_promedio
##   <chr>          <dbl>      <dbl>         <dbl>
## 1 media          5944.    317144.       161544.
## 2 mice           2770.    306605.       154687.
## 3 knn            5357.    210051.       107704.
## 4 rf             4722.    222455.       113588.

En los resultados del experimento MCAR, los métodos no presentan un comportamiento homogéneo entre las dos variables analizadas. Para personal, el error es relativamente más bajo en MICE (≈ 2770) comparado con k-NN y Random Forest. Esto tiene sentido porque personal es una variable con menor dispersión y un comportamiento más lineal, lo cual se ajusta bien al mecanismo de imputación por regresiones iterativas de MICE.

En cambio, para m000a (ingresos) el patrón cambia: k-NN y Random Forest obtienen los menores errores (≈ 210 mil–222 mil) mientras que MICE y media simple quedan rezagados. La variable de ingresos presenta una distribución altamente asimétrica, con valores extremos y relaciones no lineales respecto al personal. En este contexto, métodos basados en proximidad (k-NN) y árboles (Random Forest) capturan mejor la estructura compleja de la variable.

# grafica comparativa
res_final_long <- res_final %>% pivot_longer(cols = starts_with("rmse_"), names_to = "variable", values_to = "rmse")

ggplot(res_final_long, aes(x = metodo, y = rmse, fill = variable)) +
  geom_col(position = "dodge") +
  labs(title = "comparativa final de rmse por metodo", y = "rmse") + theme_minimal()

Considerando el RMSE promedio, el método con mejor desempeño global es k-NN (k = 20), seguido de cerca por Random Forest. Esto indica que, cuando la relación entre unidades es más importante que la estructura paramétrica de los datos, la imputación por vecinos cercanos puede recuperar mejor los valores verdaderos.

No obstante, MICE conserva una ventaja relevante: mantiene adecuadamente la correlación entre personal e ingresos. De hecho, el coeficiente imputado por MICE (0.724) es casi idéntico al real (0.725), lo cual es deseable para análisis estructurales y modelado. mejor la estructura compleja de la variable.

9 Evaluación adicional

• preservación de correlaciones entre personal y m000a
• comparación de distribuciones (densidad)
• coherencia temporal (si se cuenta con series históricas por unidad)

# correlaciones
cor_real <- cor(datos_reales$personal, datos_reales$m000a, use = "complete.obs")
cor_imp <- map_dbl(metodos, ~ cor(datos_reales$personal, .x$m000a, use = "complete.obs"))
cor_real; cor_imp

## [1] 0.7255031

##     media      mice       knn        rf 
## 0.7207267 0.7241018 0.7254506 0.7265438

En términos de preservación de dependencia entre variables, todos los métodos funcionan razonablemente bien, pero Random Forest y k-NN presentan una ligera tendencia a incrementar artificialmente la correlación, debido a que suavizan valores extremos de ingresos.

MICE, por el contrario, mantiene la correlación casi intacta, lo que sugiere que es preferible cuando se realizan análisis de “estructura económica” o cálculos que dependen del vínculo entre empleo e ingresos, como productividad o ingresos por trabajador.

# densidad m000a (comparativa)
df_plot <- bind_rows(
  tibble(valor = datos_reales$m000a, origen = "real"),
  tibble(valor = imp_media$m000a, origen = "media"),
  tibble(valor = imp_mice$m000a, origen = "mice"),
  tibble(valor = imp_knn$m000a, origen = "knn"),
  tibble(valor = imp_rf$m000a, origen = "rf")
)

p <- ggplot(df_plot, aes(x = valor, color = origen)) + geom_density() + xlim(0, quantile(datos_reales$m000a, 0.95)) + theme_minimal()
print(p)

Al comparar densidades, observamos que la imputación por media introduce un sesgo sustancial, eliminando variabilidad en ingresos. MICE mantiene bien la forma general, aunque tiende a subestimar ligeramente la cola superior.

k-NN y Random Forest logran preservarla mejor, lo cual se vuelve especialmente relevante dado que la distribución real de ingresos es fuertemente sesgada hacia la derecha. Su capacidad de imitar colas largas sugiere que estos métodos son más robustos cuando se trata de imputar montos económicos con alta heterogeneidad.

10 Resultados y discusión

En esta sección se muestra la tabla res_final, que resume el RMSE obtenido por cada método en las dos variables de interés: personal e ingresos (m000a). Los resultados indican que ningún método domina en ambas variables, sino que el desempeño depende de la estructura estadística:

Para personal, el método con menor RMSE fue generalmente MICE, lo cual es consistente con la naturaleza menos dispersa y más lineal de esta variable. La imputación mediante regresiones iterativas (pmm) aprovecha la correlación existente entre unidades y preserva adecuadamente sus variaciones.

Para ingresos (m000a), los métodos con mejor desempeño fueron k-NN (k = 20) y Random Forest, ya que capturan relaciones no lineales y manejan mejor la presencia de valores extremos típicos de las variables económicas monetarias. Estos métodos acercan los valores imputados a la distribución real, especialmente en la cola superior.

Además, se evaluó la preservación de la correlación entre personal e ingresos, observando que MICE mantiene prácticamente intacta la dependencia estadística entre ambas variables, mientras que k-NN y RF tienden a suavizar los valores grandes y por lo tanto elevan ligeramente la correlación.

Desde el punto de vista de la coherencia institucional, la EMEC requiere mantener continuidad histórica por empresa, especialmente en unidades de menor tamaño (estrato 1). Por ello se justifica el uso de los métodos actuales (promedios móviles, variaciones porcentuales y vecino por estrato), ya que preservan patrones propios de cada unidad. Sin embargo, los resultados muestran que los métodos multivariados pueden mejorar la precisión microestadística.

Por esta razón se recomienda un esquema híbrido:

Utilizar procedimientos institucionales (promedios móviles, variación anual, vecino ajustado por estrato) cuando la prioridad sea mantener la trayectoria histórica de la unidad económica.

Emplear métodos multivariados avanzados (MICE, Random Forest, k-NN) cuando el objetivo sea mejorar la calidad de estimaciones agregadas, capturar relaciones entre variables y reducir sesgo en series sujetas a heterogeneidad alta, como los ingresos.

Finalmente, el análisis de la distribución muestra que la imputación por media reduce artificialmente la variabilidad, mientras que MICE subestima levemente la cola superior y k-NN/RF la preservan mejor.

11 Conclusiones y recomendaciones

1. La presente versión del documento incorpora la búsqueda de hiperparámetros solicitada por el profesor (elección de k en k-NN y selección de maxit y m en MICE), asegurando que los métodos comparados operen en su configuración óptima.

2. De acuerdo con la tabla res_final, no existe un método dominante para todas las variables. Para ingresos (m000a), los métodos con menor error fueron k-NN y Random Forest, mientras que para personal el método con mejor desempeño fue MICE, especialmente en términos de preservación de correlaciones.

3. Se recomienda una política combinada:

   • Utilizar métodos institucionales cuando se requiera mantener la coherencia temporal por unidad económica.
     
   • Usar métodos multivariados (MICE, RF, k-NN) como alternativa para mejorar la precisión microestadística y la calidad de estimaciones agregadas.

4. Se aconseja incluir de manera rutinaria diagnósticos de no respuesta, tales como la prueba de Little para MCAR, exploración de mecanismos MAR, análisis de correlaciones posteriores a la imputación y evaluación de sensibilidad.

5. Las variables económicas con alta asimetría —como ingresos (m000a)— se benefician particularmente de métodos no paramétricos como Random Forest y k-NN, los cuales capturan mejor las colas largas y valores atípicos sin asumir linealidad.

res_final %>%
arrange(rmse_promedio)

## # A tibble: 4 × 4
##   metodo rmse_personal rmse_m000a rmse_promedio
##   <chr>          <dbl>      <dbl>         <dbl>
## 1 knn            5357.    210051.       107704.
## 2 rf             4722.    222455.       113588.
## 3 mice           2770.    306605.       154687.
## 4 media          5944.    317144.       161544.

A partir de los valores presentados en la tabla res_final, el ranking final según el RMSE promedio fue:

1. k-NN (k = 20)
   
2. Random Forest
   
3. MICE
   
4. Media

Este orden confirma que los métodos no paramétricos (k-NN y Random Forest) resultan más efectivos para imputar ingresos (m000a), ya que capturan relaciones no lineales y manejan adecuadamente valores extremos propios de variables económicas. Por otro lado, MICE continúa siendo una alternativa sólida para la variable personal, especialmente cuando se busca preservar relaciones estadísticas como la correlación entre empleo e ingresos. La imputación por media, como era esperable, presenta el peor desempeño debido a su incapacidad para reproducir la variabilidad real del fenómeno.

12 Anexos

• reproducibilidad: semilla, versiones de paquetes
• notas sobre ejecución en servidor (memoria para missForest y mice en datos grandes)

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] bnstruct_1.0.15  igraph_2.2.1     bitops_1.0-9     gridExtra_2.3   
##  [5] Metrics_0.1.4    missForest_1.6.1 VIM_6.2.6        colorspace_2.1-1
##  [9] mice_3.18.0      finalfit_1.1.0   naniar_1.1.0     openxlsx_4.2.8.1
## [13] lubridate_1.9.4  forcats_1.0.0    stringr_1.5.2    dplyr_1.1.4     
## [17] purrr_1.1.0      readr_2.1.5      tidyr_1.3.1      tibble_3.3.0    
## [21] ggplot2_4.0.0    tidyverse_2.0.0 
## 
## loaded via a namespace (and not attached):
##  [1] Rdpack_2.6.4         rlang_1.1.6          magrittr_2.0.3      
##  [4] e1071_1.7-16         compiler_4.4.0       vctrs_0.6.5         
##  [7] pkgconfig_2.0.3      shape_1.4.6.1        fastmap_1.2.0       
## [10] backports_1.5.0      labeling_0.4.3       utf8_1.2.6          
## [13] rmarkdown_2.29       tzdb_0.5.0           nloptr_2.2.1        
## [16] itertools_0.1-3      visdat_0.6.0         xfun_0.53           
## [19] glmnet_4.1-10        jomo_2.7-6           randomForest_4.7-1.2
## [22] cachem_1.1.0         jsonlite_2.0.0       pan_1.9             
## [25] broom_1.0.10         parallel_4.4.0       R6_2.6.1            
## [28] bslib_0.9.0          stringi_1.8.7        vcd_1.4-13          
## [31] RColorBrewer_1.1-3   ranger_0.17.0        car_3.1-3           
## [34] boot_1.3-32          rpart_4.1.24         lmtest_0.9-40       
## [37] jquerylib_0.1.4      Rcpp_1.1.0           iterators_1.0.14    
## [40] knitr_1.50           zoo_1.8-14           Matrix_1.7-4        
## [43] splines_4.4.0        nnet_7.3-20          timechange_0.3.0    
## [46] tidyselect_1.2.1     rstudioapi_0.17.1    abind_1.4-8         
## [49] yaml_2.3.10          codetools_0.2-20     doRNG_1.8.6.2       
## [52] lattice_0.22-7       withr_3.0.2          S7_0.2.0            
## [55] evaluate_1.0.5       survival_3.8-3       proxy_0.4-27        
## [58] zip_2.3.3            pillar_1.11.0        carData_3.0-5       
## [61] rngtools_1.5.2       foreach_1.5.2        reformulas_0.4.1    
## [64] generics_0.1.4       sp_2.2-0             hms_1.1.3           
## [67] scales_1.4.0         laeken_0.5.3         minqa_1.2.8         
## [70] class_7.3-23         glue_1.8.0           tools_4.4.0         
## [73] robustbase_0.99-6    data.table_1.17.8    lme4_1.1-37         
## [76] rbibutils_2.3        nlme_3.1-168         Formula_1.2-5       
## [79] cli_3.6.5            gtable_0.3.6         DEoptimR_1.1-4      
## [82] sass_0.4.10          digest_0.6.37        farver_2.1.2        
## [85] htmltools_0.5.8.1    lifecycle_1.0.4      mitml_0.4-5         
## [88] MASS_7.3-65