Introducción

Con la intención de generar un análisis reproducible, fácil y entendible, se elabora el siguiente código para poder generar un primer análisis del Cuestionario Sociodemográfico de la ENOE. El código abarca los siguientes pasos:

Sin antes, se ha de mencionar que, al ser un análisis sencillo no se profundiza teóricamente el diseño de muestra y en los supuestos de los modelos, por lo que es tarea de la persona lectora de este código buscar la información complementaría.Asimismo, este código es la segunda parte del dos publicaciones, en las cuales se análisa el tema en cuestión.

Se agradece cualquier tipo de comentario y se advierte que el código presente se mantiene en actualización.

Putos clave para la discusión de un modelo

  • Objetivo del Modelo
  • Estructura del Modelo
  • Supuesto del Modelo
  • Ajuste del Modelo
  • Selección del Modelo

Librerías de trabajo

A continuación, se enlistas las librerías a utilizar, mismas que serán mencionadas a lo largo del código para poder identificar el uso que tienen.

El siguiente código instala las paqueterías necesarias en caso de no contar previamente con ellas y posteriormente las carga.

if(!require('pacman')) install.packages('pacman')
pacman::p_load(tidyverse,#manipulación de datos
               janitor,#limpieza y tablas de frecuencia
               foreign,#abrir archivos de distintas extenciones
               survey, srvyr, # diseño de encuestas
               scales,#escalas numericas
               stargazer, #Tablas para regresiones
               GGally,#graficos simultaneos
               car,#pruebas
               MASS#Librería selección variables
               )#modelos
#Notación cientifica
options(scipen = 999)

Criterios

Objetivo del Modelo

Estimar los parametros de regresión para la determinación del ingreso en función de las variables:

  • edad
  • sexo
  • escolaridad
  • condición laboral

Función de regresión

\(y = \beta_0 + \beta_1 x + ... + \beta_n+\epsilon\)

Supuestos del modelo de regresión lineal

Linealidad:

  • La relación entre las variables independientes y la variable dependiente es lineal.

Independencia de los errores

  • Las observaciones son independientes entre si

Homocedasticidad

  • La varianza de los errores es constante a lo largo de todas las observaciones

Normalidad de los errores

  • Los errores del modelo estan normalmente distribuidos

No colinealidad y Errores de media cero

  • Ninguna variable independiente es combinación lineal de otras

  • Los errores deben seguir una distribución normal

Pruebas de los supuestos

Descarga

Con el siguiente código descargamos los datos directamente de INEGI y creamos un directorio temporal de trabajo para no cargar la descarga en el actual y liberar memoria, asimismo, se descarga en la dirección de carpetas previamente creada (data/raw) y descomprime, posteriormente se carga el Cuestionario Sociodemográfico SDEMT en nuestro ambiente de trabajo y se elimina el directorio temporal y se librera de la memoria los objetos creados. Para este ejercicio, utilizaremos los datos correspondientes al primer trimestre del 2024.

# url del archivo
url <- "https://www.inegi.org.mx/contenidos/programas/enoe/15ymas/microdatos/enoe_2024_trim1_dbf.zip"
# creación de directorio temporal
td <- tempdir()
# descarga del archivo temporal
tf <-  tempfile(tmpdir=td, fileext=".zip")
download.file(url, tf)
# unzip y crear ruta de carpetas
unzip(tf, exdir="data/raw", overwrite=TRUE)
# cargo en ambiente
sdemt124_raw <-read.dbf('data/raw/ENOE_SDEMT124.dbf')
# elimino td
unlink(td)
unlink(tf)
rm(td,tf,url)

Limpieza

Titulos a formato snake_case, minusculas y espacios con guiones bajos “_”

# Formateamos los titulos
sdemt124 <- sdemt124_raw %>% 
  janitor::clean_names()

Cambiamos el formato de las variables para facilitar la interpretación de los resultados del modelo.

