Imputación de ingresos en la Encuesta Permanente de Hogares en base a los datos del Sistema Integrado Previsional Argentino mediante modelos de machine learning
Trabajo final para la Especialización en Explotación de datos y descubrimiento del conocimiento de la Facultad de Cs. Exactas y Naturales de la Universidad de Buenos Aires
Autor: Facundo Lastra
Diciembre 2022
Resumen
El presente trabajo busca contribuir al estudio del nivel de salarios en Argentina, a partir de realizar un ejercicio de imputación de ingresos en la Encuesta Permanente de Hogares (EPH-INDEC), utilizando como fuente la Muestra Longitudinal de Empleo Registrado (MLER) del Sistema Integrado Previsional Argentino (SIPA). Este ejercicio explora una posible corrección a los problemas de subdeclaración de ingresos de la EPH, buscando otorgar a esta encuesta de información más fehaciente sobre los ingresos laborales de la población asalariada registrada. El trabajo se organiza de la siguiente manera. En primer lugar, introducimos el problema, presentamos las fuentes de datos y realizamos una revisión de la literatura que abordó los problemas de subdeclaración y captación de ingresos de las encuestas de hogares. Luego definimos las variables que utilizamos, compatibilizamos la información de ambas fuentes y realizamos un análisis exploratorio de los datos. En tercer lugar, ajustamos distintos modelos de machine learning para la predicción e imputación de ingresos, y evaluamos los resultados obtenidos. Por último, establecemos las conclusiones, recomendaciones y líneas de investigación que se desprenden de este trabajo.
Palabras clave: EPH - Mercado de Trabajo - Imputación de ingresos - SIPA - Machine learning
1.Introducción
#Cargamos paquetes utilizados
library(tidyverse)
library(Hmisc)
library(kableExtra)
library(statar)
library(eph)
library(ggplot2)
library(broom)
library(caret)
library(doParallel)
library(Metrics)
library(lightgbm)
library(formattable)
library(openxlsx)1.1. Motivación del problema
Los estadísticas socio-laborales se construyen a través de procesos que incluyen etapas de recolección, compilación, procesamiento y publicación de información. En todas estas etapas se toman decisiones metodológicas para intentar que los datos expresen de manera más representativa las situaciones concretas que se quieren describir. Los ingresos laborales, por su propia naturaleza, resultan un aspecto de las estadísticas socio-laborales muy difícil de cuantificar y estimar, principalmente por problemas para la recolección de los datos. Por su parte, las encuestas de hogares son la herramienta más difundida para estudiar las relaciones laborales y una fuente importante para la construcción de estadísticas sobre los ingresos de los hogares. Este tipo de encuestas ofrecen un amplio rango de datos demográficos, que permiten caracterizar las condiciones de vida y empleo de la población. Son, además, un instrumento que permite captar información sobre un amplio universo que incluye no sólo a la población asalariada, sino también al cuentapropismo y otras formas atípicas de empleo. Sin embargo, estas encuestas presentan problemas para la captación de los casos con ingresos más altos.
Existen otras fuentes de datos sobre ingresos que han sido usadas de forma complementaria a la EPH. En particular, se destacan los datos del Sistema Integrado Previsional Argentino (SIPA), que provienen de los registros administrativos de la seguridad social. Si bien esta fuente de información no está exenta de posibles problemas de subdeclaración o no captación de ingresos, sus datos son más fehacientes que los de la EPH, dado que están basados en declaraciones juradas que las empresas realizan todos los meses ante las autoridades nacionales para determinar los aportes y contribuciones al sistema de seguridad social de sus empleados/as. En este marco, el objetivo del presente trabajo es proponer una metodología para la imputación de ingresos en la Encuesta Permanente de Hogares, basada en el uso de herramientas de aprendizaje automático o machine learning sobre datos del SIPA. De esta manera exploramos una posible corrección a los problemas de subdeclaración, buscando otorgar a la EPH de información más fehaciente sobre los ingresos laborales de la población asalariada registrada. Por el tipo de información que utilizamos para construir los modelos de imputación, este ejercicio es de carácter exploratorio y, por ello, hacia el final del presente trabajo reflexionamos sobre la posibilidad de perfeccionar la imputación si se accede a información más detallada que los datos de SIPA disponibles para uso público.
1.2. Presentación de la información, fuentes de datos y su alcance
Como decíamos anteriormente, la EPH es una herramienta muy difundida para el estudio de las condiciones laborales y las características sociodemográficas de la población. Esta encuesta cubre actualmente a 31 aglomerados del todo país, donde habitan aproximadamente el 62% de la población. Tiene como objetivo principal relevar información sobre las características sociodemográficas y económicas de la población, permitiendo realizar investigaciones sobre la estructura social, a partir de datos sobre individuos y hogares. La EPH es una encuesta por muestreo que se realiza mediante operativos de relevamiento y se basa en las respuestas declaradas por los encuestados/as. En el marco del presente trabajo, centrado en el estudio de los ingresos laborales, vale la pena remarcar que la EPH es una fuente ampliamente utilizada para el estudio de los ingresos de la población asalariada, cuentapropista, como también de las transferencias del estado (jubilaciones, pensiones y otros subsidios o planes sociales). En este trabajo utilizaremos la ‘versión continua’ de esta encuesta, la cual provee información trimestral homogénea desde el año 2003 hasta la actualidad.
tabla1 <- read.xlsx("tablas/tabla1.xlsx")
colnames(tabla1)[1] <- ""
tabla1 %>%
kable(caption = "<b>Tabla 1. </b> Características de la EPH y la MLER") %>%
kable_styling(bootstrap_options = c("striped"), full_width = F, position = "float_right")| EPH | MLER | |
|---|---|---|
| Origen de la información | Encuestas a hogares | Declaraciones juradas de empresas |
| Tipo de fuente | Encuesta por muestreo | Muestreo aleatorio simple al azar |
| Información sobre región | Por aglomerados urbanos | Por provincia |
| Criterio de para definir ubicación | Domicilio de residencia | Lugar de radicación de la empresa |
| Problemas de captación de ingresos | Mayor | Menor |
La otra fuente de datos en la que se centra nuestro trabajo es el SIPA, que forma parte del Sistema de Seguridad Social Nacional (SSSN). El SIPA fue creado a fines del año 2008 (Ley 26.425) y es un régimen previsional público que se financia por un sistema solidario de reparto. El sistema de seguridad social argentino también está formado por otros subsistemas de menor tamaño que no forman parte del SIPA y que quedan por fuera de su alcance. Entre ellos se encuentran algunos sistemas previsionales administrados por provincias y municipios, las cajas previsionales de bancos y asociaciones profesionales, subsistemas de las fuerzas de seguridad, como también otros subsistemas más pequeños. De esta manera el SIPA incluye a todas las personas físicas mayores de 18 años de edad que trabajen en relación de dependencia en la actividad pública o privada, excluyendo a quienes estén comprendidos en los subsistemas recién mencionados (Sánchez, Pacífico y Kennedy, 2016; ANSES, 2011).
La información del SIPA a nivel de microdatos no se encuentra disponible de forma directa y pública para su uso estadístico. El Observatorio de Empleo y Dinámica Empresarial (OEDE) del Ministerio de Trabajo, Empleo y Seguridad Social (MTEySS) publicó bases de microdatos con fuente en los registros administrativos del sistema de seguridad social, adaptados para su uso estadístico mediante la Muestra Longitudinal de Empleo Registrado (MLER). Esta base de datos sigue estándares internacionales de diseminación de información de registros administrativos, que incluyen la anonimización de personas y empresas, y que contiene información de más de 500 mil trabajadores/as. El instrumento de recolección de datos en este caso son las declaraciones juradas que las empresas deben presentar ante la AFIP a fin de cada mes con la nómina de empleados y sus remuneraciones. La MLER es una muestra representativa de las personas que tuvieron una relación laboral asalariada registrada entre 1996 y 2015 (OEDE, 2018). Este es el último año con información disponible, siendo ello un límite muy importante de esta fuente de información, que no es actualizada de forma sistemática.
tabla2 <- read.xlsx("tablas/tabla2.xlsx")
colnames(tabla2)[1] <- ""
tabla2 %>%
kable(caption = "<b>Tabla 2. </b> Universo de análisis en la EPH y la MLER") %>%
kable_styling(bootstrap_options = c("striped"), full_width = F, position = "float_left")| EPH | MLER | |
|---|---|---|
| Alcance | Aglomerados Urbanos | Total país |
| Ramas primarias | Baja participación | Alta participación |
| Empleo privado registrado | Sí | Sí |
| Empleo privado no registrado | Sí | No |
| Cuentapropismo | Sí | No |
| Trabajo doméstico remunerado | Sí | No |
| Sector público | Sí | Parcial |
| Patrones | Sí | No |
| Trabajadores/as familiares | Sí | No |
La base de datos de la MLER está disponible para el total del país y para todas las ramas de actividad económica privada, excluyendo a los subsistemas mencionados anteriormente. Quedan excluidos de esta base una gran porción de los/as empleados/as del sector público y la totalidad de personas ocupadas en el trabajo doméstico remunerado (la categoría ‘servicio doméstico’). Asimismo, por el origen mismo de los datos provenientes de los registros públicos, quedan por fuera de la base de datos los vínculos no registrados o ‘en negro’ (ya que no cotizan en la seguridad social), los/as cuentapropistas, patrones y trabajadores/as familiares. Como se puede ver en la Tabla 2 el universo de análisis en la MLER es mucho más restringido en cuanto a las categorías ocupacionales que cubre en comparación con la EPH. Mientras la primera cubre sólo al empleo privado registrado, la EPH alcanza también al resto de las categorías ocupacionales, aunque con foco en las actividades que se desarrollan en los aglomerados urbanos.
1.3. El problema de subcaptación de ingresos y su posible corrección
Si bien son un instrumento de gran relevancia para la generación de estadísticas públicas, las encuestas de hogares tienen límites para captar los ingresos laborales de los perceptores de rentas, las ganancias empresariales y los cuentapropistas calificados (Roca y Pena, 2001). Aunque a un menor nivel, estos límites también se verifican para la población asalariada, y particularmente para los deciles de más altos ingresos. Ello sucede por tres motivos principales: la no respuesta de ingresos a la encuesta, la subdeclaración de ingresos y la subcaptación de perceptores (Gómez Sabaíni y Rossignolo, 2014). Estos problemas parecen haber sido profundizados desde inicios de los años 2000, dado que la brecha entre el promedio salarial que arroja la EPH y el estimado a través del SIPA ha ido en aumento a favor de este último (Sánchez, Pacífico y Kennedy, 2016). Por ello es plausible suponer que, a pesar de los posibles problemas de subcaptación que pudieran tener los datos del SIPA, éstos poseen información de mejor calidad que la EPH. Por ello sería razonable realizar algún tipo de corrección en la EPH utilizando la información proveniente de los registros públicos de la seguridad social, para así mejorar la calidad de la información de dicha encuesta.
Desde el INDEC se realizaron esfuerzos para intentar corregir los problemas de captación de ingresos en la EPH. Hasta el año 2015, la no respuesta de ingresos fue tratada por el método hot deck, imputando valores mediante un algoritmo que busca registros similares, los cuales se utilizan como ‘donantes’ del valor a imputar. Desde el 2016 en adelante, los microdatos se corrigen asignando ponderación igual a cero a los casos sin respuesta y reasignando la ponderación a otro hogar. También la EPH incorporó la corrección de ingresos para los cuentapropistas y empleadores que realizan tareas técnicas y profesionales, pero sin modificar los ingresos declarados por los asalariados registrados. Sin embargo, el aumento de la brecha con respecto a los datos de SIPA muestra que esas correcciones no lograron revertir el problema de fondo, en particular para la población asalariada y registrada.
A pesar de que existe un consenso generalizado sobre la existencia de problemas de subcaptación de ingresos en las encuestas de hogares, aún no se aplican técnicas estadísticas novedosas como el machine learning para intentar suplir aunque sea parcialmente este problema. Esta falta de innovación por parte de los institutos de estadística públicos se verifica en distintas áreas de vacancia en las estadísticas demográficas y socio-laborales. Un ejemplo de ello es el trabajo realizado por Rosati (2021), quien aplica distintos modelos de árboles de decisión y redes neuronales para realizar ejercicios de imputación de datos de ingresos sin respuesta en la EPH. El autor demuestra que todos los algoritmos entrenados tienen una mejor performance que el método de hot deck utilizado actualmente por el organismo que construye los datos.
El presente trabajo busca realizar un aporte en un sentido similar, proponiendo una metodología de corrección de datos que podría aplicarse en la construcción de las estadísticas públicas del mercado de trabajo. Así es que a continuación presentamos una aplicación de modelos de regresión basados en árboles de decisión para la imputación de ingresos en la EPH, aplicando las técnicas de aprendizaje supervisado que proliferaron en los últimos años. En el próximo apartado seguimos analizando el problema de la brecha existente entre los ingresos estimados por la EPH y los datos provenientes del SIPA en base al procesamiento de estas dos fuentes, para exponer luego la aplicación de estos métodos en el tercer apartado.
2. Creación de variables y análisis exploratorio de los datos
Realizamos nuestro análisis y procesamiento de los datos para el período 2003-2014. Este recorte temporal se debe a que el 2003 es el primer año con datos disponibles en la versión de la EPH que estamos utilizando, mientras que el último año con datos completos en la MLER es el 2014. Dadas estas restricciones, el ejercicio de imputación de ingresos se realizará para este último año. En el análisis exploratorio para el período 2003-2014, utilizamos los microdatos correspondientes a la onda del 3er trimestre en la EPH y los datos correspondientes al mes de julio en la MLER.
###Preprocesamiento MLER ###
#Cargo base cruda
MLER_raw <- readRDS("Fuentes/mler.rds")
#Con este loop definimos los meses y años que quiero tomar
vector <- c()
for (anio in 2003:2014) {
for (mes in 7) {
vector <- c(vector, paste0("rt", anio, "0", mes))}
}
variables_MLER <- c("fnac_anu", "ide_trabajador", "letra","pondera", "provi", "r32", "r34",
"relacion", "tr_0414", "tr_0404", "sexo", vector)
MLER <- MLER_raw %>%
select(variables_MLER)
remove(MLER_raw)
MLER <- MLER %>%
gather(.,
key = Periodo_remunera,
value= Remun_bruta,
12:ncol(.)) %>%
filter(Remun_bruta>0) %>%
mutate(ANO4 = as.numeric(substr(Periodo_remunera,3,6)),
Mes = as.numeric(substr(Periodo_remunera,7,8)),
TRIMESTRE = case_when(Mes %in% c(1:3) ~1,
Mes %in% c(4:6) ~2,
Mes %in% c(7:9) ~3,
Mes %in% c(10:12) ~4),
ide_trabajador = as.character(ide_trabajador))
saveRDS(MLER, file="Fuentes/mler_2003_2014.rds")
###Preprocesamiento EPH###
#Descarga de bases con el paquete 'eph'
EPHs_3ro <- get_microdata(year = 2003:2015, trimester = 3, type = "individual") # (3er trimestre, excepto 2007 que usamos 2do)
EPHs_2do <- get_microdata(year = 2007, trimester = 2, type = "individual")
EPH <- bind_rows(EPHs_2do, EPHs_3ro)
variables_EPH <- c("P21", "CAT_OCUP", "ESTADO", "PP07H", "PP07I",
"PP04B1", "PP04A", "REGION", "ANO4", "PONDERA",
"REGION", "CH04", "CH06", "PP04C", "PP04C99",
"PP04B_COD", "PP04B_CAES")
EPH <- EPH %>%
select(variables_EPH)
EPH$P21[EPH$P21<0|is.na(EPH$P21)]<-0
saveRDS(EPH, file="Fuentes/eph_2003_2014.rds")2.1. Definición del universo de análisis y de la variable de ingresos
Como indicamos en la presentación de las fuentes utilizadas, la EPH y la MLER poseen información sobre universos distintos y utilizan diferentes mecanismos de recolección. Ello hace necesario realizar un procesamiento de los datos para compatibilizar el universo de análisis con el que trabajamos. Por un lado, filtramos los datos de EPH, tomando información sobre la ocupación principal de la población asalariada registrada del sector privado, excluyendo la categoría de ‘servicio doméstico’. Por otro lado, utilizamos los datos de la MLER sobre la ocupación principal de las personas (suponiendo que su ocupación principal es la de ingreso máximo) y les descontamos los aportes patronales. A continuación resumimos los criterios adoptados:
- Criterios EPH:
- Asalariados del sector privado (CAT_OCUP==3, PP04A!=1)
- Ramas distintas a servicio doméstico (PP04B1!=1)
- Registro: Descuento jubilatorio en las variables PP07H o PP07I
- Ingresos de la ocupación principal (P21)
- Criterios MLER :
- Ingresos: Se tomará el máximo ingreso en casos de más de una relación laboral en simultáneo
- Ingresos Brutos: Se multiplican por 0.83 para descontar los aportes personales
- En ambos casos se toman como válidos todos los ingresos netos mayores a $9
#Cargamos bases preprocesadas
EPH_preprocesada <- readRDS("Fuentes/eph_2003_2014.rds")
MLER_preprocesada <- readRDS("Fuentes/mler_2003_2014.rds")
#Aplicamos filtros
EPH <- EPH_preprocesada %>%
filter(CAT_OCUP==3, ESTADO == 1, PP07H==1 | PP07I==1,
PP04B1!=1, PP04A!=1) %>%
mutate(Remun_neta = P21)%>%
mutate(Base = "EPH") %>%
filter(P21>9)
MLER <- MLER_preprocesada %>%
mutate(Remun_neta = Remun_bruta*0.83) %>%
filter(Remun_neta>9) %>%
group_by(Periodo_remunera,ide_trabajador) %>%
filter(Remun_neta==max(Remun_neta)) %>%
mutate(PONDERA = pondera,
Base = "MLER")
#Armamos tabla sobre múltiple ocupación en la MLER
MLER_multiple_ocup <- MLER %>%
group_by(ANO4) %>%
summarise(Relaciones_laborales = sum(pondera),
Trabajadores = length(unique(ide_trabajador))*33.3,
Cociente = Trabajadores/Relaciones_laborales)
colnames(MLER_multiple_ocup) <- c("Año", "Relaciones Laborales", "Trabajadores/as", "Cociente")
MLER_multiple_ocup %>%
kable(caption = "<b>Tabla 3. </b> Cantidad de trabajadores y de relaciones laborales. MLER 2003-2014") %>%
kable_styling(bootstrap_options = c("striped"), full_width = F, position = "float_right")| Año | Relaciones Laborales | Trabajadores/as | Cociente |
|---|---|---|---|
| 2003 | 3532067 | 3527069 | 0.9985852 |
| 2004 | 3993467 | 3988008 | 0.9986331 |
| 2005 | 4439900 | 4434028 | 0.9986775 |
| 2006 | 4869133 | 4862799 | 0.9986991 |
| 2007 | 5299600 | 5292769 | 0.9987110 |
| 2008 | 5647800 | 5640221 | 0.9986581 |
| 2009 | 5498933 | 5491570 | 0.9986609 |
| 2010 | 5633366 | 5625869 | 0.9986690 |
| 2011 | 5889533 | 5881479 | 0.9986325 |
| 2012 | 5954366 | 5946248 | 0.9986365 |
| 2013 | 5987933 | 5979981 | 0.9986719 |
| 2014 | 5989033 | 5981013 | 0.9986609 |
Un problema que se podría presentar a la hora de comparar los datos de la MLER y la EPH es que la primera presenta datos a nivel de relación laboral, mientras que la segunda lo hace a nivel de personas. Como consecuencia de ello, la MLER capta todas las relaciones laborales registradas de una persona, mientras que la EPH sólo contiene información detallada sobre la ocupación principal (en la que se trabaja más horas) e información con menos detalle de la ocupación secundaria. Para evaluar la magnitud de casos con múltiple ocupación en la MLER, observamos la relación entre la cantidad de relaciones laborales captadas por la muestra y la cantidad de trabajadores.
En la Tabla 3 se observa que son pocos los casos de múltiple ocupación. Para todos los años estudiados, la cantidad de trabajadores representa más del 99% de las relaciones laborales. Esto quiere decir que la relación es casi de uno a uno, por lo que son pocas las relaciones laborales captadas por la MLER que corresponden a una misma persona en un mismo periodo. Al ser tan pocos los casos de múltiple ocupación y al tener información heterogénea entre la EPH y la MLER sobre las ocupaciones secundarias, tomamos la decisión de centrar el análisis en la ocupación principal de las personas. En el caso de la EPH utilizamos la distinción entre la ocupación principal y secundaria propuesta por dicha encuesta, mientras que en la MLER consideramos como ocupación principal a aquella que le reporte mayores ingresos a la persona.
2.2. Creación de variables explicativas y homogeneización de la información
Además de la información sobre ingresos de la ocupación principal, las variables que utilizaremos como predictoras para el ejercicio de imputación son: la edad, la edad al cuadrado, el sexo, la región, la rama de actividad, el tamaño del establecimiento y el percentil de ingresos. Para ello creamos las variables en base a la información de las bases y/o realizamos una homogeneización a partir de la información original de ambas fuentes.
# Generamos variables sobre Ramas de actividad
## Cargamos informacion sobre homogeneizacion de las clasificaciones sobre ramas de actividad según como lo indicamos en la Figura 1 más adelante en este trabajo
rama_cod_cod_MLER <- read.xlsx("Fuentes/clasificacion_ramas.xlsx",
sheet = "rama_cod_cod_MLER")
rama_cod_cod_EPH2003_2010 <- read.xlsx("Fuentes/clasificacion_ramas.xlsx",
sheet = "rama_cod_cod_EPH2003_2010")
rama_cod_cod_EPH2011_2014 <- read.xlsx("Fuentes/clasificacion_ramas.xlsx",
sheet = "rama_cod_cod_EPH2011_2014")
rama_cod_descripcion <- read.xlsx("Fuentes/clasificacion_ramas.xlsx",
sheet = "rama_cod_descripcion")
##Organizamos etiquetas de variable CAES utilizando la función del paquete eph
EPH <- organize_caes(EPH)
##Generamos variable rama
EPH2003_2010 <- EPH %>%
filter(ANO4 %in% 2003:2010) %>%
left_join(rama_cod_cod_EPH2003_2010, by=c("caes_seccion_cod"="cod_EPH2003_2010")) %>%
left_join(rama_cod_descripcion, by="rama_cod")
EPH2011_2014 <- EPH %>%
filter(ANO4 %in% 2011:2014) %>%
left_join(rama_cod_cod_EPH2003_2010, by=c("caes_seccion_cod"="cod_EPH2003_2010")) %>%
left_join(rama_cod_descripcion, by="rama_cod")
EPH <- bind_rows(EPH2003_2010, EPH2011_2014)
remove(EPH2003_2010, EPH2011_2014)
MLER <- MLER %>% left_join(rama_cod_cod_MLER, by=c("letra"="cod_mler")) %>%
left_join(rama_cod_descripcion, by="rama_cod")
#Generación y homogeneización del resto de variables
MLER <- MLER %>%
ungroup() %>%
mutate(
tramo_tamanio = case_when(
tr_0414=="Hasta 9" ~ "tramo 1",
tr_0414=="de 10 a 49" ~ "tramo 2",
tr_0414=="de 50 a 200" ~ "tramo 3",
tr_0414=="más de 200" ~ "tramo 4"),
edad=2014-fnac_anu,
edad2=edad*edad,
provi=case_when(
provi=="Neuquén" ~ "Neuquén",
provi=="Cordoba" ~ "Córdoba" ,
provi=="Entre RÃos" ~ "Entre Ríos",
provi=="RÃo Negro" ~ "Río Negro",
provi=="Santa Fe" ~ "Santa Fé",
provi=="Tucumán" ~ "Tucumán",
TRUE ~ as.character(provi)),
ingreso=Remun_neta) %>%
mutate(
sexo= factor(sexo,
levels=c("Mujer", "Hombre")),
tramo_tamanio= factor(tramo_tamanio,
levels= c("tramo 1","tramo 2","tramo 3","tramo 4")),
region= factor(
case_when(
provi %in% c("Capital Federal", "Gran Buenos Aires") ~ "AMBA",
provi %in% c("Santa Fé", "Córdoba", "Buenos Aires", "La Pampa", "Entre Ríos") ~ "Centro",
provi %in% c("Tucumán", "La Rioja" , "Salta", "Catamarca", "Santiago del Estero", "Jujuy") ~ "NOA",
provi %in% c("Formosa", "Corrientes", "Chaco", "Misiones") ~ "NEA",
provi %in% c("San Juan", "San Luis", "Mendoza" ) ~ "Cuyo",
provi %in% c("Chubut", "Tierra del Fuego", "Santa Cruz", "Neuquén", "Río Negro") ~ "Patagonia"),
levels=c("AMBA", "Centro", "NEA", "Cuyo", "Patagonia", "NOA"))) %>%
group_by(ANO4) %>%
mutate(
percentil=ntile(Remun_neta, 100)) %>%
ungroup()
EPH <- EPH %>%
ungroup() %>%
mutate(
tramo_tamanio = case_when(
PP04C%in%c(1:6) ~ "tramo 1", #Hasta 10 inclusive
PP04C%in%c(7:8) ~ "tramo 2", #De 11 a 40
PP04C%in%c(9:10) ~ "tramo 3", #De 41 a 200
PP04C%in%c(11:12) ~ "tramo 4"),
edad=CH06,
edad2=edad*edad,
sexo= case_when(
CH04 ==1 ~ "Hombre",
CH04 ==2 ~ "Mujer"),
region= factor(
case_when(
REGION ==01 ~ "AMBA",
REGION ==43 ~ "Centro",
REGION ==40 ~ "NOA",
REGION ==41 ~ "NEA",
REGION ==42 ~ "Cuyo",
REGION ==44 ~ "Patagonia"),
levels=c("AMBA", "Centro", "NEA", "Cuyo", "Patagonia", "NOA")),
ingreso=Remun_neta) %>%
mutate(
sexo= factor(sexo,
levels=c("Mujer", "Hombre")),
tramo_tamanio= factor(tramo_tamanio,
levels= c("tramo 1","tramo 2","tramo 3","tramo 4"))) %>%
group_by(ANO4) %>%
mutate(
percentil=ntile(Remun_neta, 100)) %>%
ungroup()
saveRDS(EPH, "Fuentes/dataset_EPH.RDS")
saveRDS(MLER, "Fuentes/dataset_MLER.RDS")En el caso de la edad, utilizamos la edad declarada por el/la encuestado/a en la EPH, mientras que en la MLER la estimamos restando la fecha de nacimiento del/la empleado/a por el año en el que fue declarada la relación laboral. Como se suele hacer en estimaciones de ecuaciones de Mincer, incluimos la edad al cuadrado como una variable que podría aportar información sobre el plus que reciben los/as trabajadores/as por su experiencia laboral, dando valores aún más altos a los casos de mayor edad.
El tamaño de establecimiento es organizado en 4 tramos según cantidad de empleados: de 1 a 10 empleados, 11 a 40, 41 a 200 y más de 200 para la EPH; y de 1 a 9, 10 a 49, 50 a 200 y más de 200 para la MLER. Como se puede observar, los tramos no coinciden perfectamente debido a la forma en que son publicados los microdatos en ambas bases, pero igualmente guardan un alto grado de similitud.
También estimamos el percentil de ingresos, el cual aporta información sobre en qué parte de la distribución se encuentra cada caso tomando como referencia los demás casos de la misma base. La inclusión del percentil como variable explicativa permite conservar información sobre la posición de cada caso en la estructura de la distribución de ingresos en la EPH y aprovechar dicha información para la imputación de datos desde la MLER.
Por último, la información sobre región y rama de actividad se construye homogeneizando la información de las dos fuentes utilizadas, con criterios que detallamos en el análisis exploratorio de los datos del próximo subapartado.
2.3. Cantidad de casos y distribución de las variables
Luego de filtrar y transformar las bases de microdatos según los criterios mencionados, procedemos a analizar las características de éstas en los años 2003-2014. En primer lugar, observamos cuál es la evolución de la cantidad absoluta de casos y de los casos ponderados para cada una de las bases durante este período.
comp_casos_EPH <- EPH %>%
group_by(ANO4) %>%
summarise(Casos=n(),
Casos_ponderados=sum(PONDERA))
comp_casos_MLER <- MLER %>%
group_by(ANO4) %>%
summarise(Casos=n(),
Casos_ponderados=sum(PONDERA))
comp_casos <- comp_casos_EPH %>%
left_join(comp_casos_MLER, by="ANO4") %>%
mutate(ANO4=as.character(ANO4))
colnames(comp_casos) <- c("Año", "Casos", "Casos Ponderados", "Casos", "Casos ponderados")
comp_casos %>%
kable(caption = "<b>Tabla 4. </b> Cantidad de casos y de casos ponderados. EPH y MLER 2003-2014",
format.args = list(big.mark = ".")) %>%
kable_styling(bootstrap_options = c("striped"), full_width = F, position = "float_left") %>%
add_header_above(., c(" ", "EPH" = 2, "MLER" = 2)) | Año | Casos | Casos Ponderados | Casos | Casos ponderados |
|---|---|---|---|---|
| 2003 | 3.061 | 1.778.945 | 105.962 | 3.532.067 |
| 2004 | 3.689 | 2.017.652 | 119.804 | 3.993.467 |
| 2006 | 6.260 | 2.665.078 | 146.074 | 4.869.133 |
| 2007 | 7.285 | 3.371.943 | 158.988 | 5.299.600 |
| 2008 | 7.338 | 3.544.833 | 169.434 | 5.647.800 |
| 2009 | 7.330 | 3.640.148 | 164.968 | 5.498.933 |
| 2010 | 7.395 | 3.752.123 | 169.001 | 5.633.366 |
| 2011 | 7.442 | 3.947.640 | 176.686 | 5.889.533 |
| 2013 | 6.831 | 3.787.188 | 179.638 | 5.987.933 |
| 2014 | 7.518 | 3.840.748 | 179.671 | 5.989.033 |
En la Tabla 4 se observa la gran diferencia que existe entre la EPH y la MLER en cuanto al tamaño de la muestra y la población sobre que ésta se proyecta. En el inicio del período estudiado, el dataset preprocesado de la EPH cuenta con algo más de 3.000 casos relevados para el año 2003, correspondientes a la muestra de los registros de la EPH puntual para el 3er trimestre de ese año. Ponderados, estos casos representan a aproximadamente 1,78 millones de trabajadores/as. Los datos de la MLER, en cambio, están compuestos por más de 300.000 casos que se corresponden a la muestra del 3% de la población asalariada registrada del mes de julio de cada año. Estos casos ponderados representan algo más de 3,53 millones relaciones registradas con los criterios adoptados en nuestro trabajo.
También se evidencia en los datos el fuerte crecimiento del empleo privado registrado que ocurrió durante el período estudiado, especialmente durante los años 2003-2008. Según el dato de la EPH se pasó de 1,78 millones personas asalariadas en 2003 a 3,54 millones en 2008 y, según la información de la MLER, el empleo registrado privado pasó de ser 3,53 millones a 5,64 millones para el mismo período. Ambas fuentes también coinciden en mostrar el posterior estancamiento en la creación de empleo privado registrado que se verificó desde el 2009 hasta el fin del período analizado.
Estos grandes contrastes en la cantidad de casos se deben a las diferencias que tiene una muestra basada en encuestas como lo es la EPH, en comparación con una muestra basada en registros administrativos como la MLER. Estas diferencias se basan, principal aunque no exclusivamente, en que la primera sólo nos puede aportar información para el empleo urbano de 31 aglomerados del país, mientras que la MLER abarca a las relaciones laborales declaradas por las empresas de todo el país. Detengámonos ahora en evaluar si estas diferencias en cuanto al alcance de las encuestas se traducen en una dispar representatividad de las regiones del territorio nacional.
comp_region_EPH <- EPH %>%
filter(ANO4 %in% c(2003, 2014)) %>%
filter(!is.na(region)) %>%
group_by(ANO4, region) %>%
summarise(Casos=n(),
Casos_ponderados=sum(PONDERA)) %>%
group_by(ANO4) %>%
mutate(Porcentaje=Casos_ponderados/sum(Casos_ponderados),
base="EPH") %>%
select(ANO4, region, Porcentaje, base)
comp_region_MLER <- MLER %>%
filter(ANO4 %in% c(2003, 2014)) %>%
filter(!is.na(region)) %>%
group_by(ANO4, region) %>%
summarise(Casos=n(),
Casos_ponderados=sum(PONDERA)) %>%
group_by(ANO4) %>%
mutate(Porcentaje=Casos_ponderados/sum(Casos_ponderados),
base="MLER") %>%
select(ANO4, region, Porcentaje, base)
comp_region<- bind_rows(comp_region_EPH, comp_region_MLER) %>%
pivot_wider(names_from = c("base", "ANO4"), values_from = Porcentaje) %>%
mutate(across(2:ncol(.), ~percent(.)))
colnames(comp_region) <- c("Región", "2003", "2014", "2003", "2014")
comp_region %>%
kable(caption = "<b>Tabla 5. </b> Porcentaje de casos ponderados según región. EPH y MLER. Años 2003 y 2014") %>%
kable_styling(bootstrap_options = c("striped"), full_width = F, position = "float_right") %>%
add_header_above(., c(" ", "EPH" = 2, "MLER" = 2)) | Región | 2003 | 2014 | 2003 | 2014 |
|---|---|---|---|---|
| AMBA | 65.02% | 59.73% | 46.67% | 44.30% |
| Centro | 19.14% | 21.56% | 31.49% | 32.27% |
| NEA | 2.02% | 2.81% | 3.72% | 4.39% |
| Cuyo | 5.16% | 5.80% | 5.73% | 5.91% |
| Patagonia | 3.05% | 4.12% | 5.52% | 6.19% |
| NOA | 5.60% | 5.98% | 6.87% | 6.94% |
Al respecto es importante destacar que la información sobre la provincia o región se define según distintos criterios para cada una de las bases: en la EPH los datos sobre el aglomerado de pertenencia del/la trabajador/a se construyen en función de la dirección donde se encuentra su hogar, mientras que en la MLER el dato sobre la provincia de la relación laboral se basa en el domicilio de radicación de las empresas. En la Tabla 5 mostramos el porcentaje de casos ponderados sobre el total, según las cinco regiones en que clasificamos las bases, a partir de homogeneizar los datos de la provincia de radicación de las empresas en la MLER con la información sobre los aglomerados urbanos donde residen los hogares en la EPH.
Como allí podemos observar, la estructura de los datos según región es distinta para cada base, pero igualmente guardan un importante grado de similitud. CABA concentra el 18% de los casos ponderados en la EPH y cerca del 26% en la MLER. La zona Centro (que incluye a Santa Fé, Córdoba, Provincia de Buenos Aires, La Pampa y Entre Ríos) significa el 63% de los casos ponderados en la EPH y poco menos del 51% en la MLER. Las otras cuatro zonas muestran guarismos menores al 7% para cada una de las dos bases.
Otra dimensión importante a tener en cuenta es el alcance en la captación de casos según ramas de actividad. Para la construcción de una variable para la rama de actividad de los puestos de trabajo fue necesario homogeneizar la información disponible sobre rama de actividad dentro de la EPH y, a su vez, entre EPH y MLER. Para ello proponemos el siguiente agrupamiento de ramas según las etiquetas disponibles en cada base.
Figura 1. Agrupamiento propuesto para las etiquetas sobre rama de actividad. EPH y MLER.
La Figura 1 muestra la organización de las etiquetas sobre ramas y su codificación, que homogeniza información de la MLER y la EPH. Como se observa en la Figura anterior, la información original de las bases sobre las ramas de actividad está organizada de distintas maneras en la MLER y la EPH. Asimismo, dentro de la EPH coexisten dos sistemas de clasificación, uno para el periodo 2003-2010 y otro para los años 2011-2014, que tienen pequeñas diferencias entre sí (Las etiquetas y códigos sombreados corresponden al período 2011-2014 y los no sombreados a 2003-2010).
En la composición del empleo por ramas sí se notan diferencias importantes entre la MLER y la EPH, que pueden explicarse por el distinto alcance que cada una de las bases de datos tienen. Como la EPH tiene un alcance restringido a los aglomerados urbanos de mayor tamaño, las actividades primarias agrupadas en la categoría “Agricultura, pesca, ganadería, …” poseen sólo un peso de entre el 1,82% para el 2003 y el 0,71% para el 2015 en los datos de esta encuesta. En cambio, en la MLER su porcentaje sobre el empleo total es del 8,33% y 6,67% para cada año respectivamente. Como contraparte, en términos generales se verifica que las ramas de actividad que podrían considerárselas como más relacionadas al espacio urbano tienen una mayor preponderancia en la EPH que en la MLER.
comp_rama_EPH <- EPH %>%
filter(ANO4 %in% c(2003, 2014)) %>%
filter(!is.na(rama)) %>%
group_by(ANO4, rama) %>%
summarise(Casos=n(),
Casos_ponderados=sum(PONDERA)) %>%
group_by(ANO4) %>%
mutate(Porcentaje=Casos_ponderados/sum(Casos_ponderados),
base="EPH") %>%
select(ANO4, rama, Porcentaje, base)
comp_rama_MLER <- MLER %>%
filter(ANO4 %in% c(2003, 2014)) %>%
filter(!is.na(rama)) %>%
group_by(ANO4, rama) %>%
summarise(Casos=n(),
Casos_ponderados=sum(PONDERA)) %>%
group_by(ANO4) %>%
mutate(Porcentaje=Casos_ponderados/sum(Casos_ponderados),
base="MLER") %>%
select(ANO4, rama, Porcentaje, base)
comp_rama<- bind_rows(comp_rama_EPH, comp_rama_MLER) %>%
pivot_wider(names_from = c("base", "ANO4"), values_from = Porcentaje) %>%
mutate(across(2:ncol(.), ~percent(.)))
colnames(comp_rama) <- c("Rama de actividad", "2003", "2014", "2003", "2014")
comp_rama %>%
kable(caption = "<b>Tabla 6. </b> Porcentaje de casos ponderados según rama de actividad. EPH y MLER. Años 2003 y 2014") %>%
kable_styling(bootstrap_options = c("striped"), full_width = F) %>%
add_header_above(., c(" ", "EPH" = 2, "MLER" = 2)) | Rama de actividad | 2003 | 2014 | 2003 | 2014 |
|---|---|---|---|---|
| Actividades inmobiliarias, empresariales, de alquiler y administrativas | 11.67% | 4.96% | 14.27% | 13.86% |
| Agricultura, pesca, ganadería, caza, silvicultura, minas y canteras | 1.82% | 29.54% | 8.33% | 6.67% |
| Comercio y reparaciones | 19.37% | 22.94% | 16.31% | 18.35% |
| Construcción | 2.71% | 5.36% | 4.19% | 6.88% |
| Enseñanza, actividades profesionales, científicas y técnicas | 8.69% | 3.72% | 7.08% | 6.12% |
| Hoteles y restaurantes | 3.56% | 9.65% | 3.21% | 4.02% |
| Industrias manufactureras | 21.58% | 1.15% | 22.26% | 20.68% |
| Intermediación financiera | 3.01% | 5.63% | 2.93% | 2.63% |
| Servicios comunitarios, sociales, personales y de entretenimiento | 7.30% | 1.23% | 6.72% | 6.03% |
| Servicios sociales y de salud | 6.96% | 8.83% | 4.49% | 4.64% |
| Suministro de electricidad, gas, agua y afines | 2.03% | 0.92% | 1.27% | 1.10% |
| Transporte, almacenamiento y comunicación | 11.31% | 6.08% | 8.95% | 9.00% |
Como conclusión del análisis exploratorio que realizamos en este apartado, podemos decir que la EPH y la MLER muestran grandes diferencias, aunque todas ellas explicables por los diferentes mecanismos de recolección de datos y el alcance de las muestras. Asimismo, law composiciones de las muestras ponderadas no parecen variar significativamente en cuanto a la representatividad de las 5 zonas del país, aunque sí hay diferencias en las composición por rama de actividad, dado que las actividades primarias son mejor captadas por la MLER que la EPH.
2.4. Comparación de ingresos
Veamos ahora cómo se distribuyen los ingresos en los registros de la EPH y la MLER, de forma tal que podamos entender cuáles son los fundamentos de la brecha existente entre las dos fuentes. A continuación graficamos la distribución de casos sin ponderar de las remuneraciones netas para el año 2003.
EPH_MLER <- EPH %>%
select(ANO4,Remun_neta,PONDERA,Base) %>%
bind_rows(MLER %>% select(ANO4,Remun_neta,PONDERA,Base))
Resumen_ingresos <- EPH_MLER %>%
group_by(ANO4,Base) %>%
summarise(Casos = n(),
Casos_expand = sum(PONDERA),
Ingresos_promedio = weighted.mean(Remun_neta,PONDERA),
Ingresos_mediana = median(Remun_neta))
EPH_MLER2003 <- EPH_MLER %>%
filter(ANO4 %in% c(2003))
mediaEPH2003 <- mean(EPH_MLER2003$Remun_neta[EPH_MLER2003$Base=="EPH"])
mediaMLER2003 <- mean(EPH_MLER2003$Remun_neta[EPH_MLER2003$Base=="MLER"])
ggplot(EPH_MLER2003,
aes(x=Remun_neta,group = Base, fill = Base))+
geom_density(alpha = 0.5)+
scale_x_continuous(limits = c(0,5000))+
labs(title = "Gráfico 1. Distribución de Ingresos según base. Año 2003",
subtitle = "Ingresos entre 0 y 5000 pesos. Casos muestrales")+
theme(title = element_text(size = 12),
legend.position = "bottom",
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.title.x = element_blank(),
legend.key.size = unit(x = 2,units = "mm"))+
geom_vline(xintercept = mediaEPH2003 ,
linetype='dashed', color='red', size=0.8) +
geom_vline(xintercept = mediaMLER2003,
linetype='dashed', color='blue', size=0.8) +
annotate("text",x= mediaEPH2003*1.65, y=0.001,
label= paste0("Media EPH: $", round(mediaEPH2003)), color='red') +
annotate("text",x= mediaMLER2003*1.65, y=0.0015,
label= paste0("Media MLER: $", round(mediaMLER2003)), color='blue')
Las líneas verticales punteadas azul y roja indican el promedio de la remuneración neta para la MLER y la EPH respectivamente. Como queda representado en el gráfico, el promedio de la MLER es mayor al promedio de la EPH, ilustrando el sesgo hacia remuneraciones menores que tiene la EPH. En cuanto a la distribución de las remuneraciones, la MLER muestra una mayor concentración de casos en las colas de la distribución, es decir, en las remuneraciones más bajas y en las más altas. No obstante, el sesgo hacia las remuneraciones más altas más que compensa a las remuneraciones bajas, llevando el promedio de la MLER por sobre el de la EPH.
Por otra parte, se pueden ver distintos máximos locales hacia la derecha de la distribución de la EPH. Estos se podrían explicar por la tendencia hacia el ‘redondeo’ en las respuestas de los/as encuestados/as, lo que puede interpretarse como una muestra de la inexactitud que tienen los datos basados en la respuesta autodeclarada en las encuestas de hogares. En cambio, la cola derecha de la distribución de la MLER desciende más gradualmente, ya que la información de esta fuente se corresponde con el monto exacto de la liquidación de los salarios de los/as empleados/as incluidos/as en la muestra.
EPH_MLER2014 <- EPH_MLER %>%
filter(ANO4 %in% c(2014))
mediaEPH2014 <- mean(EPH_MLER2014$Remun_neta[EPH_MLER2014$Base=="EPH"])
mediaMLER2014 <- mean(EPH_MLER2014$Remun_neta[EPH_MLER2014$Base=="MLER"])
ggplot(EPH_MLER2014,
aes(x=Remun_neta,group = Base, fill = Base))+
geom_density(alpha = 0.5)+
scale_x_continuous(limits = c(0,40000))+
labs(title = "Gráfico 2. Distribución de Ingresos según base. Año 2014",
subtitle = "Ingresos entre 0 y 40000 pesos. Casos muestrales")+
theme(title = element_text(size = 12),
legend.position = "bottom",
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.title.x = element_blank(),
legend.key.size = unit(x = 2,units = "mm"))+
geom_vline(xintercept = mediaEPH2014 ,
linetype='dashed', color='red', size=0.8) +
geom_vline(xintercept = mediaMLER2014,
linetype='dashed', color='blue', size=0.8) +
annotate("text",x= mediaEPH2014*1.65, y=0.00014 ,
label= paste0("Media EPH: $", round(mediaEPH2014)), color='red') +
annotate("text",x= mediaMLER2014*1.64, y=0.00012 ,
label= paste0("Media MLER: $", round(mediaMLER2014)), color='blue')
En la otra punta del período que estamos estudiando, el año 2014, los patrones que describimos anteriormente se encuentran aún más exacerbados. En el Gráfico 2 se nota aún más la mayor concentración de casos en las colas para la MLER. A su vez, el promedio de la MLER se aleja más del promedio de la EPH en 2014, en comparación con la distancia entre promedios registrada para el 2003. También se verifican, y de forma aún más acentuada que en el gráfico anterior, los máximos locales hacia la derecha de la distribución de las remuneraciones netas reportadas por la EPH, que da cuenta del sesgo hacia el ‘redondeo’ de los ingresos declarados.
Para tener una idea sobre cómo fue evolucionando la brecha entre los datos reportados por la EPH y por la MLER, observemos a continuación cuál es la evolución de los promedios de los ingresos laborales para cada una de las bases a lo largo de todo el período.
Resumen_ingresos <- Resumen_ingresos %>%
mutate(ANO4=as.integer(ANO4))
ggplot(Resumen_ingresos, aes(y=Ingresos_promedio, x=ANO4, col=Base)) +
geom_line(size=1.5) +
theme_minimal() +
scale_x_continuous(labels=as.character(Resumen_ingresos$ANO4),breaks=Resumen_ingresos$ANO4) +
labs(title = "Gráfico 3. Remuneración neta promedio de la ocupación principal",
subtitle = "Expresada en términos nominales. EPH y MLER 2003-2014") +
theme(title = element_text(size = 12),
legend.position = "bottom",
axis.title.y = element_blank(),
axis.title.x = element_blank())
En el gráfico se observa cómo la remuneración neta promedio para cada una de las bases se fue distanciando a lo largo del tiempo. La línea azul, que representa el promedio de las remuneraciones netas en la MLER, siempre está por encima del mismo dato calculado por la EPH. En el período 2003-2007, la distancia entre las dos líneas se fue ensanchando, pero es a partir del año 2008 que el distanciamiento se fue haciendo cada vez más pronunciado. Esto puede ser tomado como una muestra de que los problemas de captación de ingresos de la EPH se fueron profundizando a lo largo del tiempo.
La mayor diferencia entre promedios se expresa también en la brecha de ingresos que se verifica entre las dos bases de datos estudiadas, tal como lo ilustra el Gráfico 4.
#Brecha en procentaje
eph_ingresos <- Resumen_ingresos %>%
select(ANO4, Base, Ingresos_promedio) %>%
filter(Base=="EPH") %>%
mutate(Ingresos_promedio_eph=Ingresos_promedio) %>%
select(ANO4, Ingresos_promedio_eph)
mler_ingresos <- Resumen_ingresos %>%
select(ANO4, Base, Ingresos_promedio) %>%
filter(Base=="MLER") %>%
mutate(Ingresos_promedio_mler=Ingresos_promedio) %>%
select(ANO4, Ingresos_promedio_mler)
comp_brecha_ingresos <- eph_ingresos %>% left_join(mler_ingresos, by="ANO4") %>%
mutate(brecha=Ingresos_promedio_mler/Ingresos_promedio_eph)
ggplot(comp_brecha_ingresos, aes(y=brecha, x=ANO4)) +
geom_line(color = "light blue", size=2 ) +
theme_minimal() +
scale_x_continuous(labels=as.character(Resumen_ingresos$ANO4),breaks=Resumen_ingresos$ANO4) +
scale_y_continuous(limits=c(0.5, 1.5)) +
labs(title = "Gráfico 4. Brecha de la remuneración neta promedio de la ocupación principal",
subtitle = "MLER/EPH 2003-2014") +
theme(title = element_text(size = 12),
legend.position = "bottom",
axis.title.y = element_blank(),
axis.title.x = element_blank())
Como ya lo adelantamos al comienzo de este trabajo, la brecha promedio de los ingresos mantuvo una tendencia hacia el alza en todo el período. En 2003 la brecha fue de tan sólo 1,03, luego comenzó a subir hasta llegar a 1,12 en 2007 y terminó hacia el final del periodo a un valor de 1,27 para 2014. Este derrotero ilustra cómo los problemas de subcaptación de ingresos en la EPH se han ido profundizado en los 11 años con datos.
Resulta interesante saber si esta brecha en la captación de ingresos se distribuye de forma igual en todas las ramas de actividad. Con el reagrupamiento de categorías que describimos en el sub-apartado anterior, pasamos ahora analizar el promedio de ingresos para cada rama y su brecha para los dos extremos del período analizado, los años 2003 y 2014.
comp_ramas <- EPH %>% filter(ANO4 %in% c(2003, 2014)) %>%
select(ANO4,Remun_neta, rama, PONDERA, Base) %>%
bind_rows(MLER %>% ungroup() %>% filter(ANO4 %in% c(2003, 2014)) %>%
select(ANO4,Remun_neta, rama, PONDERA, Base))
Resumen_ramas <- comp_ramas %>%
group_by(ANO4, Base, rama) %>%
summarise(Casos = n(),
Casos_expand = sum(PONDERA),
Ingresos_promedio = weighted.mean(Remun_neta,PONDERA))
Resumen_ramas_comparativo <- Resumen_ramas %>%
select(ANO4, Base, rama, Ingresos_promedio) %>%
pivot_wider(names_from = c(ANO4, Base), values_from = Ingresos_promedio,
names_prefix = "ingresos_") %>%
filter(!is.na(rama)) %>%
mutate(brecha2003= ingresos_2003_MLER / ingresos_2003_EPH,
brecha2014= ingresos_2014_MLER / ingresos_2014_EPH) %>%
arrange(brecha2014)
colnames(Resumen_ramas_comparativo) <- c("Rama", "EPH 2003", "MLER 2003", "EPH 2014", "MLER 2014", "Brecha 2003", "Brecha 2014")
Resumen_ramas_comparativo %>%
kable(digits=2, caption = "<b>Tabla 6.</b> Promedios de remuneraciones netas y brecha entre bases según ramas de actividad. EPH y MLER. Años 2003 y 2014.") %>%
kable_styling("striped") %>%
column_spec(column = 1, width = "3in")| Rama | EPH 2003 | MLER 2003 | EPH 2014 | MLER 2014 | Brecha 2003 | Brecha 2014 |
|---|---|---|---|---|---|---|
| Hoteles y restaurantes | 563.09 | 508.89 | 8181.36 | 5391.62 | 0.90 | 0.66 |
| Enseñanza, actividades profesionales, científicas y técnicas | 600.70 | 516.06 | 7387.42 | 5730.71 | 0.86 | 0.78 |
| Actividades inmobiliarias, empresariales, de alquiler y administrativas | 836.91 | 717.67 | 10225.99 | 8170.88 | 0.86 | 0.80 |
| Industrias manufactureras | 860.46 | 1012.99 | 11576.94 | 11074.10 | 1.18 | 0.96 |
| Servicios comunitarios, sociales, personales y de entretenimiento | 853.69 | 761.50 | 8835.35 | 8931.85 | 0.89 | 1.01 |
| Construcción | 845.36 | 581.88 | 6751.39 | 7332.11 | 0.69 | 1.09 |
| Comercio y reparaciones | 688.70 | 702.59 | 6929.63 | 8758.08 | 1.02 | 1.26 |
| Agricultura, pesca, ganadería, caza, silvicultura, minas y canteras | 1271.73 | 821.86 | 8234.50 | 10630.56 | 0.65 | 1.29 |
| Servicios sociales y de salud | 749.03 | 766.22 | 6417.54 | 9123.04 | 1.02 | 1.42 |
| Intermediación financiera | 1113.06 | 1684.70 | 9407.47 | 16673.18 | 1.51 | 1.77 |
| Transporte, almacenamiento y comunicación | 807.07 | 1005.46 | 6000.78 | 12598.11 | 1.25 | 2.10 |
| Suministro de electricidad, gas, agua y afines | 1204.92 | 1497.05 | 7284.62 | 19862.60 | 1.24 | 2.73 |
En la Tabla se aprecia que hay distintos niveles de brechas según ramas, mostrando una gran variabilidad según la actividad que se trate. De hecho, hay ramas en donde las brechas son menores a 1, es decir, que los ingresos promedio en EPH son mayores que en MLER. Las brechas parecerían ordenarse según nivel de ingresos, ya que las ramas de mayores ingresos son las que generalmente tienen mayor brecha de ingresos reportados según la fuente de datos.
A continuación analizamos cómo varía la brecha según deciles de ingresos, para analizar las diferencias entre EPH y MLER en los distintos segmentos de la distribución.
#Calculo deciles
EPH2003_2014 <- EPH %>%
filter(ANO4 %in% c(2003, 2014)) %>%
group_by(ANO4) %>%
mutate(decil=ntile(Remun_neta, 10))
MLER2003_2014 <- MLER %>%
filter(ANO4 %in% c(2003, 2014)) %>%
group_by(ANO4) %>%
mutate(decil=ntile(Remun_neta, 10))
comp_deciles <- EPH2003_2014 %>%
select(ANO4, Remun_neta, decil, PONDERA, Base) %>%
bind_rows(MLER2003_2014 %>% ungroup() %>% select(ANO4,Remun_neta, decil, PONDERA, Base))
Resumen_deciles <- comp_deciles %>%
group_by(ANO4, Base, decil) %>%
summarise(Casos = n(),
Casos_expand = sum(PONDERA),
Ingresos_promedio = weighted.mean(Remun_neta,PONDERA))
Resumen_deciles_comparativo <- Resumen_deciles %>%
select(ANO4, Base, decil, Ingresos_promedio) %>%
pivot_wider(names_from = c(ANO4, Base), values_from = Ingresos_promedio,
names_prefix = "ingresos_") %>%
filter(!is.na(decil)) %>%
mutate(brecha2003= ingresos_2003_MLER / ingresos_2003_EPH,
brecha2014= ingresos_2014_MLER / ingresos_2014_EPH) %>%
arrange(brecha2014)
colnames(Resumen_deciles_comparativo) <- c("Decil", "EPH 2003", "MLER 2003", "EPH 2014", "MLER 2014", "Brecha 2003", "Brecha 2014")
Resumen_deciles_comparativo %>%
kable(digits=2, caption = "<b>Tabla 7.</b> Promedios de remuneraciones netas y brecha entre bases según deciles de ingresos. EPH y MLER. Años 2003 y 2014.")%>%
kable_styling("striped")| Decil | EPH 2003 | MLER 2003 | EPH 2014 | MLER 2014 | Brecha 2003 | Brecha 2014 |
|---|---|---|---|---|---|---|
| 1 | 220.28 | 156.75 | 2482.58 | 1640.98 | 0.71 | 0.66 |
| 2 | 347.20 | 287.70 | 4044.49 | 3561.03 | 0.83 | 0.88 |
| 3 | 409.07 | 374.19 | 4918.68 | 4753.71 | 0.91 | 0.97 |
| 4 | 480.91 | 445.70 | 5756.00 | 5948.08 | 0.93 | 1.03 |
| 5 | 529.02 | 514.41 | 6309.90 | 7134.56 | 0.97 | 1.13 |
| 6 | 616.18 | 588.76 | 7084.07 | 8279.00 | 0.96 | 1.17 |
| 7 | 724.94 | 698.06 | 7960.16 | 9363.99 | 0.96 | 1.18 |
| 8 | 851.42 | 871.45 | 9215.58 | 11132.65 | 1.02 | 1.21 |
| 9 | 1146.71 | 1190.97 | 11164.47 | 14222.75 | 1.04 | 1.27 |
| 10 | 2198.84 | 3193.34 | 17141.83 | 29614.76 | 1.45 | 1.73 |
En el cuadro se observan distintas características de la brecha entre los ingresos reportados por la MLER y la EPH. En primer lugar, se observa que para 2003 la brecha entre bases era un problema exclusivo del decil de más altos ingresos, ya que los otros 9 deciles presentaban brechas menores a 1,04, es decir que la brecha era casi inexistente o incluso los ingresos de EPH superan a los de la MLER. En el año 2014, la situación cambia y se registran brechas importantes a partir del quinto decil, las cuales muestran un progresivo aumento a medida que se va subiendo de decil de ingresos. Esto se relaciona con la distribución de ingresos que estudiamos al comienzo de este apartado: como la distribución de la remuneración bruta en la MLER tiene un sesgo hacia los extremos, es de esperar que la brecha de ingresos para los deciles más bajos de ingresos sea menor a 1 en comparación a la EPH y que sea mayor a 1 en los deciles de ingresos más altos.
3. Entrenamiento de modelos para la imputación
En esta sección realizamos el ajuste de regresiones y el entrenamiento de árboles de decisión para la posterior imputación de ingresos. Realizamos el ejercicio de imputación de ingresos para el año 2014, último período con datos en la MLER para el mes utilizado de referencia. En primer lugar, realizamos una limpieza de datos, eliminando casos con valores faltantes y outliers. Luego particionamos los datos de la MLER en datasets de entrenamiento y testeo, ajustamos los modelos de imputación de ingresos, para después comparar la performance de los modelos analizando el RMSE y el MAE para los datasets de entrenamiento y testeo. Por último, elegimos el mejor modelo para la imputación de ingresos en la EPH.
3.1. Preparación de los datasets de entrenamiento e imputación
variables_dataset <- c("ingreso", "edad", "edad2", "sexo", "region", "tramo_tamanio", "percentil", "rama")
EPH2014 <- readRDS("Fuentes/dataset_EPH.RDS") %>%
filter(ANO4==2014) %>%
select(variables_dataset)
MLER2014 <- readRDS("Fuentes/dataset_MLER.RDS") %>%
filter(ANO4==2014) %>%
select(variables_dataset)3.1.1. Limpieza de datos
Antes de entrenar los modelos realizamos una limpieza de las bases preprocesadas en el análisis exploratorio de la sección anterior. En primer lugar quitamos los casos con valores faltantes, dejando en las bases tan sólo a los registros con datos completos en todas sus variables. En segundo lugar, estudiamos las variables ingreso y edad de la MLER para el año 2014, con el objetivo de verificar si hay datos con valores atípicos para quitarlos de la muestra y excluirlos de los datasets de entrenamiento y testeo. Eliminamos los 373 casos con ingresos mayores a la suma de la media de ingresos más 5 desvíos estándar, que representan menos del 0,3% de los casos presentes en la base. Realizamos el mismo procedimiento para la variable edad y no encontramos datos atípicos en edad según el criterio de 5 desvíos estándar.
#Chequeamos cuantos casos con valores faltantes hay por variable
na_count_MLER2014 <-sapply(MLER2014, function(y) sum(length(which(is.na(y)))))
na_count_MLER2014 <- data.frame(na_count_MLER2014)
na_count_EPH2014 <-sapply(EPH2014, function(y) sum(length(which(is.na(y)))))
na_count_EPH2014 <- data.frame(na_count_EPH2014)
#Nos quedamos con casos sin valores faltantes
MLER2014 <- MLER2014[complete.cases(MLER2014), ]
EPH2014 <- EPH2014[complete.cases(EPH2014), ]
#Eliminamos outliers en ingresos y edad segun el criterio de 5 desvios estandar
ingreso_promediomas5desvios = mean(MLER2014$ingreso)+5*sd(MLER2014$ingreso)
resumen_outliers_ingresos <- MLER2014 %>% summarise(arriba=sum(ingreso>ingreso_promediomas5desvios),
abajo=sum(ingreso<=ingreso_promediomas5desvios)) %>%
mutate(por_arriba=arriba/(arriba+abajo))
MLER2014 <- MLER2014 %>% filter(ingreso<ingreso_promediomas5desvios)3.1.2. Separación en datasets de train y test
set.seed(45848748)
indices <- sample(1:nrow(MLER2014), size = 0.80 * nrow(MLER2014))
df_train <- MLER2014[indices,]
df_test <- MLER2014[-indices,]
saveRDS(MLER2014, "Fuentes/dataset_MLER2014_limpio.RDS")
saveRDS(df_train, "Fuentes/dataset_MLER_train.RDS")
saveRDS(df_test, "Fuentes/dataset_MLER_test.RDS")
saveRDS(EPH2014, "Fuentes/dataset_EPH2014.RDS")
comp_casos_train <- df_train %>%
summarise(Casos=n(),
'Casos ponderados'=Casos*33.33333)
comp_casos_test <- df_test %>%
summarise(Casos=n(),
'Casos ponderados'=Casos*33.33333)
col <- c('Dataset de entrenamiento', 'Dataset de testeo') %>% as.data.frame()
comp_casos_df <- bind_rows(comp_casos_train, comp_casos_test)
comp_casos <- bind_cols(col, comp_casos_df)
colnames(comp_casos)[1] <- ""
comp_casos %>%
kable(caption = "<b>Tabla 8. </b> Tamaño de los datasets de train y test") %>%
kable_styling(bootstrap_options = c("striped"), full_width = F, position = "float_right")| Casos | Casos ponderados | |
|---|---|---|
| Dataset de entrenamiento | 133859 | 4461966 |
| Dataset de testeo | 33465 | 1115500 |
Separamos el dataset de la MLER en particiones de train y test, en una proporción de 80-20 respectivamente. El dataset de entrenamiento está compuesto por 133.859 casos, que ponderados representan a 4.461.966 relaciones laborales. Por su parte, la partición de testeo es el 20% del dataset original preprocesado e incluye 33.465 casos, que representan a unos 1.115.500 puestos de trabajo.
3.1.3. Ajuste de regresiones lineales
#Cargamos datasets de entrenamiento y testeo
df <- readRDS("Fuentes/dataset_MLER2014_limpio.RDS")
df_train <- readRDS("Fuentes/dataset_MLER_train.RDS")
df_test <- readRDS("Fuentes/dataset_MLER_test.RDS")En primer lugar estimamos dos modelos lineales simples para utilizarlos de base de comparación con los árboles de decisión. El primer modelo inicial predice el ingreso por medio de una regresión lineal múltiple con las siguientes variables predictoras: edad, sexo y percentil.
#Ajustamos un modelo lineal multiple con 3 variables basicas
modelo1 <- lm(ingreso ~ edad + sexo + percentil, data=df_train)
saveRDS(modelo1, file='./modelos/modelo1.rds', compress=FALSE)
tidy_modelo1 <- tidy(modelo1, conf.int = TRUE)
tidy_modelo1 %>%
kable(caption = "<b>Tabla 9. </b> Resumen Modelo Lineal I", digits=3) %>%
kable_styling(bootstrap_options = c("striped")) | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | -2048.196 | 45.707 | -44.812 | 0 | -2137.781 | -1958.612 |
| edad | 19.263 | 1.044 | 18.452 | 0 | 17.217 | 21.310 |
| sexoHombre | 380.020 | 26.183 | 14.514 | 0 | 328.703 | 431.338 |
| percentil | 203.154 | 0.434 | 468.471 | 0 | 202.304 | 204.004 |
Las tres variables predictoras poseen coeficientes de regresión estadísticamente significativos y sus intervalos de confianza no incluyen al cero. Según este modelo, ser de sexo masculino agrega $380 al ingreso esperado con respecto a ser mujer ceteris paribus las demás variables. De manera análoga, cada año de edad suma $19 al ingreso esperado y cada percentil suma $203.
#Observamos más información sobre el modelo
glance(modelo1)%>%
kable(digits=3) %>%
kable_styling(bootstrap_options = c("striped"))El R al cuadrado del Modelo I nos indica que éste explica menos del 61% de la variabilidad total del problema, por lo que podría concluirse que no es un modelo muy potente para predecir los ingresos de la MLER. Más adelante en este trabajo, nos centraremos en las métricas de performance sobre el dataset de train para evaluar cuál es el desempeño de este modelo base en comparación a las regresiones basadas en árboles de decisión.
El segundo modelo utiliza todas las variables del dataset de entrenamiento como predictoras.
#Ajustamos un modelo lineal multiple con todas las variables disponibles.
modelo2 <- lm(ingreso ~ edad + edad2 + sexo + region + tramo_tamanio + percentil + rama, data=df_train)
saveRDS(modelo2, file='./modelos/modelo2.rds', compress=FALSE)
tidy_modelo2 <- tidy(modelo2, conf.int = TRUE)
tidy_modelo2 %>%
kable(caption = "<b>Tabla 10. </b> Resumen Modelo Lineal II", digits=3) %>%
kable_styling(bootstrap_options = c("striped")) %>%
scroll_box( height = "300px")| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| (Intercept) | -1786.250 | 129.192 | -13.826 | 0.000 | -2039.463 | -1533.036 |
| edad | 13.263 | 6.288 | 2.109 | 0.035 | 0.939 | 25.586 |
| edad2 | 0.078 | 0.074 | 1.053 | 0.293 | -0.067 | 0.222 |
| sexoHombre | 182.620 | 28.121 | 6.494 | 0.000 | 127.503 | 237.737 |
| regionCentro | -367.935 | 27.799 | -13.236 | 0.000 | -422.420 | -313.450 |
| regionNEA | -197.163 | 61.150 | -3.224 | 0.001 | -317.016 | -77.310 |
| regionCuyo | -181.931 | 53.279 | -3.415 | 0.001 | -286.356 | -77.506 |
| regionPatagonia | 1615.928 | 53.976 | 29.938 | 0.000 | 1510.136 | 1721.720 |
| regionNOA | -354.107 | 50.799 | -6.971 | 0.000 | -453.672 | -254.542 |
| tramo_tamaniotramo 2 | -367.244 | 36.007 | -10.199 | 0.000 | -437.817 | -296.672 |
| tramo_tamaniotramo 3 | -254.020 | 38.933 | -6.525 | 0.000 | -330.329 | -177.712 |
| tramo_tamaniotramo 4 | 456.359 | 36.778 | 12.408 | 0.000 | 384.274 | 528.444 |
| percentil | 200.682 | 0.516 | 389.026 | 0.000 | 199.671 | 201.693 |
| ramaAgricultura, pesca, ganadería, caza, silvicultura, minas y canteras | 2500.993 | 58.197 | 42.974 | 0.000 | 2386.928 | 2615.058 |
| ramaComercio y reparaciones | -578.471 | 42.848 | -13.501 | 0.000 | -662.451 | -494.490 |
| ramaConstrucción | 433.727 | 55.652 | 7.794 | 0.000 | 324.650 | 542.803 |
| ramaEnseñanza, actividades profesionales, científicas y técnicas | 495.744 | 59.357 | 8.352 | 0.000 | 379.404 | 612.083 |
| ramaHoteles y restaurantes | 522.050 | 66.298 | 7.874 | 0.000 | 392.108 | 651.992 |
| ramaIndustrias manufactureras | 58.867 | 42.222 | 1.394 | 0.163 | -23.888 | 141.621 |
| ramaIntermediación financiera | 651.473 | 80.702 | 8.073 | 0.000 | 493.298 | 809.648 |
| ramaServicios comunitarios, sociales, personales y de entretenimiento | -135.950 | 57.678 | -2.357 | 0.018 | -248.997 | -22.903 |
| ramaServicios sociales y de salud | -486.409 | 63.248 | -7.691 | 0.000 | -610.374 | -362.444 |
| ramaSuministro de electricidad, gas, agua y afines | 4315.565 | 143.964 | 29.977 | 0.000 | 4033.397 | 4597.733 |
| ramaTransporte, almacenamiento y comunicación | -236.871 | 51.403 | -4.608 | 0.000 | -337.620 | -136.123 |
En este caso la edad al cuadrado no es una variable explicativa con coeficientes estadísticamente significativo, a diferencia del resto de las variables que sí lo tienen. En el caso de la variable categórica rama, si bien el regresor de la categoría ‘Servicios comunitarios, sociales…’ no es significativo, la totalidad de las categorías sí lo son según el test F de significatividad.
#Observamos más información sobre el modelo
glance(modelo2)%>%
kable(digits=3) %>%
kable_styling(bootstrap_options = c("striped"))No obstante, este segundo modelo base tampoco parece mejorar mucho su desempeño con respecto al modelo anterior en lo que respecta a la porción de la variabilidad explicada. El R2 y el R2 ajustado son de 0,658.
3.2. Entrenamiento de modelos de machine learning
Para este trabajo utilizaremos árboles de decisión. Estos son modelos
predictivos basados en algoritmos de entrenamiento supervisado, que
trabajan dividiendo recursivamente los datos de entrenamiento y
estableciendo reglas lógicas que luego pueden aplicarse a nuevos casos
para realizar predicciones. Se trata de un método ampliamente utilizado
para problemas de clasificación, pero que también puede utilizarse para
problemas de regresión, como en el caso del presente trabajo. Entrenamos
los modelos en lenguaje R utilizando el paquete caret para
los algoritmos Random Forest y XGBoost, y el paquete
lightgbm para el entrenamiento del modelo
LightGBM. La aplicación de los dos primeros está basada en el
trabajo de Rosatti (2021), mientras que la implementación de
LightGBM fue realizada en base a la documentación del paquete
mencionado.
Estos algoritmos se basan en técnicas de ensamble de modelos, también llamados ‘clasificadores basados en comités’. La idea de estos ensambles es generar distintos modelos de clasificación o regresión y decidir la predicción final por algún sistema de votación o elección entre los outputs de los distintos modelos (de ahí la metáfora de un ‘comité’ integrado por distintos especialistas). De esta manera se procura reforzar la capacidad predictiva de los árboles de decisión, generando árboles sobre muestras o submuestras del dataset de entrenamiento que tomen en cuenta distintos aspectos del problema a resolver.
Random Forest trabaja realizando submuestras en términos de casos y de atributos. Estima iterativamente árboles utilizando un subconjunto de todos los predictores disponibles y una submuestra del dataset mediante la técnica de bootstrap. Una vez entrenado el modelo, las predicciones se realizan por medio de voto mayoritario entre los árboles generados. XGBoost también se basa en metodologías de remuestreo, pero utiliza la técnica de boosting, la cual no realiza muestreos del dataset aleatoriamente, sino que se focaliza en aquellos casos en los que el modelo funciona peor. Por último, LightGBM es una mejora del modelo anterior que utiliza histogramas de distribución para la división de las variables discretas en bins o compartimientos. Esto hace mucho más eficiente la estimación de árboles y permite tunear hiperparámetros con menor costo computacional. También se diferencia de XGBoost en que la construcción de árboles no se realiza por niveles sino por hojas.
Para el entrenamiento de los modelos utilizamos la partición de entrenamiento de los datos preprocesados de la MLER para el mes de agosto del 2014, sobre la cual realizamos un tunning de hiperparámetros mediante la metodología de 5-fold cross-validation.
## Defino índices para CV
### Index para cross validation
set.seed(4584)
cv_index <- createFolds(y = df_train$ingreso,
k=5,
list=TRUE,
returnTrain=TRUE)
fitControl <- trainControl(
index=cv_index,
method="cv",
number=5,
verbose = TRUE,
allowParallel=TRUE)
### Esquema de validación cruzada para estimación de error
set.seed(8748)
cv_index2 <- createFolds(y = df_train$ingreso,
k=5,
list=TRUE,
returnTrain=TRUE)
fitControl2 <- trainControl(
index=cv_index2,
method="cv",
number=5,
allowParallel=TRUE)3.2.1. Random Forest
### Definición de grid de hiperparámetros
rfGrid <- expand.grid(mtry=c(1:(ncol(df_train)-1)),
splitrule='variance',
min.node.size=c(100, 500, 1000, 2000))
fitControl$allowParallel<-TRUE
### Activar procesamientos paralelos
cl <- makeCluster(6, type='PSOCK', outfile='')
registerDoParallel(cl)
### Reloj
t0<-proc.time()
### Entrenamiento - Tuning hiperparámetros
rf_model<- train(ingreso ~ ., data = df_train,
method = "ranger",
trControl = fitControl,
tuneGrid = rfGrid,
metric='RMSE')
### Outputs
reloj <- proc.time() - t0
print(paste0("Este árbol tardó ", round(reloj[3]/60, 0) , " minutos en correr"))
# Selecting tuning parameters
# Fitting mtry = 7, splitrule = variance, min.node.size = 100 on full training set
# "Este árbol tardó 54 minutos en correr"
### Entrenamiento - Mejor modelo entrenado sobre fitControl2
### Reloj
t0<-proc.time()
### Entrenamiento - Tuning hiperparámetros
rf_model<-train(ingreso ~ ., data = df_train,
method = "ranger",
trControl = fitControl2,
verbose = TRUE,
tuneGrid = rf_model$bestTune,
metric='RMSE')
stopCluster(cl)
### Outputs
reloj <- proc.time() - t0
print(paste0("Este arbol tardo ", round(reloj[3]/60, 0) , " minutos en correr"))
saveRDS(rf_model, file='./modelos/rf.rds', compress=FALSE)
print(rf_model$bestTune)Para el tuneo de hiperparámetros en Random Forest se utilizó
un grid con el parámetro mtry, que setea la
cantidad de atributos que son tomados para la construcción de cada
árbol, y el parámetro min.node.size, que determina la
cantidad de mínima de casos para cada nodo terminal del árbol (de modo
que un menor número en este parámetro resulta en árboles de mayor
profundidad). El mejor modelo encontrado por Random Forest es
el que toma las 7 variables predictoras para la estimación de árboles y
100 observaciones como mínimo número para los nodos terminales. El
algoritmo tomó 56 minutos en correr sobre la base de datos de
entrenamiento utilizada.
3.2.2. XGBoost
### Definición de grid de hiperparámetros
xgbGrid <- expand.grid(nrounds = c(10, 50, 100, 200),
max_depth = c(5, 10, 15, 20),
colsample_bytree = c(0.5, 0.75),
eta = c(0.3, 0.1, 0.05, 0.01),
gamma=c(0, 0.01),
min_child_weight = c(1, 0),
subsample = 1)
fitControl$allowParallel<-TRUE
### Activar procesamientos paralelos
cl <- makeCluster(6, type='PSOCK', outfile='')
registerDoParallel(cl)
### Reloj
t0<-proc.time()
### Entrenamiento - Tuning hiperparámetros
xgb_model<-train(ingreso ~ ., data = df_train,
method = "xgbTree",
trControl = fitControl,
verbose = TRUE,
tuneGrid = xgbGrid,
metric='RMSE')
stopCluster(cl)
### Entrenamiento - Mejor modelo entrenado sobre fitControl2
### Activar procesamientos paralelos
cl <- makeCluster(6, type='PSOCK', outfile='')
registerDoParallel(cl)
fitControl2$allowParallel<-TRUE
### Entrenamiento - Tuning hiperparámetros
xgb_model<-train(ingreso ~ ., data = df_train,
method = "xgbTree",
trControl = fitControl2,
verbose = TRUE,
tuneGrid = xgb_model$bestTune,
metric='RMSE')
stopCluster(cl)
### Outputs
reloj <- proc.time() - t0
print(paste0("Este arbol tardo ", round(reloj[3]/60, 0) , " minutos en correr"))
saveRDS(xgb_model, file='./modelos/xgb.rds', compress=FALSE)
print(xgb_model$bestTune)
# Selecting tuning parameters
# Fitting nrounds = 100, max_depth = 5, eta = 0.1, gamma = 0.01, colsample_bytree = 0.75, min_child_weight = 0, subsample = 1 on full training set
# [1] "Este arbol tardo 207 minutos en correr"Para el modelo XGBoost utilizamos un grid con 6
hiperparámetros: el número de iteraciones, la profundidad máxima, los
parámetros eta y gamma, y la cantidad mínima
de observaciones requeridas para crear un nuevo nodo. El mejor árbol
encontrado utiliza 100 iteraciones, tiene una profundidad máxima de 5,
sus parámetros eta y gamma son iguales a 0,1 y
no utiliza una cantidad mínima de observaciones para crear un nuevo
nodo. Este árbol tomó 207 minutos en correr sobre el dataset de
entrenamiento.
3.2.3. LightGBM
#Creamos un encoding compatible con LightGBM para las variables categóricas
encode_ordinal <- function(x, order = unique(x)) {
x <- as.numeric(factor(x, levels = order, exclude = NULL))
x
}
df_train_lgbm <- df_train %>%
mutate(sexo=encode_ordinal(sexo),
region= encode_ordinal(region),
tramo_tamanio=encode_ordinal(tramo_tamanio),
rama=encode_ordinal(rama)) %>%
data.matrix()
df_test_lgbm <- df_test %>%
mutate(sexo=encode_ordinal(sexo),
region= encode_ordinal(region),
tramo_tamanio=encode_ordinal(tramo_tamanio),
rama=encode_ordinal(rama)) %>%
data.matrix()
#Creo bases en formato lgb
y_train <- df_train_lgbm[,1]
y_test <- df_test_lgbm[,1]
train_lgb <- lgb.Dataset(as.matrix(df_train_lgbm[,2:ncol(df_train_lgbm)]),
categorical_feature = c("sexo", "region", "tramo_tamanio", "rama"),
label=y_train)
test_lgb <- lgb.Dataset.create.valid(train_lgb,
as.matrix(df_test_lgbm[,2:ncol(df_test_lgbm)]),
categorical_feature = c("sexo", "region", "tramo_tamanio", "rama"),
label = y_test)
lgbmGrid <- expand.grid(max_depth = seq(2,16,2),
num_leaves =seq(8,24,4),
num_iterations = seq(200, 1200, 200),
early_stopping_rounds=seq(50,300,50),
learning_rate = seq(.40, .60, .5))
rmse_fit <- NULL
rmse_predict <- NULL
### Reloj
t0<-proc.time()
for (j in 1:nrow(lgbmGrid)) {
set.seed(4584)
light_gbn_tuned <- lgb.train(
params = list(
objective = "regression",
metric = "rmse",
max_depth = lgbmGrid$max_depth[j],
num_leaves =lgbmGrid$num_leaves[j],
num_iterations = lgbmGrid$num_iterations[j],
early_stopping_rounds=lgbmGrid$early_stopping_rounds[j],
learning_rate = lgbmGrid$learning_rate[j]
),
valids = list(test = test_lgb),
data = train_lgb
)
yhat_fit_tuned <- predict(light_gbn_tuned,as.matrix(df_train_lgbm[,2:ncol(df_train_lgbm)]))
yhat_predict_tuned <- predict(light_gbn_tuned,as.matrix(df_test_lgbm[,2:ncol(df_test_lgbm)]))
rmse_fit[j] <- RMSE(y_train,yhat_fit_tuned)
rmse_predict[j] <- RMSE(y_test,yhat_predict_tuned)
cat(j, "\n")
}
min(rmse_fit)
min(rmse_predict)
lgbmGrid[which.min(rmse_fit),]
lgbmGrid[which.min(rmse_predict),]
rmse_diff <- rmse_fit - rmse_predict
rmse_models <- data.frame(rmse_fit,rmse_predict, rmse_diff)
rmse_models_sort <- rmse_models[order(rmse_diff),]
print(lgbmGrid[which.min(rmse_predict),])
reloj <- proc.time() - t0
print(paste0("Este árbol tardó ", round(reloj[3]/60, 0) , " minutos en correr"))
#[1] "Este árbol tardó 49 minutos en correr"
#max_depth num_leaves num_iterations early_stopping_rounds learning_rate
# 4 12 200 50 0.4
#Modelo final con mejores hiperparámetros
lgbm_model <- lgb.train(
params = list(
objective = "regression",
metric = "rmse",
max_depth = 4,
num_leaves =12,
num_iterations = 200,
early_stopping_rounds=50,
learning_rate = 0.4
),
valids = list(test = test_lgb),
data = train_lgb
)
lgb.save(lgbm_model, file='./modelos/lgbm.rds')En el caso de LightGBM, el grid contiene los
hiperparámetros de máxima profundidad, cantidad de hojas, número de
iteraciones, la tasa de aprendizaje y el parámetro
early_stopping_rounds, que le indica al algoritmo a partir
de qué punto se deben detener las iteraciones, una vez que ya no se
tienen ganancias en la función objetivo. El mejor modelo encontrado
tiene una profundidad máxima de 4, 12 hojas como máximo, realiza 200
iteraciones, el parámetro early_stopping_rounds es igual a
50 y su learning_rate es 0.4. Este árbol fue el que menos
tiempo tomó en correr sobre el dataset de entrenamiento, tardando 49
minutos.
3.3. Comparación de peformances de los modelos estimados
# Cargo modelos
modelo1 <- readRDS("~/GitHub/imputacion_MLER_EPH/modelos/modelo1.rds")
modelo2 <- readRDS("~/GitHub/imputacion_MLER_EPH/modelos/modelo2.rds")
rf_model <- readRDS("~/GitHub/imputacion_MLER_EPH/modelos/rf.rds")
xgb_model <- readRDS("~/GitHub/imputacion_MLER_EPH/modelos/xgb.rds")
lgbm_model <- lgb.load("~/GitHub/imputacion_MLER_EPH/modelos/lgbm.rds")
y_train <- df_train %>%
select(c("ingreso"))
y_test <- df_test %>%
select(c("ingreso"))
x_train <- df_train %>%
select(-c("ingreso"))
x_test <- df_test %>%
select(-c("ingreso"))
#Preparo base para lightgbm
encode_ordinal <- function(x, order = unique(x)) {
x <- as.numeric(factor(x, levels = order, exclude = NULL))
x
}
x_test_lgbm <- x_test %>%
mutate(sexo=encode_ordinal(sexo),
region= encode_ordinal(region),
tramo_tamanio=encode_ordinal(tramo_tamanio),
rama=encode_ordinal(rama)) %>%
data.matrix()
x_train_lgbm <- x_train %>%
mutate(sexo=encode_ordinal(sexo),
region= encode_ordinal(region),
tramo_tamanio=encode_ordinal(tramo_tamanio),
rama=encode_ordinal(rama)) %>%
data.matrix()
#Generación de predicciones en train y medicion de performance
lista_modelos <- c("ingreso", "Modelo lineal I", "Modelo lineal II", "Random Forest", "XGBoost", "LightGBM")
pred.modelo1 <- predict(modelo1, x_train)
pred.modelo2 <- predict(modelo2, x_train)
pred.rf <- predict(rf_model , x_train)
pred.xgb <- predict(xgb_model, x_train)
pred.lgbm <- predict(lgbm_model, x_train_lgbm)
predicciones.fit <- data.frame(y_train, pred.modelo1, pred.modelo2, pred.rf, pred.xgb, pred.lgbm )
colnames(predicciones.fit) <- lista_modelos
rmse <- predicciones.fit %>%
summarise( across(2:6, .fns = ~rmse(ingreso, .)))
mae <- predicciones.fit %>%
summarise( across(2:6, .fns = ~mae(ingreso, .)))
medidas.preformance <- c("RMSE", "MAE")
performance_modelos_train <- bind_rows(rmse, mae)
row.names(performance_modelos_train) <- c("RMSE", "MAE")
performance_modelos_train <- t.data.frame(performance_modelos_train)
performance_modelos_train %>%
kable(caption = "<b>Tabla 11. </b> Performance sobre el dataset de entrenamiento",
digits=0) %>%
kable_styling(bootstrap_options = c("striped"), full_width = F, position = "float_left")| RMSE | MAE | |
|---|---|---|
| Modelo lineal I | 4438 | 2284 |
| Modelo lineal II | 4321 | 2240 |
| Random Forest | 1056 | 291 |
| XGBoost | 917 | 213 |
| LightGBM | 937 | 188 |
En este apartado comparamos la performance de los modelos entrenados. Las métricas que utilizamos para evaluar el desempeño de los modelos son la raíz cuadrada del error cuadrático medio (RMSE por sus siglas en inglés) y el error absoluto medio (o MAE). La RMSE mide la raíz del promedio de las diferencias entre valores observados y predichos, dando una mayor importancia a los errores de mayor magnitud. Por su parte, el MAE mide la magnitud promedio de los errores sin importar el signo de la diferencia entre valores observados y predichos. Cuanto mejor sea la performance del modelo, menor son la RMSE y el MAE. Comencemos analizando cuál es la performance sobre los dataset de entrenamiento.
Como es de esperar, los modelos basados en árboles de decisión presentan un desempeño mucho mejor que las regresiones lineales simples que utilizamos de base. El RMSE de las regresiones simples es superior a los $4.300 y la MAE superior a los $2.200, mientras que los ensambles de árboles se ubican entre los $1100 y $900 de RMSE y entre los $300 y $100 de MAE. Entre los ensambles de árboles, el que menor RMSE presenta es el XGBoost con $917, mientras que el LightGBM muestra el menor MAE con $188.
#Generación de predicciones en test y medicion de performance
pred.modelo1 <- predict(modelo1, x_test)
pred.modelo2 <- predict(modelo2, x_test)
pred.rf <- predict(rf_model , x_test)
pred.xgb <- predict(xgb_model, x_test)
pred.lgbm <- predict(lgbm_model, x_test_lgbm)
predicciones.test <- data.frame(y_test, pred.modelo1, pred.modelo2, pred.rf, pred.xgb, pred.lgbm )
colnames(predicciones.test) <- lista_modelos
rmse <- predicciones.test %>%
summarise( across(2:6, .fns = ~rmse(ingreso, .)))
mae <- predicciones.test %>%
summarise( across(2:6, .fns = ~mae(ingreso, .)))
performance_modelos_test <- bind_rows(rmse, mae)
row.names(performance_modelos_test) <- c("RMSE", "MAE")
performance_modelos_test <- t.data.frame(performance_modelos_test)
performance_modelos_test %>%
kable(caption = "<b>Tabla 12. </b> Performance sobre el dataset de testeo",
digits=0) %>%
kable_styling(bootstrap_options = c("striped"), full_width = F, position = "float_right")| RMSE | MAE | |
|---|---|---|
| Modelo lineal I | 4342 | 2264 |
| Modelo lineal II | 4235 | 2223 |
| Random Forest | 1044 | 300 |
| XGBoost | 917 | 213 |
| LightGBM | 911 | 185 |
La Tabla 12 muestra cuál es la performance de cada uno de los modelos sobre el dataset de testeo. Allí LightGBM aparece como el modelo de mejor performance, tanto si se considera el RMSE como la MAE. Por lo tanto, éste será el modelo que tomaremos para realizar las predicciones sobre la EPH. Es además importante destacar que el costo computacional para correr de LightGBM fue menor que para los otros dos modelos, lo que se expresó en el tiempo que demoró en correr el agloritmo en la misma computadora. El modelo LightGBM tardó 49 minutos en correr, lo que significó un poco menos que el tiempo que tomó Random Forest (56 minutos) y mucho menos que los 207 minutos que demoró XGBoost.
También resulta interesante analizar gráficamente las predicciones de cada uno de los modelos y comprarlas con la distribución de los datos de ingresos reportados por la MLER, para evaluar así los posibles sesgos de los ingresos que estos modelos imputen. A continuación realizamos esta tarea con los datos del dataset de entrenamiento.
predicciones_tidy <- predicciones.test
predicciones_tidy <- predicciones_tidy %>%
gather(key="variable", value=value)
ggplot(predicciones_tidy %>% filter(variable %in% c("ingreso", "Modelo lineal I", "Modelo lineal II")) ,
aes(x=value,group = variable, fill = variable))+
geom_density(alpha = 0.5)+
scale_x_continuous(limits = c(-10000,60000))+
labs(title = "Gráfico 5. Distribucion de ingresos predichos. Modelos lineales",
subtitle = "Ingresos predichos entre -10.000 y 60.000 pesos. Casos muestrales")+
theme(title = element_text(size = 12),
legend.position = "bottom",
legend.title = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.title.x = element_blank(),
legend.key.size = unit(x = 2,units = "mm")) 
Los modelos lineales base ajustados no logran captar la forma de la distribución de ingresos, que presenta dos máximos locales y una larga cola hacia los ingresos de mayor monto. Las regresiones lineales tienen un sesgo importante hacia los valores bajos (incluso presentando muchas predicciones con signo negativo) y también hacia los ingresos altos.
ggplot(predicciones_tidy %>% filter(variable %in% c("ingreso", "Random Forest", "XGBoost", "LightGBM")) ,
aes(x=value,group = variable, fill = variable))+
geom_density(alpha = 0.4)+
scale_x_continuous(limits = c(-10000,60000))+
labs(title = "Gráfico 6. Distribución de ingresos predichos en la base de testeo. Modelos lineales",
subtitle = "Ingresos entre -10.000 y 60.000 pesos. Casos muestrales")+
theme(title = element_text(size = 12),
legend.position = "bottom",
legend.title = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.title.x = element_blank(),
legend.key.size = unit(x = 2,units = "mm")) 
En cambio, los modelos basados en ensambles de árboles de decisión logan captar con mayor detalle la forma que toma la distribución de los ingresos en la MLER. Los tres modelos entrenados logran reproducir los picos de la distribución y ubican casos en la cola derecha de más altos ingresos.
3.4. Realización de imputación con el mejor modelo y comparación con datos originales de la EPH
Ahora utilizamos nuestro mejor modelo, el LightGBM, para imputar ingresos a la EPH. El siguiente gráfico de densidad de distribución compara la distribución de la variable original de ingreso con la distribución de los valores imputados.
#Transformamos dataset EPH a formato lgbm
EPH2014 <- readRDS("Fuentes/dataset_EPH2014.RDS")
EPH2014_lgbm <- EPH2014[, 2:8] %>%
mutate(sexo=encode_ordinal(sexo),
region= encode_ordinal(region),
tramo_tamanio=encode_ordinal(tramo_tamanio),
rama=encode_ordinal(rama)) %>%
data.matrix()
EPH2014$lgbm.imputado <- predict(lgbm_model, EPH2014_lgbm)
EPH2014_tidy <- EPH2014 %>%
select(ingreso, lgbm.imputado) %>%
gather(value = valor, key = variable)
mediaEPH <- mean(EPH2014_tidy$valor[EPH2014_tidy$variable=="ingreso"])
mediaEPHimputado <- mean(EPH2014_tidy$valor[EPH2014_tidy$variable=="lgbm.imputado"])
ggplot(EPH2014_tidy,
aes(x=valor,group = variable, fill = variable))+
geom_density(alpha = 0.5)+
scale_x_continuous(limits = c(0,60000))+
labs(title = "Gráfico 7. Distribución de ingresos de la ocupación principal e ingresos imputados con el mejor modelo encontrado. EPH 2014.",
subtitle = "Ingresos entre 0 y 60000 pesos. Casos muestrales")+
theme(title = element_text(size = 12),
legend.position = "bottom",
legend.title = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.title.x = element_blank(),
legend.key.size = unit(x = 2,units = "mm")) +
geom_vline(xintercept = mean(EPH2014_tidy$valor[EPH2014_tidy$variable=="ingreso"]),
linetype='dashed', color='red') +
geom_vline(xintercept = mean(EPH2014_tidy$valor[EPH2014_tidy$variable=="lgbm.imputado"]), linetype='dashed', color='blue') +
annotate("text",x= mediaEPH *1.85, y=0.00015,
label= paste0("Media EPH: $", round(mediaEPH)), color='red') +
annotate("text",x= mediaEPHimputado *2, y=0.00012,
label= paste0("Media EPH imputado: $", round(mediaEPHimputado)), color='blue')
En el gráfico se observa que los valores imputados concentran menos casos en los dos picos que presenta la distribución de ingresos y ubica más casos en los extremos, tanto a valores altos de la distribución como a valores bajos. También allí se puede ver que la media de los datos imputados es de $9830, mayor a la de los datos originales de la EPH, la cual se ubica en $8025. Otra característica de los valores imputados es que la distribución de casos hacia la cola derecha disminuye de forma más gradual en los datos imputados que en los valores declarados originalmente en la EPH. Esto puede asociarse a que la declaración de ingresos tiende a tener sesgos por el redondeo que realizan los encuestados al responder, los cuales podrían corregirse también por medio de esta imputación, que determina un descenso más gradual de los casos a medida que aumenta el nivel de ingresos.
Mediante esta metodología podríamos entonces reemplazar los datos originales de la EPH con información proveniente de los registros administrativos del SIPA. Esta información, que se encuentra restringida al universo de la fuerza de trabajo asalariada y registrada, podría utilizarse de manera complementaria con los microdatos de la EPH sobre el resto de las categorías ocupacionales y sectores que quedaron fuera del estudio, o bien podría ser utilizada para estudiar exclusivamente las características de la población asalariada registrada. Claro está que la utilización conjunta de datos imputados con datos originales, requeriría a su vez la corrección o evaluación de los datos originales, ya que los problemas de subdeclaración de ingresos también afectan a cuentapropistas, patrones y a quienes reciben otros tipos de transferencias monetarias.
4. Conclusiones y reflexiones finales
En este trabajo realizamos un ejercicio exploratorio de imputación de datos sobre ingresos laborales en la EPH para la población asalariada ocupada en el sector privado, excluyendo al trabajo doméstico remunerado. La imputación se realizó entrenando modelos de machine learning sobre un dataset de entrenamiento con origen en la MLER-SIPA. Esta última es una muestra de registros declarados por los empleadores en la seguridad social, que, por su carácter de declaración jurada, puede considerarse como una fuente de datos más fehaciente que los ingresos autodeclarados por los encuestados de la EPH. Nuestra propuesta de imputación tuvo como objetivo brindar de mayor robustez a los datos sobre ingresos de la EPH y aportar ideas para una posible solución al problema de subdeclaración de ingresos en dicha encuesta.
Los modelos entrenados fueron los ensambles de árboles decisión Random Forest, XGBoost y LightGBM. Realizamos el tuneo de hiperparámetros aplicando 5-fold cross-validation sobre un dataset con datos de la MLER para el mes de agosto del 2014, que previamente fue particionado en datasets de entrenamiento y de testeo. Las medidas de error elegidas para el entrenamiento y para la evaluación de performance fueron el RMSE y la MAE. A partir de la comparación del desempeño de los modelos en el dataset de test, se determinó que LightGBM brindaba el ensamble de árboles con mejor capacidad predictiva. Utilizando dicho modelo procedimos a imputar los datos en la EPH. Esto dio como resultado una distribución de ingresos imputados con más casos hacia las colas de la distribución, principalmente hacia los ingresos más altos. También la distribución de ingresos imputados mostró una distribución más gradual hacia la derecha que en los datos reportados por la EPH.
Como lo adelantamos al comienzo del presente trabajo, nuestro ejercicio de imputación no puede considerarse una solución completa a esta problemática por varios motivos. En primer lugar, la imputación sólo abarca a un subconjunto muy pequeño de la población ocupada, por lo que es necesario completar la imputación asignando o corrigiendo los ingresos laborales del resto de la población ocupada mediante distintas posibles soluciones. Ello podría hacerse complementando los ingresos imputados a la población asalariada registrada con información proveniente de estudios complementarios sobre ingresos de asalariados no registrados, cuentapropistas, trabajadores del sector público y el trabajo doméstico remunerado.
En segundo lugar, nuestro ejercicio tiene el límite de haber sido realizado con la restringida información disponible en la MLER sobre el puesto de trabajo y los atributos de los/as empleados/as. Una imputación de este tipo podría realizarse con un mayor grado de precisión si se contara con información más detallada, que se encuentra disponible en las bases de datos del sistema previsional argentino. En particular, una información que brindaría mayor robustez a nuestros modelos sería la cantidad de horas trabajadas de los/as empleados/as. Esa variable podría incluirse como predictora en los modelos entrenados, o bien se podría realizar el entrenamiento y la imputación en términos de ingresos horarios, en lugar de ingresos totales.
Tercero, es también necesario remarcar que nuestro ejercicio se realizó tan sólo con los datos de agosto de 2014, dado que es el último año con datos disponibles en las dos fuentes. La implementación completa de una metodología de imputación de este tipo debería abarcar la mayor cantidad de casos posible, analizando información para todo el período con disponibilidad de datos. En este caso, el período de recolección del dato pasaría también a ser parte de las variables predictoras y además se haría necesario evaluar la aplicación de algún tipo de corrección por inflación.
A pesar de estos límites, el ejercicio de imputación que realizamos ha mostrado la utilidad de las técnicas de machine learning para la construcción y mejora de las estadísticas socio-laborales, en este caso complementando la información auto-declarada por los encuestados con datos provenientes de los registros administrativos de la seguridad social. De esta manera, el presente trabajo intenta ser un aporte para una necesaria incorporación de este tipo de técnicas en los institutos públicos de estadística.
Bibliografía
ANSES (2011). “Marco conceptual del sistema de estadísticas e indicadores del Sistema Integrado Previsional Argentino”. Observatorio de la Seguridad Social, 2º edición, septiembre 2011.
Diego Kozlowski, Pablo Tiscornia, Guido Weksler, German Rosati and Natsumi Shokida (2020). eph: Argentina’s Permanent Household Survey Data and Manipulation Utilities. R package version https://doi.org/10.5281/zenodo.3462677
Gómez Sabaíni, J. y Rossignolo, D. (2014). “La tributación sobre las altas rentas en América Latina”. Serie Estudios y Perspectivas Nº13, Oficina de la CEPAL en Montevideo. Impreso en Naciones Unidas, Santiago de Chile.
Llach, J. J., & Montoya, S. (1999). En pos de la equidad: La pobreza y la distribución del ingreso en el Area Metropolitana de Buenos Aires: diagnóstico y alternativas políticas. IERA.
Roca, E. y Pena, H. (2001). “La declaración de ingresos en las encuestas de hogares”. Congreso Nacional de Estudios del Trabajo, ASET. Buenos Aires, 1, 2 y 3 de agosto de 2001.
Rosati, G. F. (2021). Métodos de Machine Learning como alternativa para la imputación de datos perdidos. Estudios del Trabajo. Revista de la Asociación Argentina de Especialistas en Estudios del Trabajo (ASET), 61.
Sánchez, M., Pacífico, L., & Kennedy, D. (2016). La participación asalariada en el ingreso y su composición según el vínculo laboral: fuentes de información, metodologías y alternativas de estimación (No. 21). Documentos de Trabajo.