#Variables Numericas
#Edad
sdemt124$edad <- as.numeric(as.character(sdemt124$eda))
#Variables factor
# Sexo
sdemt124$sex <- factor(sdemt124$sex, labels = c("Hombre", "Mujer"), levels = c(1,2))
# Condición de Empleo formar e informal
sdemt124$cond <- factor(sdemt124$emp_ppal, labels = c("informal", "formal"), levels = c(1,2))
#Persona Migrante
sdemt124$mig <- factor(ifelse(as.numeric(sdemt124$l_nac_c) < 33,0,1))

Cabe recalcar que los datos para los cuales se clasifican las variables y no se encuentran coincidencias serán manejados como NA´s.

Diseño de Encuesta

Debido a que el diseño de la ENOE corresponde a una muestra compleja, se deben aplicar un diseño muestral a la base de datos, para ello usamos el paquete srvyr.

# Diseño Muestral
diseño <- sdemt124 %>%
  as_survey_design(
    ids = upm, # Unidades Primarias de Muestreo
    strata = est_d_tri, # Estrato del diseño
    weights = fac_tri, # Factor de expansión
    nest = TRUE) # Diseño anidado

# Confirmación INEGI: 59,120,905    población ocupada
diseño %>% 
  filter(r_def == "00",
        (c_res == "1" | c_res == "3"),
        edad >= 15 & edad <=98)%>% 
  filter(clase2 == 1) %>%
  summarise(Total = survey_total())
## # A tibble: 1 × 2
##      Total Total_se
##      <dbl>    <dbl>
## 1 59120905  274979.

Nos quedamos con las variables que vamos a utilizar para realizar la regresión.

#Con diseño
sdemt_srvy <- diseño %>% 
  filter(r_def == "00",
        (c_res == "1" | c_res == "3"),
        edad >= 15 & edad <=98)%>%
  filter(clase2 == 1,
         ingocup > 0,
         anios_esc < 99)

#Data base (diseño y filtros)
datos <- diseño %>% 
  filter(r_def == "00",
        (c_res == "1" | c_res == "3"),
        edad >= 15 & edad <=98)%>% 
  filter(clase2 == 1,
         ingocup > 0,
         anios_esc < 99) %>% 
  dplyr::select(ingocup, anios_esc, edad, sex, cond, mig)

#Data base de Pob Migrante
dta <- diseño %>% 
  filter(r_def == "00",
        (c_res == "1" | c_res == "3"),
        edad >= 15 & edad <=98)%>% 
  filter(clase2 == 1) %>% 
  filter(mig == 1,
         ingocup > 0,
         anios_esc < 99) %>%
  dplyr::select(ingocup, anios_esc, edad, sex, cond, mig)

#Recortada Pob Mig sin diseño
dt <- sdemt124 %>%
  filter(r_def == "00", 
        (c_res == "1" | c_res == "3"),
        edad >= 15 & edad <=98)%>% 
  filter(clase2 == 1) %>%
  filter(mig == 1,
         ingocup > 0,
         anios_esc < 99) %>%
  dplyr::select(ingocup, anios_esc, edad, sex, cond, mig)

Valores perdidos de la base de datos

#Valores perdidos de mi base de datos
#Base original
sdemt124 %>%
  dplyr::select(ingocup, anios_esc, edad, sex, cond, mig) %>%
  summarise(across(everything(), ~ sum(is.na(.)))) %>% 
  pivot_longer(everything(), names_to = "variable", values_to = "na_count") %>% 
  adorn_totals() %>% 
  as.data.frame()
##    variable na_count
## 1   ingocup        0
## 2 anios_esc        0
## 3      edad     7651
## 4       sex     7651
## 5      cond   229605
## 6       mig     7651
## 7     Total   252558
#Recortada sin diseño y con filtros de la ENOE
dt %>%
  summarise(across(everything(), ~ sum(is.na(.)))) %>%
  pivot_longer(everything(), names_to = "variable", values_to = "na_count") %>%
  adorn_totals() %>%
  as.data.frame()
##    variable na_count
## 1   ingocup        0
## 2 anios_esc        0
## 3      edad        0
## 4       sex        0
## 5      cond        0
## 6       mig        0
## 7     Total        0
str(dt)
## 'data.frame':    754 obs. of  6 variables:
##  $ ingocup  : int  8000 4000 60000 50000 9030 3440 9000 7740 6450 8600 ...
##  $ anios_esc: int  0 18 16 18 9 9 6 8 6 6 ...
##  $ edad     : num  48 73 40 35 25 46 47 57 19 21 ...
##  $ sex      : Factor w/ 2 levels "Hombre","Mujer": 1 2 1 1 1 2 1 1 2 2 ...
##  $ cond     : Factor w/ 2 levels "informal","formal": 2 1 2 2 1 1 1 1 1 1 ...
##  $ mig      : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...
##  - attr(*, "data_types")= chr [1:114] "C" "C" "C" "C" ...

Al aplicar los filtros de la ENOE, se eliminan los NA de las variables y obtenemos variables sin datos perdidos.

Gráfica de distribución de las variables

# Crear un gráfico de pares (pairs plot) usando la función ggpairs de GGally
ggpairs(dt,#datos
        columns = names(dt),# Selecciona todas las columnas de 'dt'
        lower = list(continuous = "smooth"), # Para variables continuas inferiores, usar un ajuste suave (smooth)
        upper = list(continuous = wrap("cor", size = 3)), # En la parte superior del gráfico, para las variables continuas, mostrar la correlación (cor)
        axisLabels = "show")                 # Mostrar etiquetas en los ejes

#Estructura
skimr::skim(dta$variables)
Data summary
Name dta$variables
Number of rows 754
Number of columns 6
_______________________
Column type frequency:
factor 3
numeric 3
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
sex 0 1 FALSE 2 Hom: 460, Muj: 294
cond 0 1 FALSE 2 inf: 476, for: 278
mig 0 1 FALSE 1 1: 754, 0: 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
ingocup 0 1 12499.17 14478.73 167 5160 8600 14000 166667 ▇▁▁▁▁
anios_esc 0 1 11.18 4.94 0 8 12 16 22 ▂▃▇▆▁
edad 0 1 36.49 14.00 15 25 34 47 83 ▇▇▅▂▁

Regresión Lineal Simple

De acuerdo a la escolaridad, se espera que las personas ocupadas tengan un ingreso proporcional a los años de educación con los que cuentan.

#Regresión Lineal Simple de la Población
reg_lin <- svyglm(log(ingocup) ~ anios_esc, design = dta)

#Regresión Lineal Simple de la Muestra
reg_lin.2 <- lm(log(ingocup) ~ anios_esc, data = dt)

#Tabla con summary
stargazer(reg_lin,reg_lin.2, type="text", df=FALSE)
## 
## ================================================
##                         Dependent variable:     
##                     ----------------------------
##                             log(ingocup)        
##                     survey-weighted      OLS    
##                          normal                 
##                           (1)            (2)    
## ------------------------------------------------
## anios_esc               0.074***      0.078***  
##                         (0.009)        (0.006)  
##                                                 
## Constant                8.212***      8.172***  
##                         (0.118)        (0.072)  
##                                                 
## ------------------------------------------------
## Observations              754            754    
## R2                                      0.187   
## Adjusted R2                             0.186   
## Log Likelihood         -1,070.363               
## Akaike Inf. Crit.      2,144.727                
## Residual Std. Error                     0.803   
## F Statistic                          172.713*** 
## ================================================
## Note:                *p<0.1; **p<0.05; ***p<0.01

Resultados de la regresión simple

  • Por cada año adicional de escolaridad, el ingreso mensual predicho aumenta en un 7.8%

Como se observa, los resultados de estimar el modelo sin el diseño muestral se obtienen valores sobreestimados, mientras que, al aplicar el diseño de muestra, a parte de que podemos inferir valroes poblacionales, tenemos valores más ajustados y menores a la regresión sin diseño. Apartir de aquí interpretaremos los resultados para la regresión con el diseño de muestra aplicado.

Intervalo de confianza

confint(reg_lin)
##                  2.5 %     97.5 %
## (Intercept) 7.98069729 8.44326809
## anios_esc   0.05534684 0.09240363

De acuerdo a los intervalos de confianza el cambio en el ingreso se estima entre 5.53% y 9.24% por cada año de educación.

Regresión Lineal Multiple

Apartir de las caracteristicas individuales, determinaremos los estimadores de las variables de estudio para el cambio en el ingreso.

# Regresión Pob Mig y no Mig
reg_mul <- svyglm(log(ingocup) ~ anios_esc + edad + sex + cond,
                 design = datos)

#Pob Migrante
reg_mul.2 <- svyglm(log(ingocup) ~ anios_esc + edad + sex + cond,
                 design = dta)

stargazer(reg_mul,reg_mul.2, type="text", df=FALSE,
          column.labels = c("Población Ocupada", "Población Migrante Ocupada"))
## 
## ==============================================================
##                               Dependent variable:             
##                   --------------------------------------------
##                                   log(ingocup)                
##                   Población Ocupada Población Migrante Ocupada
##                          (1)                   (2)            
## --------------------------------------------------------------
## anios_esc             0.052***               0.056***         
##                        (0.001)               (0.008)          
##                                                               
## edad                  0.002***               0.009***         
##                       (0.0003)               (0.003)          
##                                                               
## sexMujer              -0.331***             -0.448***         
##                        (0.007)               (0.074)          
##                                                               
## condformal            0.525***               0.520***         
##                        (0.008)               (0.069)          
##                                                               
## Constant              8.183***               8.085***         
##                        (0.019)               (0.152)          
##                                                               
## --------------------------------------------------------------
## Observations           131,832                 754            
## Log Likelihood      -165,978.100             -979.284         
## Akaike Inf. Crit.    331,966.100            1,968.569         
## ==============================================================
## Note:                              *p<0.1; **p<0.05; ***p<0.01

Intervalos de confianza para los coeficientes de las variables

confint(reg_mul.2)
##                   2.5 %      97.5 %
## (Intercept)  7.78526349  8.38472374
## anios_esc    0.03927032  0.07192403
## edad         0.00336338  0.01465685
## sexMujer    -0.59246256 -0.30320328
## condformal   0.38444477  0.65558765

Resultados

  • Cada año adicional de educación está asociado con un incremento del 5.6% en el ingreso
  • Cada año adicional de edad está asociado con un incremento del 0.09% en el ingreso laboral
  • Ser mujer está asociado con una disminución del 44.8% (~5 mil pesos) en el ingreso ocupacional en comparación a los hombres
  • Tener un empleo formal está asociado con un incremento del 52% en el ingreso ocupacional
#Analisis de significancia estadística
anova(reg_mul.2, design = dta)
## Anova table:  (Rao-Scott LRT)
## svyglm(formula = log(ingocup) ~ anios_esc, design = dta)
##            stats     DEff       df ddf                   p    
## anios_esc 87.698  1.42672  1.00000 372 0.00000000000005361 ***
## edad      19.174  1.51602  1.00000 371            0.000455 ***
## sex       34.336  1.02713  1.00000 370 0.00000001732644193 ***
## cond      39.405  0.69263  1.00000 369 0.00000000000040403 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Las variables son significativas para explicar el ingreso.

# Con MASS- Librería selección variables
stepAIC(reg_mul.2)
## Start:  AIC=1558.81
## log(ingocup) ~ anios_esc + edad + sex + cond
## 
##             Df Deviance    AIC
## <none>           340.01 1558.8
## - edad       1   351.24 1581.3
## - sex        1   375.56 1631.8
## - cond       1   379.42 1639.5
## - anios_esc  1   382.09 1644.8
## Stratified 1 - level Cluster Sampling design (with replacement)
## With (577) clusters.
## Called via srvyr
## Sampling variables:
##  - ids: upm
##  - strata: est_d_tri
##  - weights: fac_tri
## 
## Call:  svyglm(formula = log(ingocup) ~ anios_esc + edad + sex + cond, 
##     design = dta)
## 
## Coefficients:
## (Intercept)    anios_esc         edad     sexMujer   condformal  
##     8.08499      0.05560      0.00901     -0.44783      0.52002  
## 
## Degrees of Freedom: 753 Total (i.e. Null);  369 Residual
## Null Deviance:       520.6 
## Residual Deviance: 340   AIC: 1969

El criterio de Akaike no hay multicolinealidad.

Diagnóstico de Residuos

#Prueba analítica de normalidad
shapiro.test(reg_mul.2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  reg_mul.2$residuals
## W = 0.96546, p-value = 0.000000000002384

El test de Shapiro concluye que rechazamos la hipótesis nula y concluimos que los datos no están distribuidos normalmente.

library(lmtest)
bptest(reg_mul.2) # prueba analítica de homocedasticidad 
## 
##  studentized Breusch-Pagan test
## 
## data:  reg_mul.2
## BP = 5.3516, df = 4, p-value = 0.2531

\[H_{0}: Homocedasticidad\] \[H_{1}: Heterocedasticiad\]

# Residuos del modelo
residuals <- residuals(reg_mul.2, type = "response")

# Gráficos de residuos
plot(residuals)

hist(residuals)

La distribución de los residuos se mantiene dentro de la media, sin embargo, hay minimos valores en lo largo de las colas de la curva.

ggplot(reg_mul.2, aes(x = reg_mul.2$fitted.values, y = reg_mul.2$residuals)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", col = "red") +
  labs(title = "Gráfico de Residuos vs Valores Ajustados",
       x = "Valores Ajustados",
       y = "Residuos") +
  theme_minimal()

#Tres desviaciones estandar
plot(reg_mul.2, design = dta, which = 3)

  • El patrón en forma de embudo indica heterocedasticidad, donde la varianza de los errores cambia con el nivel de los valores ajustados.
  • De forma gráfica se observan valores sobre la media, sin embargo, existen valores dispersos.
library(car)
crPlots(reg_mul.2, design = dta)

Los residuos por cada regresor vs el ingreso, permite observar la heterogeneidad.

Colinealidad

Una manera de detectar colinealidad es realizar un análisis de componentes principales de las variables independientes. Esta técnica es matemáticamente compleja y aquí se hace sólo un resumen de la misma, es necesario que para entender esta técnica, debe profundizarse en ella, sin embargo se realiza aquí como mero ejercicio.

vif(reg_mul.2)
## anios_esc      edad       sex      cond 
##  1.038090  1.009184  1.113556  1.086077

Considerando los resultados anteriores, se observan valores VIF por debajo de 4, así que las variables independientes ya no presentan multicolinealidad.

####Distribución normal de los residuos:

# Crear un data frame con los residuos y los cuantiles teóricos
qq_data <- data.frame(
  sample = reg_mul.2$residuals,
  theoretical = qqnorm(reg_mul.2$residuals, plot.it = FALSE)$x
)

# Crear el gráfico Q-Q con ggplot2
ggplot(qq_data, aes(sample = sample)) +
  stat_qq() +
  stat_qq_line() +
  labs(title = "Q-Q Plot of Residuals",
       x = "Theoretical Quantiles",
       y = "Sample Quantiles") +
  theme_minimal()

#qqnorm(reg_mul.2$residuals)
#qqline(reg_mul.2$residuals)

La condición de normalidad no se satisface, posiblemente debido a un datos atípicos. Se repite el análisis excluyendo la observación a la que pertenece el residuo atípico.

library(lmtest)
bptest(reg_mul)
## 
##  studentized Breusch-Pagan test
## 
## data:  reg_mul
## BP = 4241.4, df = 4, p-value < 0.00000000000000022

Con un valor menor a 0.05, rechazamos la hipotesis nula y asumimos la existencia de heterocedasticidad.

#Paquete Car: Durbin-Watson (DW)
dwt(reg_mul.2,alternative = "two.sided")
##  lag Autocorrelation D-W Statistic p-value
##    1     -0.02145209      2.042901    0.55
##  Alternative hypothesis: rho != 0

En este caso, el valor es 2.042901, lo que sugiere que no hay autocorrelación significativa.

#Libreria car
outlierTest(reg_mul.2)
##      rstudent  unadjusted p-value     Bonferroni p
## 230 -7.572553 0.00000000000029818 0.00000000022483
## 529 -5.181196 0.00000036410000000 0.00027453000000
## 548  4.465571 0.00001063700000000 0.00802000000000

Para los valores más extremos se obtiene que el valor p ajustado usando la corrección de Bonferroni para los tres casos, es menor a 0.05 por lo que las observaciones son identificadas como outliers significativos.