Determinantes del ingreso laboral en el Cauca. Ecuación de Mincer, GEIH 2025-IV
José Luis
2026-11-03
Parte I. Preparación
1. Paquetes
library(haven) # lectura de archivos .dta
library(dplyr) # manejo de datos
library(ggplot2) # gráficos
library(scales) # formato de ejes
library(sandwich) # errores estándar robustos y agrupados
library(lmtest) # coeftest, waldtest, pruebas de diagnóstico2. Lectura de la base
ruta <- "C:/Users/josec/OneDrive/Desktop/Documentos_Investigacion/Documentos_Mincer"
cauca <- read_dta(file.path(ruta, "GEIH_2025_T4_Cauca.dta"))
dim(cauca)## [1] 6847 5393. Funciones auxiliares
# ---- Formato para texto y gráficos ----
num <- function(x, d = 1) format(round(x, d), big.mark = ".", decimal.mark = ",",
nsmall = d, trim = TRUE)
pesos <- label_number(prefix = "$", big.mark = ".", decimal.mark = ",")
# ---- Formato para tablas de LaTeX: coma decimal {,} y signo menos ----
fx <- function(x, d = 3) {
s <- formatC(x, format = "f", digits = d, big.mark = ".", decimal.mark = ",")
s <- gsub(",", "{,}", s, fixed = TRUE)
sub("^-", "$-$", s)
}
fp <- function(p) ifelse(p < 0.001, "$<${}0{,}001", fx(p, 3))
fi <- function(n) format(n, big.mark = ".", trim = TRUE)
estrellas <- function(p) ifelse(p < 0.01, "***", ifelse(p < 0.05, "**", ifelse(p < 0.1, "*", "")))
fila <- function(etq, val) paste0(etq, " & ", paste(val, collapse = " & "), " \\\\")
grupo <- function(titulo, ncol) paste0("\\multicolumn{", ncol, "}{l}{\\textit{", titulo, "}} \\\\")
# ---- De log-puntos a porcentaje exacto ----
pct <- function(b) 100 * (exp(b) - 1)
# ---- Cuantil ponderado (descriptivos) ----
wq <- function(x, w, p) {
ok <- is.finite(x) & is.finite(w); x <- x[ok]; w <- w[ok]
if (length(x) == 0) return(NA_real_)
o <- order(x); x <- x[o]; w <- w[o]
x[which(cumsum(w) / sum(w) >= p)[1]]
}
# ---- Combinación lineal de coeficientes con su error (estimación, error, p-valor) ----
comb <- function(m, V, pesos_c) {
g <- setNames(rep(0, length(coef(m))), names(coef(m)))
g[names(pesos_c)] <- pesos_c
est <- sum(g * coef(m))
se <- sqrt(as.numeric(t(g) %*% V %*% g))
c(estimacion = est, error = se, p_valor = 2 * pt(-abs(est / se), df.residual(m)))
}
# ---- Cajas ponderadas para violín + caja (bigotes 1,5 RIC acotados a p1-p99) ----
cajas_pond <- function(d, var, ..., log = FALSE) {
d %>%
group_by(...) %>%
summarise(n = n(),
q1 = wq({{ var }}, w, 0.25),
mediana = wq({{ var }}, w, 0.50),
q3 = wq({{ var }}, w, 0.75),
min_d = wq({{ var }}, w, 0.01),
max_d = wq({{ var }}, w, 0.99),
media = if (log) exp(weighted.mean(log({{ var }}), w)) else weighted.mean({{ var }}, w),
.groups = "drop") %>%
mutate(ric = if (log) log(q3) - log(q1) else q3 - q1,
ymin = if (log) pmax(min_d, exp(log(q1) - 1.5 * ric)) else pmax(min_d, q1 - 1.5 * ric),
ymax = if (log) pmin(max_d, exp(log(q3) + 1.5 * ric)) else pmin(max_d, q3 + 1.5 * ric))
}
violin_caja <- function(d, cajas, x, y, colores) {
list(
geom_violin(data = d, aes(x = {{ x }}, y = {{ y }}, fill = {{ x }}, weight = w),
color = NA, alpha = 0.45, scale = "width", width = 0.8, adjust = 1.2),
geom_errorbar(data = cajas, inherit.aes = FALSE,
aes(x = {{ x }}, ymin = ymin, ymax = ymax), width = 0.08, linewidth = 0.6),
geom_boxplot(data = cajas, inherit.aes = FALSE, stat = "identity",
aes(x = {{ x }}, ymin = ymin, lower = q1, middle = mediana, upper = q3, ymax = ymax),
width = 0.16, fill = "white", color = "black", linewidth = 0.6),
geom_point(data = cajas, inherit.aes = FALSE, aes(x = {{ x }}, y = media),
shape = 23, size = 4.5, fill = "gold", color = "black"),
geom_text(data = cajas, inherit.aes = FALSE, aes(x = {{ x }}, y = media, label = etq_media),
nudge_x = 0.2, hjust = 0, size = 3.4, fontface = "bold"),
scale_fill_manual(values = colores)
)
}4. Parámetros
smmlv <- 1423500 # salario mínimo mensual 2025
sem_mes <- 4.285 # semanas por mes (el mismo factor de la variable horatmes del DANE)
jornada_sem <- 44 # jornada legal semanal desde julio de 2025 (Ley 2101 de 2021)
jornada_mes <- jornada_sem * sem_mes # 188,5 horas al mes
smmlv_hora <- smmlv / jornada_mes # salario mínimo por hora
colores_sexo <- c(Hombres = "steelblue", Mujeres = "firebrick")
colores_zona <- c(Urbano = "#7570b3", Rural = "#1b9e77")
colores_emp <- c("Cotiza a pensión" = "#2c7fb8", "No cotiza a pensión" = "#d95f0e")
fuente <- "Fuente: GEIH 2025-IV, DANE."
nota_violin <- "Violín: densidad ponderada. Caja: cuartiles ponderados. Rombo: media ponderada."
tema <- theme_minimal(base_size = 11) +
theme(plot.title = element_text(face = "bold", size = 12),
plot.caption = element_text(size = 9, color = "grey40", hjust = 0),
strip.text = element_text(face = "bold", size = 11),
panel.grid.minor = element_blank(), legend.position = "none")
carpeta <- file.path(ruta, "figuras") # figuras y tablas de salida
dir.create(carpeta, showWarnings = FALSE, recursive = TRUE)5. Verificación de códigos
table(cauca$MES) # tres meses del trimestre##
## 10 11 12
## 2225 2328 2294table(sexo = cauca$bsex) # 1 = hombre, 0 = mujer## sexo
## 0 1
## 3727 3120table(zona = cauca$area) # 1 = urbano, 0 = rural## zona
## 0 1
## 1191 5656table(pension = cauca$P6920, useNA = "ifany") # 1 = cotiza, 2 = no cotiza, 3 = pensionado## pension
## 1 2 3 <NA>
## 990 1997 44 3816table(posicion = cauca$P6430, useNA = "ifany") # posición ocupacional## posicion
## 1 2 3 4 5 6 7 8 <NA>
## 1027 180 84 1578 64 75 19 4 3816c(edad_P6040 = "P6040" %in% names(cauca), vivienda = "DIRECTORIO" %in% names(cauca))## edad_P6040 vivienda
## TRUE TRUEParte II. Muestra y análisis descriptivo
6. Base de trabajo y flujo de depuración
base0 <- cauca %>%
transmute(Inglab = as.numeric(Inglab),
Horatsem = as.numeric(Horatsem),
educ = as.numeric(anoeduca),
exper = as.numeric(exper),
edad = as.numeric(P6040),
estrato = as.numeric(estrato),
niv = as.numeric(NivEduc),
w = as.numeric(FEX_C18) / 3, # factor mensual llevado al trimestre (descriptivos)
P6920 = as.numeric(P6920),
P6430 = as.numeric(P6430),
quiere_mas = as.numeric(P7090),
etnia = as.numeric(etnia),
bsex = as.numeric(bsex),
area = as.numeric(area),
vivienda = as.numeric(DIRECTORIO)) %>%
mutate(sexo = factor(bsex, levels = c(1, 0), labels = c("Hombres", "Mujeres")),
zona = factor(area, levels = c(1, 0), labels = c("Urbano", "Rural")),
mujer = as.integer(bsex == 0),
rural = as.integer(area == 0),
no_cotiza = case_when(P6920 == 2 ~ 1, P6920 == 1 ~ 0, TRUE ~ NA_real_),
nivel = factor(case_when(niv <= 2 ~ "Primaria o menos", niv == 3 ~ "Secundaria",
niv == 4 ~ "Media", niv == 5 ~ "Técnica o tecnológica",
niv == 6 ~ "Universitaria", niv == 7 ~ "Posgrado"),
levels = c("Primaria o menos", "Secundaria", "Media",
"Técnica o tecnológica", "Universitaria", "Posgrado")),
zona_estrato = case_when(rural == 1 ~ "Rural",
estrato == 1 ~ "Urbano E1", estrato == 2 ~ "Urbano E2",
estrato == 3 ~ "Urbano E3", estrato %in% 4:6 ~ "Urbano E4-6",
TRUE ~ NA_character_)) # urbanos sin estrato
# ---- Depuración paso a paso ----
f1 <- base0 %>% filter(Inglab > 0)
f2 <- f1 %>% filter(!is.na(educ), !is.na(exper), !is.na(nivel))
f3 <- f2 %>% filter(edad >= 15, edad <= 65, exper >= 0)
f4 <- f3 %>% filter(Horatsem >= 5)
f5 <- f4 %>% filter(!is.na(zona_estrato))
df <- f5 %>%
mutate(horas_mes = Horatsem * sem_mes,
lingreso = log(Inglab), # variable dependiente principal
lhoras = log(horas_mes), # robustez: control por horas mensuales
lw = log(Inglab / horas_mes), # robustez: ingreso por hora
exper2 = exper^2 / 100,
empleo = factor(no_cotiza, levels = c(0, 1),
labels = c("Cotiza a pensión", "No cotiza a pensión")),
zona_estrato = factor(zona_estrato, levels = c("Urbano E1", "Urbano E2", "Urbano E3",
"Urbano E4-6", "Rural")),
w_norm = w / mean(w)) # pesos de media 1 (solo para robustez)
dm <- df # una sola base para descriptivos y modelos
flujo <- c("Registros del trimestre" = nrow(base0),
"Ocupados con ingreso laboral positivo" = nrow(f1),
"Con educación, experiencia y nivel educativo" = nrow(f2),
"Edad de 15 a 65 años y experiencia no negativa" = nrow(f3),
"Al menos 5 horas de trabajo semanal" = nrow(f4),
"Con dato de estrato (urbanos)" = nrow(f5))
tabla_flujo <- data.frame(paso = names(flujo), observaciones = flujo,
excluidas = c(NA, NA, -diff(flujo)[-1]), row.names = NULL)
tabla_flujo## paso observaciones excluidas
## 1 Registros del trimestre 6847 NA
## 2 Ocupados con ingreso laboral positivo 2875 NA
## 3 Con educación, experiencia y nivel educativo 2875 0
## 4 Edad de 15 a 65 años y experiencia no negativa 2682 193
## 5 Al menos 5 horas de trabajo semanal 2668 14
## 6 Con dato de estrato (urbanos) 2654 14# Detalle de exclusiones del paso 4
c(mayores_65 = sum(f2$edad > 65), menores_15 = sum(f2$edad < 15),
exper_negativa = sum(f2$exper < 0 & f2$edad <= 65 & f2$edad >= 15))## mayores_65 menores_15 exper_negativa
## 189 0 4stopifnot(all(is.finite(df$w)), all(df$w > 0), all(is.finite(df$lingreso)))
cat("Total ponderado:", fi(round(sum(df$w))), "ocupados\n")## Total ponderado: 566.107 ocupadostable(df$zona_estrato)##
## Urbano E1 Urbano E2 Urbano E3 Urbano E4-6 Rural
## 743 867 438 194 4127. Cifras de contexto (ponderadas)
indicadores <- function(d) c(
"Ocupados (n)" = nrow(d),
"Total ponderado" = sum(d$w),
"% no cotiza a pensión" = 100 * weighted.mean(d$no_cotiza, d$w, na.rm = TRUE),
"% rural" = 100 * weighted.mean(d$rural, d$w),
"% mujeres" = 100 * weighted.mean(d$mujer, d$w),
"% grupo étnico" = 100 * weighted.mean(d$etnia == 1, d$w, na.rm = TRUE),
"% bajo el mínimo" = 100 * weighted.mean(d$Inglab < smmlv, d$w),
"% bajo el mínimo por hora" = 100 * weighted.mean(d$Inglab / (d$Horatsem * sem_mes) < smmlv_hora, d$w),
"Ingreso mediano" = wq(d$Inglab, d$w, 0.5),
"Años de educación" = weighted.mean(d$educ, d$w))
round(data.frame(todos_ocupados = indicadores(f1), muestra_analisis = indicadores(df)), 1)## todos_ocupados muestra_analisis
## Ocupados (n) 2875.0 2654.0
## Total ponderado 612050.7 566107.3
## % no cotiza a pensión 81.2 80.6
## % rural 58.6 58.7
## % mujeres 35.7 36.0
## % grupo étnico 46.9 47.6
## % bajo el mínimo 72.2 71.5
## % bajo el mínimo por hora 70.9 70.2
## Ingreso mediano 825000.0 850000.0
## Años de educación 8.4 8.6# Proporción rural en la muestra y en la población (relevante al no ponderar los modelos)
c(pct_muestra = 100 * mean(df$rural), pct_poblacion = 100 * weighted.mean(df$rural, df$w))## pct_muestra pct_poblacion
## 15.52374 58.68223# Bajo el mínimo por zona
df %>% group_by(zona) %>% summarise(bajo_min = 100 * weighted.mean(Inglab < smmlv, w))## # A tibble: 2 × 2
## zona bajo_min
## <fct> <dbl>
## 1 Urbano 53.9
## 2 Rural 83.9# Estrato en la zona rural (justifica la variable zona-estrato)
f4 %>% filter(rural == 1) %>%
summarise(estrato1 = 100 * weighted.mean(estrato == 1, w, na.rm = TRUE),
sin_estrato = 100 * weighted.mean(is.na(estrato), w))## # A tibble: 1 × 2
## estrato1 sin_estrato
## <dbl> <dbl>
## 1 99.5 12.38. Comparación por sexo (Anexo)
Medias ponderadas y diferencia mujeres - hombres con errores robustos (HC1).
vars_sexo <- c(lingreso = "Log del ingreso mensual", lw = "Log del ingreso por hora",
horas_mes = "Horas mensuales", educ = "Años de educación",
exper = "Experiencia potencial", rural = "Zona rural (proporción)",
no_cotiza = "No cotiza a pensión (proporción)")
prueba_sexo <- function(v) {
d <- df %>% filter(!is.na(.data[[v]]))
m <- lm(d[[v]] ~ d$mujer, weights = d$w)
b <- coef(m); V <- vcovHC(m, type = "HC1")
c(hombres = unname(b[1]), mujeres = unname(b[1] + b[2]), diferencia = unname(b[2]),
p_valor = unname(2 * pnorm(-abs(b[2] / sqrt(V[2, 2])))))
}
tabla_sexo <- t(sapply(names(vars_sexo), prueba_sexo))
rownames(tabla_sexo) <- vars_sexo
round(tabla_sexo, 3)## hombres mujeres diferencia p_valor
## Log del ingreso mensual 13.695 13.505 -0.190 0.002
## Log del ingreso por hora 8.618 8.622 0.003 0.954
## Horas mensuales 171.096 148.273 -22.823 0.000
## Años de educación 8.015 9.708 1.693 0.000
## Experiencia potencial 24.877 22.986 -1.891 0.079
## Zona rural (proporción) 0.644 0.485 -0.159 0.000
## No cotiza a pensión (proporción) 0.825 0.772 -0.054 0.022# Brecha bruta en % (mensual y por hora)
c(mensual = pct(tabla_sexo[1, "diferencia"]), por_hora = pct(tabla_sexo[2, "diferencia"]))## mensual por_hora
## -17.2859658 0.3159989table(df$sexo)##
## Hombres Mujeres
## 1419 12359. Ocupados según posición ocupacional
Justifica hablar de ingreso laboral y no de salario.
tab_pos <- df %>%
mutate(grupo = case_when(P6430 %in% c(1, 2, 3, 8) ~ "Asalariados",
P6430 == 4 ~ "Cuenta propia", P6430 == 5 ~ "Patrones", TRUE ~ "Otro")) %>%
group_by(grupo) %>%
summarise(n = n(), pct_ocup = 100 * sum(w) / sum(df$w),
no_cotiza = 100 * weighted.mean(no_cotiza, w, na.rm = TRUE),
mediana = wq(Inglab, w, 0.5),
bajo_min = 100 * weighted.mean(Inglab < smmlv, w),
horas_mes = weighted.mean(horas_mes, w),
educ = weighted.mean(educ, w), mujeres = 100 * weighted.mean(mujer, w))
tab_pos## # A tibble: 4 × 9
## grupo n pct_ocup no_cotiza mediana bajo_min horas_mes educ mujeres
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Asalariados 1229 30.0 50.0 1423500 42.4 180. 11.0 48.1
## 2 Cuenta prop… 1356 65.4 93.6 600000 86.3 152. 7.61 31.9
## 3 Otro 18 2.65 96.7 1200000 69.4 194. 6.63 11.6
## 4 Patrones 51 1.94 92.0 2000000 25.4 210. 8.69 20.210. Gráficas descriptivas
10.1 Ingreso mensual por sexo: niveles frente a logaritmo
lab_niv <- "Ingreso mensual (pesos)"; lab_log <- "Log del ingreso mensual"
largo <- bind_rows(df %>% transmute(sexo, w, escala = lab_niv, valor = Inglab),
df %>% transmute(sexo, w, escala = lab_log, valor = lingreso)) %>%
mutate(escala = factor(escala, levels = c(lab_niv, lab_log)))
cajas_s <- cajas_pond(largo, valor, escala, sexo) %>%
mutate(etq_media = ifelse(escala == lab_niv, paste0("$", num(media / 1000, 0), " mil"), num(media, 2)))
largo_vis <- largo %>% group_by(escala) %>%
filter(valor >= wq(valor, w, 0.01), valor <= wq(valor, w, 0.99)) %>% ungroup()
ref_min <- data.frame(escala = factor(c(lab_niv, lab_log), levels = c(lab_niv, lab_log)),
valor = c(smmlv, log(smmlv)))
g_ingreso_sexo <- ggplot() +
geom_hline(data = ref_min, aes(yintercept = valor), linetype = "dashed", color = "gray40") +
violin_caja(largo_vis, cajas_s, sexo, valor, colores_sexo) +
facet_wrap(~ escala, scales = "free_y") +
scale_y_continuous(labels = label_number(big.mark = ".", decimal.mark = ",")) +
labs(title = "El logaritmo reduce la asimetría del ingreso", x = NULL, y = NULL,
caption = paste("Violín entre los percentiles 1 y 99. Caja: cuartiles ponderados. Rombo: media ponderada.",
"\nLínea punteada: salario mínimo 2025.", fuente)) +
tema
g_ingreso_sexo
10.2 Horas mensuales según cotización a pensión y sexo
df_cot <- df %>% filter(!is.na(empleo))
cajas_h <- cajas_pond(df_cot, horas_mes, sexo, empleo) %>%
left_join(df_cot %>% group_by(sexo, empleo) %>%
summarise(corta = weighted.mean(Horatsem < 20, w), # < 85,7 h al mes
quiere = weighted.mean(quiere_mas == 1, w, na.rm = TRUE), .groups = "drop"),
by = c("sexo", "empleo")) %>%
mutate(etq_media = paste0(num(media), " h"),
etiqueta = paste0("Menos de 86 h al mes: ", num(100 * corta, 0), "%",
"\nQuiere más horas: ", num(100 * quiere, 0), "%"))
cajas_h %>% select(sexo, empleo, mediana, media, corta, quiere)## # A tibble: 4 × 6
## sexo empleo mediana media corta quiere
## <fct> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Hombres Cotiza a pensión 189. 195. 0.00567 0.0850
## 2 Hombres No cotiza a pensión 171. 166. 0.0698 0.137
## 3 Mujeres Cotiza a pensión 189. 188. 0.00104 0.0741
## 4 Mujeres No cotiza a pensión 137. 137. 0.232 0.216g_horas <- ggplot() +
geom_hline(yintercept = jornada_mes, linetype = "dashed", color = "gray40") +
violin_caja(df_cot, cajas_h, empleo, horas_mes, colores_emp) +
geom_text(data = cajas_h, aes(x = empleo, y = 400, label = etiqueta), size = 3, lineheight = 0.9) +
facet_wrap(~ sexo) +
coord_cartesian(ylim = c(0, 412)) +
labs(title = "Horas trabajadas al mes según cotización a pensión y sexo", x = NULL, y = "Horas al mes",
caption = paste(nota_violin, "Línea punteada: jornada legal (44 h semanales = 188,5 h al mes).", fuente)) +
tema
g_horas
10.3 Perfil del ingreso por tramos de experiencia
Promedios ponderados sin controles; las etiquetas muestran la educación media de cada tramo (efecto cohorte).
perfil <- df %>%
mutate(tramo = cut(exper, c(-1, 10, 20, 30, 40, 50, 60),
labels = c("0-10", "11-20", "21-30", "31-40", "41-50", "51-60"))) %>%
group_by(sexo, tramo) %>%
summarise(l = weighted.mean(lingreso, w), educ_m = weighted.mean(educ, w), n = n(), .groups = "drop")
perfil## # A tibble: 12 × 5
## sexo tramo l educ_m n
## <fct> <fct> <dbl> <dbl> <int>
## 1 Hombres 0-10 13.7 10.8 304
## 2 Hombres 11-20 13.9 9.44 334
## 3 Hombres 21-30 13.9 8.63 296
## 4 Hombres 31-40 13.6 6.25 265
## 5 Hombres 41-50 13.4 4.78 157
## 6 Hombres 51-60 13.1 2.34 63
## 7 Mujeres 0-10 13.5 12.3 266
## 8 Mujeres 11-20 13.6 11.3 299
## 9 Mujeres 21-30 13.7 9.91 284
## 10 Mujeres 31-40 13.5 8.02 228
## 11 Mujeres 41-50 13.3 5.17 120
## 12 Mujeres 51-60 12.8 3.54 38g_exp_perfil <- ggplot(perfil, aes(tramo, l, color = sexo, group = sexo)) +
geom_line(linewidth = 1.2) + geom_point(size = 3) +
geom_text(aes(label = paste0(num(educ_m, 1), " a"), vjust = ifelse(sexo == "Hombres", -1.2, 2)),
size = 3, show.legend = FALSE) +
scale_color_manual(values = colores_sexo) +
labs(title = "Perfil del ingreso mensual por tramos de experiencia",
x = "Años de experiencia potencial", y = "Log del ingreso mensual", color = NULL,
caption = paste("Promedios ponderados sin controles. Etiquetas: años promedio de educación.", fuente)) +
tema + theme(legend.position = "top")
g_exp_perfil
10.4 Educación e ingreso
medias_edu <- df %>% group_by(sexo, educ) %>%
summarise(n = n(), l = weighted.mean(lingreso, w), .groups = "drop") %>% filter(n >= 10)
g_educ_ingreso <- ggplot(df, aes(educ, lingreso, color = sexo)) +
geom_smooth(aes(weight = w), method = "lm", formula = y ~ x, se = FALSE, linetype = "dashed") +
geom_line(data = medias_edu, aes(y = l), linewidth = 1.1) +
geom_point(data = medias_edu, aes(y = l, size = n), shape = 21, fill = "white", stroke = 1.2) +
scale_size_area(max_size = 5, guide = "none") +
scale_color_manual(values = colores_sexo) +
scale_x_continuous(breaks = c(0, 5, 9, 11, 16, 21)) +
coord_cartesian(ylim = range(medias_edu$l) + c(-0.3, 0.3)) +
labs(title = "El ingreso crece con la educación, sobre todo con el título universitario",
x = "Años de educación", y = "Log del ingreso mensual", color = NULL,
caption = paste("Línea continua: media ponderada por año (grupos con al menos 10 personas).",
"Línea punteada: ajuste lineal ponderado.", fuente)) +
tema + theme(legend.position = "top")
g_educ_ingreso
10.5 Ingreso por zona y sexo
cajas_z <- cajas_pond(df, Inglab, sexo, zona, log = TRUE) %>%
left_join(df %>% group_by(sexo, zona) %>%
summarise(bajo_min = weighted.mean(Inglab < smmlv, w),
nc = weighted.mean(no_cotiza, w, na.rm = TRUE), .groups = "drop"),
by = c("sexo", "zona")) %>%
mutate(etq_media = paste0("$", num(media / 1000, 0), " mil"),
etiqueta = paste0("Bajo el mínimo: ", num(100 * bajo_min, 0), "%\nNo cotiza: ", num(100 * nc, 0), "%"))
cajas_z %>% select(sexo, zona, mediana, media, bajo_min, nc)## # A tibble: 4 × 6
## sexo zona mediana media bajo_min nc
## <fct> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Hombres Urbano 1423500 1311400. 0.498 0.678
## 2 Hombres Rural 800000 713515. 0.827 0.906
## 3 Mujeres Urbano 1260000 1071913. 0.591 0.632
## 4 Mujeres Rural 400000 489854. 0.869 0.919g_ingreso_zona <- ggplot() +
geom_hline(yintercept = smmlv, linetype = "dashed", color = "gray40") +
violin_caja(df, cajas_z, zona, Inglab, colores_zona) +
geom_text(data = cajas_z, aes(x = zona, y = 4e7, label = etiqueta), size = 3, lineheight = 0.9) +
facet_wrap(~ sexo) +
scale_y_log10(breaks = c(1e5, 2.5e5, 5e5, 1e6, 2.5e6, 5e6, 1e7), labels = pesos) +
coord_cartesian(ylim = c(3e4, 6e7)) +
labs(title = "Ingreso laboral mensual por zona y sexo (escala logarítmica)", x = NULL,
y = "Ingreso laboral mensual",
caption = paste("Rombo: media geométrica. Línea punteada: salario mínimo 2025.", fuente)) +
tema
g_ingreso_zona
10.6 Ingreso por zona y estrato, centrado en la media, por sexo
l_total <- weighted.mean(df$lingreso, df$w)
centrado <- bind_rows(
df %>% group_by(zona_estrato) %>% summarise(n = n(), l = weighted.mean(lingreso, w), .groups = "drop") %>%
mutate(grupo = "Total"),
df %>% group_by(zona_estrato, grupo = as.character(sexo)) %>%
summarise(n = n(), l = weighted.mean(lingreso, w), .groups = "drop")) %>%
mutate(grupo = factor(grupo, levels = c("Total", "Hombres", "Mujeres")), dif_pct = pct(l - l_total))
centrado## # A tibble: 15 × 5
## zona_estrato n l grupo dif_pct
## <fct> <int> <dbl> <fct> <dbl>
## 1 Urbano E1 743 13.8 Total 15.9
## 2 Urbano E2 867 14.1 Total 54.1
## 3 Urbano E3 438 14.4 Total 109.
## 4 Urbano E4-6 194 15.0 Total 301.
## 5 Rural 412 13.4 Total -22.9
## 6 Urbano E1 383 13.9 Hombres 32.6
## 7 Urbano E1 360 13.6 Mujeres -1.69
## 8 Urbano E2 451 14.1 Hombres 67.9
## 9 Urbano E2 416 13.9 Mujeres 36.8
## 10 Urbano E3 217 14.4 Hombres 110.
## 11 Urbano E3 221 14.4 Mujeres 108.
## 12 Urbano E4-6 89 15.0 Hombres 304.
## 13 Urbano E4-6 105 15.0 Mujeres 299.
## 14 Rural 279 13.5 Hombres -13.8
## 15 Rural 133 13.1 Mujeres -40.8g_estrato <- ggplot(centrado, aes(zona_estrato, dif_pct, fill = grupo)) +
geom_hline(yintercept = 0, color = "gray30") +
geom_col(position = position_dodge(width = 0.8), width = 0.75) +
geom_text(aes(label = paste0(ifelse(dif_pct > 0, "+", ""), num(dif_pct, 0), "%"),
vjust = ifelse(dif_pct > 0, -0.4, 1.3)),
position = position_dodge(width = 0.8), size = 3) +
scale_fill_manual(values = c(Total = "gray55", colores_sexo)) +
scale_y_continuous(labels = function(x) paste0(x, "%"), expand = expansion(mult = 0.12)) +
labs(title = "El ingreso crece con el estrato urbano y es menor en la zona rural",
subtitle = paste0("Diferencia frente a la media geométrica total (", pesos(round(exp(l_total), -3)), " al mes)"),
x = NULL, y = "Diferencia respecto a la media", fill = NULL,
caption = paste("Medias geométricas ponderadas del ingreso mensual.", fuente)) +
tema + theme(legend.position = "top")
g_estrato
10.7 Exportar figuras
figuras <- list(ingreso_sexo = g_ingreso_sexo, horas_cotizacion = g_horas,
exp_perfil = g_exp_perfil, educ_ingreso = g_educ_ingreso,
ingreso_zona = g_ingreso_zona, estrato_centrado = g_estrato)
for (nombre in names(figuras))
ggsave(file.path(carpeta, paste0(nombre, ".png")),
figuras[[nombre]] + labs(title = NULL, subtitle = NULL),
width = 6.5, height = 4.2, units = "in", dpi = 300, bg = "white")
list.files(carpeta)## [1] "composicion.jpg" "educ_ingreso.jpg" "educ_ingreso.pdf"
## [4] "educ_ingreso.png" "estrato_centrado.jpg" "estrato_centrado.pdf"
## [7] "estrato_centrado.png" "exp_perfil.jpg" "exp_perfil.pdf"
## [10] "exp_perfil.png" "horas_cotizacion.png" "horas_informalidad.jpg"
## [13] "horas_informalidad.pdf" "horas_informalidad.png" "ingreso_sexo.jpg"
## [16] "ingreso_sexo.pdf" "ingreso_sexo.png" "ingreso_zona.jpg"
## [19] "ingreso_zona.pdf" "ingreso_zona.png" "residuales.png"
## [22] "residuales_m3.png" "residuales_m5.jpg" "residuales_niveles.png"
## [25] "tabla_brechas.tex" "tabla_correccion.tex" "tabla_diagnostico.tex"
## [28] "tabla_flujo.tex" "tabla_hipotesis.tex" "tabla_m4.tex"
## [31] "tabla_modelos.tex" "tabla_posicion.tex" "tabla_resultados.tex"
## [34] "tabla_robustez.tex" "tabla_sexo.tex"Parte III. Estimación
11. Modelos principales
Variable dependiente: log del ingreso laboral mensual. MCO sin ponderar: los coeficientes describen la relación en la muestra. La pertinencia de no ponderar se evalúa con la prueba de DuMouchel y Duncan (sección 21) y la estimación ponderada se reporta como robustez.
- (1) Años: retorno promedio a la educación y brecha por sexo.
- (2) Años x sexo: retorno a la educación por sexo.
- (3) Niveles x sexo: primas por nivel y brecha por sexo en cada nivel.
m1 <- lm(lingreso ~ educ + exper + exper2 + mujer + zona_estrato, data = dm)
m2 <- lm(lingreso ~ educ * mujer + exper + exper2 + zona_estrato, data = dm)
m3 <- lm(lingreso ~ nivel * mujer + exper + exper2 + zona_estrato, data = dm)
m_niv <- lm(lingreso ~ nivel + mujer + exper + exper2 + zona_estrato, data = dm) # auxiliar
modelos <- list("(1) Años" = m1, "(2) Años x sexo" = m2, "(3) Niveles x sexo" = m3)
stopifnot(length(unique(sapply(modelos, nobs))) == 1) # misma muestra en los tres
for (k in names(modelos)) { cat("\n=====", k, "(errores convencionales) =====\n"); print(coeftest(modelos[[k]])) }##
## ===== (1) Años (errores convencionales) =====
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.6506715 0.0685465 184.5560 < 2.2e-16 ***
## educ 0.0908802 0.0043362 20.9585 < 2.2e-16 ***
## exper 0.0253859 0.0034730 7.3094 3.536e-13 ***
## exper2 -0.0397707 0.0066387 -5.9908 2.372e-09 ***
## mujer -0.2304170 0.0271081 -8.4999 < 2.2e-16 ***
## zona_estratoUrbano E2 0.1379201 0.0348157 3.9614 7.647e-05 ***
## zona_estratoUrbano E3 0.2502004 0.0435509 5.7450 1.024e-08 ***
## zona_estratoUrbano E4-6 0.7684681 0.0594291 12.9308 < 2.2e-16 ***
## zona_estratoRural -0.1439131 0.0436099 -3.3000 0.0009796 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## ===== (2) Años x sexo (errores convencionales) =====
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.8455034 0.0752095 170.7963 < 2.2e-16 ***
## educ 0.0725139 0.0052559 13.7968 < 2.2e-16 ***
## mujer -0.6461070 0.0733059 -8.8139 < 2.2e-16 ***
## exper 0.0259575 0.0034508 7.5222 7.338e-14 ***
## exper2 -0.0413048 0.0065985 -6.2597 4.486e-10 ***
## zona_estratoUrbano E2 0.1354999 0.0345824 3.9182 9.148e-05 ***
## zona_estratoUrbano E3 0.2496889 0.0432563 5.7723 8.730e-09 ***
## zona_estratoUrbano E4-6 0.7599767 0.0590433 12.8715 < 2.2e-16 ***
## zona_estratoRural -0.1695728 0.0435187 -3.8965 1.000e-04 ***
## educ:mujer 0.0378960 0.0062158 6.0968 1.241e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## ===== (3) Niveles x sexo (errores convencionales) =====
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.2753995 0.0649429 204.4166 < 2.2e-16 ***
## nivelSecundaria 0.1486332 0.0684051 2.1728 0.029881 *
## nivelMedia 0.2304166 0.0553252 4.1648 3.217e-05 ***
## nivelTécnica o tecnológica 0.4576772 0.0694075 6.5941 5.152e-11 ***
## nivelUniversitaria 0.6659750 0.0682204 9.7621 < 2.2e-16 ***
## nivelPosgrado 1.2920630 0.0962340 13.4263 < 2.2e-16 ***
## mujer -0.3570682 0.0605661 -5.8955 4.212e-09 ***
## exper 0.0336860 0.0033994 9.9094 < 2.2e-16 ***
## exper2 -0.0623461 0.0064613 -9.6491 < 2.2e-16 ***
## zona_estratoUrbano E2 0.1459337 0.0337174 4.3281 1.560e-05 ***
## zona_estratoUrbano E3 0.2091503 0.0424465 4.9274 8.850e-07 ***
## zona_estratoUrbano E4-6 0.5872606 0.0594237 9.8826 < 2.2e-16 ***
## zona_estratoRural -0.2299405 0.0424849 -5.4123 6.783e-08 ***
## nivelSecundaria:mujer -0.1442318 0.1067446 -1.3512 0.176752
## nivelMedia:mujer 0.0290216 0.0768291 0.3777 0.705653
## nivelTécnica o tecnológica:mujer 0.0730179 0.0892411 0.8182 0.413311
## nivelUniversitaria:mujer 0.4093061 0.0858214 4.7693 1.950e-06 ***
## nivelPosgrado:mujer 0.3072643 0.1166299 2.6345 0.008475 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Parte IV. Heterocedasticidad: diagnóstico, fuente y corrección
12. Diagnóstico
pruebas_het <- t(sapply(modelos, function(m) {
bp <- bptest(m) # Breusch-Pagan (Koenker)
wh <- bptest(m, ~ fitted(m) + I(fitted(m)^2)) # White (forma especial)
c(BP = unname(bp$statistic), BP_gl = unname(bp$parameter), BP_p = unname(bp$p.value),
White = unname(wh$statistic), White_p = unname(wh$p.value))
}))
round(pruebas_het, 3)## BP BP_gl BP_p White White_p
## (1) Años 62.186 8 0 15.501 0
## (2) Años x sexo 65.664 9 0 19.708 0
## (3) Niveles x sexo 89.573 17 0 31.099 013. Fuente: ¿qué regresores generan la heterocedasticidad?
e2 <- residuals(m2)^2
aux <- lm(e2 ~ model.matrix(m2)[, -1])
coeftest(aux, vcov. = vcovHC(aux, type = "HC1"))##
## t test of coefficients:
##
## Estimate Std. Error t value
## (Intercept) 0.3373674 0.0960673 3.5118
## model.matrix(m2)[, -1]educ 0.0168910 0.0064315 2.6263
## model.matrix(m2)[, -1]mujer 0.3091498 0.1237382 2.4984
## model.matrix(m2)[, -1]exper -0.0133850 0.0050058 -2.6739
## model.matrix(m2)[, -1]exper2 0.0365019 0.0106661 3.4223
## model.matrix(m2)[, -1]zona_estratoUrbano E2 -0.1514026 0.0516583 -2.9308
## model.matrix(m2)[, -1]zona_estratoUrbano E3 -0.1733308 0.0589688 -2.9394
## model.matrix(m2)[, -1]zona_estratoUrbano E4-6 -0.2487999 0.0827472 -3.0067
## model.matrix(m2)[, -1]zona_estratoRural 0.0404704 0.0605989 0.6678
## model.matrix(m2)[, -1]educ:mujer -0.0130310 0.0095519 -1.3642
## Pr(>|t|)
## (Intercept) 0.0004526 ***
## model.matrix(m2)[, -1]educ 0.0086823 **
## model.matrix(m2)[, -1]mujer 0.0125352 *
## model.matrix(m2)[, -1]exper 0.0075437 **
## model.matrix(m2)[, -1]exper2 0.0006305 ***
## model.matrix(m2)[, -1]zona_estratoUrbano E2 0.0034094 **
## model.matrix(m2)[, -1]zona_estratoUrbano E3 0.0033174 **
## model.matrix(m2)[, -1]zona_estratoUrbano E4-6 0.0026655 **
## model.matrix(m2)[, -1]zona_estratoRural 0.5042931
## model.matrix(m2)[, -1]educ:mujer 0.1726119
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Coeficiente positivo: la varianza del error aumenta con esa variable.14. Corrección 1: errores robustos (inferencia)
vcovs <- lapply(modelos, vcovHC, type = "HC1")
cbind(coef = coef(m2), ee_clasico = sqrt(diag(vcov(m2))), ee_robusto = sqrt(diag(vcovs[[2]])),
razon = sqrt(diag(vcovs[[2]])) / sqrt(diag(vcov(m2))))## coef ee_clasico ee_robusto razon
## (Intercept) 12.84550336 0.075209480 0.076291681 1.0143892
## educ 0.07251386 0.005255856 0.005167478 0.9831847
## mujer -0.64610697 0.073305857 0.079914679 1.0901541
## exper 0.02595754 0.003450808 0.003722002 1.0785886
## exper2 -0.04130478 0.006598532 0.007256635 1.0997348
## zona_estratoUrbano E2 0.13549995 0.034582370 0.035636962 1.0304951
## zona_estratoUrbano E3 0.24968890 0.043256251 0.043971403 1.0165329
## zona_estratoUrbano E4-6 0.75997674 0.059043341 0.056140110 0.9508288
## zona_estratoRural -0.16957283 0.043518742 0.046784423 1.0750408
## educ:mujer 0.03789601 0.006215770 0.006520169 1.048971915. Corrección 2: MCGF (eficiencia)
Función de varianza exponencial en los regresores (Wooldridge, cap. 8): se estima log(u²) sobre los regresores y se reestima el modelo por mínimos cuadrados ponderados con pesos 1 / ĥ.
mcgf <- function(m) {
d <- model.frame(m)
X <- model.matrix(m)
aux <- lm(log(pmax(residuals(m)^2, 1e-12)) ~ X - 1)
d$peso_h <- 1 / exp(fitted(aux)) # peso = 1 / varianza estimada del error
lm(formula(m), data = d, weights = peso_h)
}
gls <- lapply(modelos, mcgf)
vcovs_g <- lapply(gls, vcovHC, type = "HC1")
# ¿Queda heterocedasticidad en el modelo transformado? (aproximación)
bp_transformado <- function(m, g) {
e_s <- residuals(g) * sqrt(weights(g) / mean(weights(g)))
X <- model.matrix(m)[, -1]
BP <- length(e_s) * summary(lm(e_s^2 ~ X))$r.squared
c(BP = BP, p_valor = pchisq(BP, df = ncol(X), lower.tail = FALSE))
}
bp_mcgf <- t(mapply(bp_transformado, modelos, gls))
round(bp_mcgf, 4)## BP p_valor
## (1) Años 18.3263 0.0189
## (2) Años x sexo 16.2816 0.0612
## (3) Niveles x sexo 29.2333 0.0324# Comparación de coeficientes y precisión en el modelo (2)
round(data.frame(MCO = coef(m2), ee_MCO = sqrt(diag(vcovs[[2]])),
MCGF = coef(gls[[2]]), ee_MCGF = sqrt(diag(vcovs_g[[2]])),
razon_ee = sqrt(diag(vcovs_g[[2]])) / sqrt(diag(vcovs[[2]]))), 4)## MCO ee_MCO MCGF ee_MCGF razon_ee
## (Intercept) 12.8455 0.0763 12.8958 0.0757 0.9923
## educ 0.0725 0.0052 0.0694 0.0051 0.9804
## mujer -0.6461 0.0799 -0.7861 0.0794 0.9931
## exper 0.0260 0.0037 0.0270 0.0037 1.0034
## exper2 -0.0413 0.0073 -0.0459 0.0073 1.0099
## zona_estratoUrbano E2 0.1355 0.0356 0.1288 0.0339 0.9526
## zona_estratoUrbano E3 0.2497 0.0440 0.2283 0.0420 0.9561
## zona_estratoUrbano E4-6 0.7600 0.0561 0.7421 0.0546 0.9717
## zona_estratoRural -0.1696 0.0468 -0.1624 0.0459 0.9820
## educ:mujer 0.0379 0.0065 0.0500 0.0064 0.9844Parte V. Resultados
16. Cantidades de interés
V1 <- vcovs[[1]]; V2 <- vcovs[[2]]; V3 <- vcovs[[3]]
# Brecha por sexo y su descomposición (bruta y con controles parciales)
brecha <- c(bruta = coef(lm(lingreso ~ mujer, data = dm))[["mujer"]],
con_horas = coef(lm(lingreso ~ mujer + lhoras, data = dm))[["mujer"]],
con_educ = coef(lm(lingreso ~ mujer + educ, data = dm))[["mujer"]],
condicionada = coef(m1)[["mujer"]])
rbind(log_puntos = brecha, porcentaje = pct(brecha))## bruta con_horas con_educ condicionada
## log_puntos -0.0910402 0.02950864 -0.2094202 -0.230417
## porcentaje -8.7018992 2.99483368 -18.8945673 -20.579766# Retorno promedio y brecha condicionada (modelo 1)
rbind(retorno = c(pct = pct(coef(m1)[["educ"]]), comb(m1, V1, c(educ = 1))),
brecha = c(pct = pct(coef(m1)[["mujer"]]), comb(m1, V1, c(mujer = 1))))## pct estimacion error p_valor
## retorno 9.513778 0.09088018 0.004628474 3.063185e-80
## brecha -20.579766 -0.23041702 0.027399238 6.631864e-17# Retorno por sexo (modelo 2)
r_h <- comb(m2, V2, c(educ = 1)); r_m <- comb(m2, V2, c(educ = 1, "educ:mujer" = 1))
rbind(hombres = c(pct = pct(r_h[["estimacion"]]), r_h),
mujeres = c(pct = pct(r_m[["estimacion"]]), r_m),
diferencia = c(pct = NA, comb(m2, V2, c("educ:mujer" = 1))))## pct estimacion error p_valor
## hombres 7.520771 0.07251386 0.005167478 3.353272e-43
## mujeres 11.673570 0.11040988 0.006118138 9.525536e-69
## diferencia NA 0.03789601 0.006520169 6.908358e-09# Brecha según años de educación (modelo 2)
t(sapply(c(5, 9, 11, 16), function(s) {
r <- comb(m2, V2, c(mujer = 1, "educ:mujer" = s)); c(anios = s, brecha_pct = pct(r[["estimacion"]]), r)
}))## anios brecha_pct estimacion error p_valor
## [1,] 5 -36.658338 -0.45662691 0.05048521 2.830819e-19
## [2,] 9 -26.290821 -0.30504286 0.03176681 1.740294e-21
## [3,] 11 -20.487093 -0.22925083 0.02737642 8.908287e-17
## [4,] 16 -3.899029 -0.03977077 0.03985291 3.184003e-01# Primas por nivel para hombres y brecha por sexo en cada nivel (modelo 3)
niveles_k <- levels(dm$nivel)[-1]
primas <- t(sapply(niveles_k, function(k) {
r <- comb(m3, V3, setNames(1, paste0("nivel", k))); c(prima_pct = pct(r[["estimacion"]]), r)
}))
round(primas, 4)## prima_pct estimacion error p_valor
## Secundaria 16.0247 0.1486 0.0649 0.0222
## Media 25.9124 0.2304 0.0545 0.0000
## Técnica o tecnológica 58.0399 0.4577 0.0637 0.0000
## Universitaria 94.6387 0.6660 0.0664 0.0000
## Posgrado 264.0289 1.2921 0.0773 0.0000brechas_nivel <- t(sapply(levels(dm$nivel), function(k) {
p <- c(mujer = 1); inter <- paste0("nivel", k, ":mujer")
if (inter %in% names(coef(m3))) p[inter] <- 1
r <- comb(m3, V3, p); c(brecha_pct = pct(r[["estimacion"]]), r)
}))
round(brechas_nivel, 4)## brecha_pct estimacion error p_valor
## Primaria o menos -30.0275 -0.3571 0.0634 0.0000
## Secundaria -39.4257 -0.5013 0.0999 0.0000
## Media -27.9671 -0.3280 0.0504 0.0000
## Técnica o tecnológica -24.7271 -0.2841 0.0612 0.0000
## Universitaria 5.3626 0.0522 0.0616 0.3963
## Posgrado -4.8584 -0.0498 0.0711 0.4836# ¿Cambia el retorno al controlar por zona-estrato?
m1_sin <- lm(lingreso ~ educ + exper + exper2 + mujer, data = dm)
c(educ_sin_zona_estrato = coef(m1_sin)[["educ"]], educ_con_zona_estrato = coef(m1)[["educ"]])## educ_sin_zona_estrato educ_con_zona_estrato
## 0.12014004 0.09088018# Zona y estrato (modelo 1), prueba conjunta y diferencia por sexo
round(t(sapply(grep("zona_estrato", names(coef(m1)), value = TRUE), function(k) {
r <- comb(m1, V1, setNames(1, k)); c(pct = pct(r[["estimacion"]]), r)
})), 4)## pct estimacion error p_valor
## zona_estratoUrbano E2 14.7884 0.1379 0.0359 0.0001
## zona_estratoUrbano E3 28.4283 0.2502 0.0442 0.0000
## zona_estratoUrbano E4-6 115.6460 0.7685 0.0567 0.0000
## zona_estratoRural -13.4037 -0.1439 0.0467 0.0021waldtest(update(m1, . ~ . - zona_estrato), m1, vcov = V1, test = "F")## Wald test
##
## Model 1: lingreso ~ educ + exper + exper2 + mujer
## Model 2: lingreso ~ educ + exper + exper2 + mujer + zona_estrato
## Res.Df Df F Pr(>F)
## 1 2649
## 2 2645 4 55.078 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m1_int <- update(m1, . ~ . + zona_estrato:mujer)
waldtest(m1, m1_int, vcov = vcovHC(m1_int, type = "HC1"), test = "F")## Wald test
##
## Model 1: lingreso ~ educ + exper + exper2 + mujer + zona_estrato
## Model 2: lingreso ~ educ + exper + exper2 + mujer + zona_estrato + mujer:zona_estrato
## Res.Df Df F Pr(>F)
## 1 2645
## 2 2641 4 5.449 0.0002284 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Educación por zona-estrato (ponderado, descriptivo)
dm %>% group_by(zona_estrato) %>%
summarise(univ = 100 * weighted.mean(educ >= 16, w), educ_m = weighted.mean(educ, w))## # A tibble: 5 × 3
## zona_estrato univ educ_m
## <fct> <dbl> <dbl>
## 1 Urbano E1 10.4 9.90
## 2 Urbano E2 13.5 10.9
## 3 Urbano E3 40.6 13.3
## 4 Urbano E4-6 71.2 15.4
## 5 Rural 1.68 7.03# Perfil de experiencia (modelo 1): efecto marginal en la media y máximo
c(efecto_media = 100 * (coef(m1)[["exper"]] + 2 * coef(m1)[["exper2"]] / 100 * mean(dm$exper)),
maximo = -100 * coef(m1)[["exper"]] / (2 * coef(m1)[["exper2"]]))## efecto_media maximo
## 0.6559406 31.9153386Parte VI. Validación
17. Residuales
diagnostico <- function(m) {
e <- residuals(m); n <- length(e); s <- (e - mean(e)) / sd(e)
asim <- mean(s^3); curt <- mean(s^4); jb <- n / 6 * (asim^2 + (curt - 3)^2 / 4)
cook <- cooks.distance(m)
c(asimetria = asim, curtosis = curt, JB = jb, JB_p = pchisq(jb, 2, lower.tail = FALSE),
atipicos = sum(abs(rstudent(m)) > 3), influyentes = sum(cook > 4 / n), cook_max = max(cook))
}
tabla_resid <- t(sapply(modelos, diagnostico))
round(tabla_resid, 3)## asimetria curtosis JB JB_p atipicos influyentes
## (1) Años -0.623 4.910 574.995 0 22 140
## (2) Años x sexo -0.626 5.078 651.029 0 25 144
## (3) Niveles x sexo -0.608 5.001 606.604 0 29 139
## cook_max
## (1) Años 0.018
## (2) Años x sexo 0.029
## (3) Niveles x sexo 0.016m_d <- m3
d_res <- data.frame(ajustado = fitted(m_d), residuo = residuals(m_d), estud = rstudent(m_d),
cook = cooks.distance(m_d), nivel = dm$nivel, obs = seq_len(nrow(dm)))
r1 <- ggplot(d_res, aes(ajustado, residuo)) + geom_point(alpha = 0.25, color = "steelblue") +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_smooth(method = "loess", se = FALSE, color = "firebrick") +
labs(title = "Residuales frente a ajustados", x = "Valor ajustado", y = "Residual") + tema
r2 <- ggplot(d_res, aes(nivel, residuo)) + geom_jitter(alpha = 0.15, width = 0.2, color = "steelblue") +
geom_boxplot(fill = NA, outlier.shape = NA, width = 0.5) + geom_hline(yintercept = 0, linetype = "dashed") +
scale_x_discrete(labels = function(x) gsub(" o ", "\no ", x)) +
labs(title = "Residuales por nivel educativo", x = NULL, y = "Residual") + tema
r3 <- ggplot(d_res, aes(sample = estud)) + stat_qq(alpha = 0.25, color = "steelblue") +
stat_qq_line(color = "firebrick") + labs(title = "Gráfico Q-Q normal", x = "Teóricos", y = "Estudentizados") + tema
r4 <- ggplot(d_res, aes(obs, cook)) + geom_segment(aes(xend = obs, yend = 0), color = "steelblue", alpha = 0.5) +
geom_hline(yintercept = 4 / nrow(d_res), linetype = "dashed", color = "firebrick") +
labs(title = "Distancia de Cook", x = "Observación", y = "Cook") + tema
g_residuales <- function() {
grid::grid.newpage(); grid::pushViewport(grid::viewport(layout = grid::grid.layout(2, 2)))
pos <- list(c(1, 1), c(1, 2), c(2, 1), c(2, 2)); graf <- list(r1, r2, r3, r4)
for (i in 1:4) print(graf[[i]], vp = grid::viewport(layout.pos.row = pos[[i]][1], layout.pos.col = pos[[i]][2]))
}
g_residuales()
png(file.path(carpeta, "residuales.png"), width = 10, height = 8, units = "in", res = 300)
g_residuales(); invisible(dev.off())18. Forma funcional, significancia y ajuste
validar <- function(m, V) {
f_glo <- waldtest(lm(lingreso ~ 1, data = dm), m, vcov = V, test = "F")
reset <- resettest(m, power = 2:3, type = "fitted")
c(F_global = f_glo$F[2], F_p = f_glo$`Pr(>F)`[2], R2_ajust = summary(m)$adj.r.squared,
RESET = unname(reset$statistic), RESET_p = unname(reset$p.value))
}
tabla_valid <- t(mapply(validar, modelos, vcovs))
round(tabla_valid, 3)## F_global F_p R2_ajust RESET RESET_p
## (1) Años 204.022 0 0.370 29.349 0.000
## (2) Años x sexo 182.293 0 0.378 59.232 0.000
## (3) Niveles x sexo 139.040 0 0.410 1.313 0.269# Pruebas F entre modelos anidados (restringido primero, completo después)
anidados <- list(
"Educación x sexo (modelo 2 vs 1)" = waldtest(m1, m2, vcov = V2, test = "F"),
"Niveles x sexo (modelo 3 vs niveles)" = waldtest(m_niv, m3, vcov = V3, test = "F"),
"Zona-estrato en conjunto (modelo 1)" = waldtest(update(m1, . ~ . - zona_estrato), m1, vcov = V1, test = "F"))
round(t(sapply(anidados, function(x) c(gl = abs(x$Df[2]), F = x$F[2], p = x$`Pr(>F)`[2]))), 4)## gl F p
## Educación x sexo (modelo 2 vs 1) 1 33.7808 0
## Niveles x sexo (modelo 3 vs niveles) 5 8.6863 0
## Zona-estrato en conjunto (modelo 1) 4 55.0778 019. Multicolinealidad: VIF generalizado (GVIF)
gvif <- function(m) {
X <- model.matrix(m); asg <- attr(X, "assign")[-1]; X <- X[, -1, drop = FALSE]
R <- cor(X)
etq <- attr(terms(m), "term.labels")
res <- t(sapply(seq_along(etq), function(j) {
idx <- which(asg == j)
g <- det(R[idx, idx, drop = FALSE]) * det(R[-idx, -idx, drop = FALSE]) / det(R)
c(GVIF = g, gl = length(idx), GVIF_aj = g^(1 / (2 * length(idx))))
}))
rownames(res) <- etq; res
}
gvif_tabs <- lapply(modelos, gvif)
lapply(gvif_tabs, round, 2)## $`(1) Años`
## GVIF gl GVIF_aj
## educ 1.99 1 1.41
## exper 13.47 1 3.67
## exper2 14.26 1 3.78
## mujer 1.02 1 1.01
## zona_estrato 1.41 4 1.04
##
## $`(2) Años x sexo`
## GVIF gl GVIF_aj
## educ 2.96 1 1.72
## mujer 7.59 1 2.76
## exper 13.48 1 3.67
## exper2 14.28 1 3.78
## zona_estrato 1.43 4 1.05
## educ:mujer 8.98 1 3.00
##
## $`(3) Niveles x sexo`
## GVIF gl GVIF_aj
## nivel 52.09 5 1.48
## mujer 5.46 1 2.34
## exper 13.79 1 3.71
## exper2 14.43 1 3.80
## zona_estrato 1.58 4 1.06
## nivel:mujer 112.65 5 1.60gvif_max <- sapply(gvif_tabs, function(g) max(g[!grepl("exper|mujer", rownames(g)), "GVIF_aj"]))
gvif_max## (1) Años (2) Años x sexo (3) Niveles x sexo
## 1.410585 1.721401 1.48482120. Hipótesis específicas
prueba <- function(m, V, pesos_c, valor = 0) {
r <- comb(m, V, pesos_c); t <- (r[["estimacion"]] - valor) / r[["error"]]
c(estimacion = r[["estimacion"]], t = t, p_valor = 2 * pt(-abs(t), df.residual(m)))
}
m_horas <- update(m1, . ~ . + lhoras) # elasticidad del ingreso a las horas
m_chow <- lm(lingreso ~ (educ + exper + exper2 + zona_estrato) * mujer, data = dm)
m_comun <- lm(lingreso ~ educ + exper + exper2 + zona_estrato, data = dm)
chow <- waldtest(m_comun, m_chow, vcov = vcovHC(m_chow, type = "HC1"), test = "F")
tabla_hip <- data.frame(
hipotesis = c("Igual retorno a la educación por sexo",
"Brecha por sexo entre universitarios = 0",
"Prima de la media = prima de la secundaria",
"Elasticidad del ingreso a las horas = 1",
"Chow: igual ecuación por sexo"),
rbind(prueba(m2, V2, c("educ:mujer" = 1)),
prueba(m3, V3, c(mujer = 1, "nivelUniversitaria:mujer" = 1)),
prueba(m3, V3, c(nivelMedia = 1, nivelSecundaria = -1)),
prueba(m_horas, vcovHC(m_horas, type = "HC1"), c(lhoras = 1), valor = 1),
c(estimacion = NA, t = chow$F[2], p_valor = chow$`Pr(>F)`[2])))
tabla_hip %>% mutate(across(where(is.numeric), ~ round(.x, 4)))## hipotesis estimacion t p_valor
## 1 Igual retorno a la educación por sexo 0.0379 5.8121 0.0000
## 2 Brecha por sexo entre universitarios = 0 0.0522 0.8484 0.3963
## 3 Prima de la media = prima de la secundaria 0.0818 1.3609 0.1737
## 4 Elasticidad del ingreso a las horas = 1 0.7862 -5.3463 0.0000
## 5 Chow: igual ecuación por sexo NA 13.1249 0.0000Parte VII. Robustez
21. Especificaciones alternativas y pertinencia de no ponderar
lims <- quantile(dm$lingreso, c(0.01, 0.99))
dm$lingreso_w <- pmin(pmax(dm$lingreso, lims[1]), lims[2])
cook2 <- cooks.distance(m2)
alt <- list(
"Principal" = m2,
"Con log(horas)" = update(m2, . ~ . + lhoras),
"Ingreso por hora" = update(m2, lw ~ .),
"MCGF" = gls[[2]],
"Ponderado" = lm(formula(m2), data = dm, weights = w),
"Winsorizado 1-99" = update(m2, lingreso_w ~ ., data = dm),
"Sin influyentes" = update(m2, data = dm[cook2 <= 4 / nrow(dm), ]))
alt_V <- lapply(alt, vcovHC, type = "HC1")
# Errores agrupados por vivienda (sensibilidad a la dependencia dentro de la vivienda)
alt[["Agrupado vivienda"]] <- m2
alt_V[["Agrupado vivienda"]] <- vcovCL(m2, cluster = ~ vivienda, type = "HC1")
cantidades_rob <- list(
"Retorno hombres (%)" = c(educ = 1),
"Retorno mujeres (%)" = c(educ = 1, "educ:mujer" = 1),
"Brecha con 11 años (%)" = c(mujer = 1, "educ:mujer" = 11),
"Estratos 4-6 vs E1 (%)" = c("zona_estratoUrbano E4-6" = 1),
"Rural vs urbano E1 (%)" = c(zona_estratoRural = 1))
celda <- function(m, V, p) { r <- comb(m, V, p); paste0(fx(pct(r[["estimacion"]]), 1), "$^{", estrellas(r[["p_valor"]]), "}$") }
tabla_rob <- sapply(names(alt), function(k) sapply(cantidades_rob, function(p) celda(alt[[k]], alt_V[[k]], p)))
tabla_rob## Principal Con log(horas)
## Retorno hombres (%) "7{,}5$^{***}$" "8{,}0$^{***}$"
## Retorno mujeres (%) "11{,}7$^{***}$" "10{,}6$^{***}$"
## Brecha con 11 años (%) "$-$20{,}5$^{***}$" "$-$10{,}9$^{***}$"
## Estratos 4-6 vs E1 (%) "113{,}8$^{***}$" "101{,}6$^{***}$"
## Rural vs urbano E1 (%) "$-$15{,}6$^{***}$" "$-$6{,}4$^{}$"
## Ingreso por hora MCGF
## Retorno hombres (%) "8{,}1$^{***}$" "7{,}2$^{***}$"
## Retorno mujeres (%) "10{,}2$^{***}$" "12{,}7$^{***}$"
## Brecha con 11 años (%) "$-$7{,}9$^{***}$" "$-$21{,}0$^{***}$"
## Estratos 4-6 vs E1 (%) "98{,}1$^{***}$" "110{,}0$^{***}$"
## Rural vs urbano E1 (%) "$-$3{,}5$^{}$" "$-$15{,}0$^{***}$"
## Ponderado Winsorizado 1-99
## Retorno hombres (%) "5{,}6$^{***}$" "7{,}5$^{***}$"
## Retorno mujeres (%) "8{,}6$^{***}$" "11{,}6$^{***}$"
## Brecha con 11 años (%) "$-$27{,}2$^{***}$" "$-$20{,}0$^{***}$"
## Estratos 4-6 vs E1 (%) "143{,}1$^{***}$" "108{,}7$^{***}$"
## Rural vs urbano E1 (%) "$-$24{,}2$^{***}$" "$-$15{,}6$^{***}$"
## Sin influyentes Agrupado vivienda
## Retorno hombres (%) "7{,}7$^{***}$" "7{,}5$^{***}$"
## Retorno mujeres (%) "13{,}5$^{***}$" "11{,}7$^{***}$"
## Brecha con 11 años (%) "$-$21{,}6$^{***}$" "$-$20{,}5$^{***}$"
## Estratos 4-6 vs E1 (%) "111{,}8$^{***}$" "113{,}8$^{***}$"
## Rural vs urbano E1 (%) "$-$13{,}2$^{***}$" "$-$15{,}6$^{***}$"# Prueba de DuMouchel y Duncan (1983): ¿los pesos cambian los coeficientes?
# H0: las estimaciones ponderada y sin ponderar coinciden. Si no se rechaza, no ponderar es adecuado.
dumouchel <- function(m) {
X <- model.matrix(m); Xw <- X * dm$w_norm
colnames(Xw) <- paste0("w_", make.names(colnames(X)))
aum <- lm(dm$lingreso ~ X + Xw - 1)
base <- lm(dm$lingreso ~ X - 1)
if (any(is.na(coef(aum)))) return(c(F = NA, gl = NA, p_valor = NA)) # columnas colineales
r <- waldtest(base, aum, vcov = vcovHC(aum, type = "HC1"), test = "F")
c(F = r$F[2], gl = abs(r$Df[2]), p_valor = r$`Pr(>F)`[2])
}
tabla_dd <- t(sapply(modelos, dumouchel))
round(tabla_dd, 4)## F gl p_valor
## (1) Años 8.0126 9 0
## (2) Años x sexo 6.2537 10 0
## (3) Niveles x sexo 3.3135 18 0Parte VIII. Tablas para LaTeX
guardar <- function(lineas, archivo) {
writeLines(lineas, file.path(carpeta, archivo), useBytes = TRUE); cat(lineas, sep = "\n"); cat("\n")
}
filas_coef <- function(mods, vcs, vars, etiquetas) {
out <- c()
for (i in seq_along(vars)) {
est <- se <- rep("", length(mods))
for (j in seq_along(mods)) {
b <- coef(mods[[j]])
if (vars[i] %in% names(b)) {
s <- sqrt(diag(vcs[[j]]))[[vars[i]]]
p <- 2 * pt(-abs(b[[vars[i]]] / s), df.residual(mods[[j]]))
est[j] <- paste0(fx(b[[vars[i]]]), "$^{", estrellas(p), "}$"); se[j] <- paste0("(", fx(s), ")")
}
}
out <- c(out, fila(etiquetas[i], est), fila("", se))
}
out
}
# ---------- Tabla de resultados ----------
vars_r <- c("educ", "educ:mujer", paste0("nivel", niveles_k), "mujer", paste0("nivel", niveles_k, ":mujer"),
"exper", "exper2", "zona_estratoUrbano E2", "zona_estratoUrbano E3",
"zona_estratoUrbano E4-6", "zona_estratoRural", "(Intercept)")
etq_r <- c("Años de educación", "Educación $\\times$ mujer", niveles_k, "Mujer",
paste0(niveles_k, " $\\times$ mujer"), "Experiencia", "Experiencia$^2$/100",
"Estrato 2 urbano", "Estrato 3 urbano", "Estratos 4 a 6 urbanos", "Rural", "Constante")
guardar(c("\\begin{table}[H]", "\\centering",
"\\caption{Determinantes del ingreso laboral mensual. MCO}", "\\label{tab:resultados}", "\\small",
"\\begin{tabular}{lccc}", "\\toprule",
fila("", c("(1) Años", "(2) Años $\\times$ sexo", "(3) Niveles $\\times$ sexo")), "\\midrule",
filas_coef(modelos, vcovs, vars_r, etq_r), "\\midrule",
fila("Observaciones", fi(sapply(modelos, nobs))),
fila("$R^2$", fx(sapply(modelos, function(m) summary(m)$r.squared))),
"\\bottomrule", "\\end{tabular}", "\\par\\vspace{2pt}",
"\\begin{minipage}{0.95\\textwidth}",
paste("\\footnotesize Variable dependiente: logaritmo del ingreso laboral mensual. Estimaciones por MCO sin",
"ponderar. Errores estándar robustos (HC1) entre paréntesis.",
"*** $p<0{,}01$, ** $p<0{,}05$, * $p<0{,}1$. Categorías base: hombre, primaria o menos y estrato 1",
"urbano. En (2), \\textit{Mujer} es la brecha con cero años de educación; en (3), en primaria o menos.",
"Fuente: GEIH 2025-IV, DANE."),
"\\end{minipage}", "\\end{table}"), "tabla_resultados.tex")## \begin{table}[H]
## \centering
## \caption{Determinantes del ingreso laboral mensual. MCO}
## \label{tab:resultados}
## \small
## \begin{tabular}{lccc}
## \toprule
## & (1) Años & (2) Años $\times$ sexo & (3) Niveles $\times$ sexo \\
## \midrule
## Años de educación & 0{,}091$^{***}$ & 0{,}073$^{***}$ & \\
## & (0{,}005) & (0{,}005) & \\
## Educación $\times$ mujer & & 0{,}038$^{***}$ & \\
## & & (0{,}007) & \\
## Secundaria & & & 0{,}149$^{**}$ \\
## & & & (0{,}065) \\
## Media & & & 0{,}230$^{***}$ \\
## & & & (0{,}054) \\
## Técnica o tecnológica & & & 0{,}458$^{***}$ \\
## & & & (0{,}064) \\
## Universitaria & & & 0{,}666$^{***}$ \\
## & & & (0{,}066) \\
## Posgrado & & & 1{,}292$^{***}$ \\
## & & & (0{,}077) \\
## Mujer & $-$0{,}230$^{***}$ & $-$0{,}646$^{***}$ & $-$0{,}357$^{***}$ \\
## & (0{,}027) & (0{,}080) & (0{,}063) \\
## Secundaria $\times$ mujer & & & $-$0{,}144$^{}$ \\
## & & & (0{,}118) \\
## Media $\times$ mujer & & & 0{,}029$^{}$ \\
## & & & (0{,}081) \\
## Técnica o tecnológica $\times$ mujer & & & 0{,}073$^{}$ \\
## & & & (0{,}088) \\
## Universitaria $\times$ mujer & & & 0{,}409$^{***}$ \\
## & & & (0{,}088) \\
## Posgrado $\times$ mujer & & & 0{,}307$^{***}$ \\
## & & & (0{,}095) \\
## Experiencia & 0{,}025$^{***}$ & 0{,}026$^{***}$ & 0{,}034$^{***}$ \\
## & (0{,}004) & (0{,}004) & (0{,}004) \\
## Experiencia$^2$/100 & $-$0{,}040$^{***}$ & $-$0{,}041$^{***}$ & $-$0{,}062$^{***}$ \\
## & (0{,}007) & (0{,}007) & (0{,}007) \\
## Estrato 2 urbano & 0{,}138$^{***}$ & 0{,}135$^{***}$ & 0{,}146$^{***}$ \\
## & (0{,}036) & (0{,}036) & (0{,}035) \\
## Estrato 3 urbano & 0{,}250$^{***}$ & 0{,}250$^{***}$ & 0{,}209$^{***}$ \\
## & (0{,}044) & (0{,}044) & (0{,}042) \\
## Estratos 4 a 6 urbanos & 0{,}768$^{***}$ & 0{,}760$^{***}$ & 0{,}587$^{***}$ \\
## & (0{,}057) & (0{,}056) & (0{,}056) \\
## Rural & $-$0{,}144$^{***}$ & $-$0{,}170$^{***}$ & $-$0{,}230$^{***}$ \\
## & (0{,}047) & (0{,}047) & (0{,}046) \\
## Constante & 12{,}651$^{***}$ & 12{,}846$^{***}$ & 13{,}275$^{***}$ \\
## & (0{,}072) & (0{,}076) & (0{,}067) \\
## \midrule
## Observaciones & 2.654 & 2.654 & 2.654 \\
## $R^2$ & 0{,}371 & 0{,}380 & 0{,}414 \\
## \bottomrule
## \end{tabular}
## \par\vspace{2pt}
## \begin{minipage}{0.95\textwidth}
## \footnotesize Variable dependiente: logaritmo del ingreso laboral mensual. Estimaciones por MCO sin ponderar. Errores estándar robustos (HC1) entre paréntesis. *** $p<0{,}01$, ** $p<0{,}05$, * $p<0{,}1$. Categorías base: hombre, primaria o menos y estrato 1 urbano. En (2), \textit{Mujer} es la brecha con cero años de educación; en (3), en primaria o menos. Fuente: GEIH 2025-IV, DANE.
## \end{minipage}
## \end{table}# ---------- Tabla de corrección de la heterocedasticidad (modelo 2) ----------
vars_c <- vars_r[vars_r %in% names(coef(m2)) & vars_r != "(Intercept)"]
etq_c <- etq_r[match(vars_c, vars_r)]
guardar(c("\\begin{table}[H]", "\\centering",
"\\caption{Corrección de la heterocedasticidad en el modelo (2)}", "\\label{tab:correccion}", "\\small",
"\\begin{tabular}{lccc}", "\\toprule",
fila("", c("MCO", "MCO", "MCGF")), fila("", c("(errores clásicos)", "(errores robustos)", "(errores robustos)")),
"\\midrule",
filas_coef(list(m2, m2, gls[[2]]), list(vcov(m2), vcovs[[2]], vcovs_g[[2]]), vars_c, etq_c),
"\\midrule",
fila("Breusch-Pagan", c(fx(pruebas_het[2, "BP"], 1), fx(pruebas_het[2, "BP"], 1), fx(bp_mcgf[2, "BP"], 1))),
fila("\\quad $p$-valor", c(fp(pruebas_het[2, "BP_p"]), fp(pruebas_het[2, "BP_p"]), fp(bp_mcgf[2, "p_valor"]))),
"\\bottomrule", "\\end{tabular}", "\\par\\vspace{2pt}",
"\\begin{minipage}{0.92\\textwidth}",
paste("\\footnotesize Errores estándar entre paréntesis. Estimaciones sin ponderar. MCGF: mínimos cuadrados",
"ponderados con peso igual al inverso de la varianza estimada del error (función exponencial de los",
"regresores). En MCGF, Breusch-Pagan sobre los residuos del modelo transformado (aproximación).",
"Fuente: GEIH 2025-IV, DANE."),
"\\end{minipage}", "\\end{table}"), "tabla_correccion.tex")## \begin{table}[H]
## \centering
## \caption{Corrección de la heterocedasticidad en el modelo (2)}
## \label{tab:correccion}
## \small
## \begin{tabular}{lccc}
## \toprule
## & MCO & MCO & MCGF \\
## & (errores clásicos) & (errores robustos) & (errores robustos) \\
## \midrule
## Años de educación & 0{,}073$^{***}$ & 0{,}073$^{***}$ & 0{,}069$^{***}$ \\
## & (0{,}005) & (0{,}005) & (0{,}005) \\
## Educación $\times$ mujer & 0{,}038$^{***}$ & 0{,}038$^{***}$ & 0{,}050$^{***}$ \\
## & (0{,}006) & (0{,}007) & (0{,}006) \\
## Mujer & $-$0{,}646$^{***}$ & $-$0{,}646$^{***}$ & $-$0{,}786$^{***}$ \\
## & (0{,}073) & (0{,}080) & (0{,}079) \\
## Experiencia & 0{,}026$^{***}$ & 0{,}026$^{***}$ & 0{,}027$^{***}$ \\
## & (0{,}003) & (0{,}004) & (0{,}004) \\
## Experiencia$^2$/100 & $-$0{,}041$^{***}$ & $-$0{,}041$^{***}$ & $-$0{,}046$^{***}$ \\
## & (0{,}007) & (0{,}007) & (0{,}007) \\
## Estrato 2 urbano & 0{,}135$^{***}$ & 0{,}135$^{***}$ & 0{,}129$^{***}$ \\
## & (0{,}035) & (0{,}036) & (0{,}034) \\
## Estrato 3 urbano & 0{,}250$^{***}$ & 0{,}250$^{***}$ & 0{,}228$^{***}$ \\
## & (0{,}043) & (0{,}044) & (0{,}042) \\
## Estratos 4 a 6 urbanos & 0{,}760$^{***}$ & 0{,}760$^{***}$ & 0{,}742$^{***}$ \\
## & (0{,}059) & (0{,}056) & (0{,}055) \\
## Rural & $-$0{,}170$^{***}$ & $-$0{,}170$^{***}$ & $-$0{,}162$^{***}$ \\
## & (0{,}044) & (0{,}047) & (0{,}046) \\
## \midrule
## Breusch-Pagan & 65{,}7 & 65{,}7 & 16{,}3 \\
## \quad $p$-valor & $<${}0{,}001 & $<${}0{,}001 & 0{,}061 \\
## \bottomrule
## \end{tabular}
## \par\vspace{2pt}
## \begin{minipage}{0.92\textwidth}
## \footnotesize Errores estándar entre paréntesis. Estimaciones sin ponderar. MCGF: mínimos cuadrados ponderados con peso igual al inverso de la varianza estimada del error (función exponencial de los regresores). En MCGF, Breusch-Pagan sobre los residuos del modelo transformado (aproximación). Fuente: GEIH 2025-IV, DANE.
## \end{minipage}
## \end{table}# ---------- Tabla de diagnóstico ----------
nm <- names(modelos)
guardar(c("\\begin{table}[H]", "\\centering",
"\\caption{Pruebas de especificación y diagnóstico}", "\\label{tab:diagnostico}", "\\small",
"\\begin{tabular}{lccc}", "\\toprule", fila("", c("(1)", "(2)", "(3)")), "\\midrule",
grupo("Heterocedasticidad", 4),
fila("Breusch-Pagan", fx(pruebas_het[nm, "BP"], 1)), fila("\\quad $p$-valor", fp(pruebas_het[nm, "BP_p"])),
fila("White", fx(pruebas_het[nm, "White"], 1)), fila("\\quad $p$-valor", fp(pruebas_het[nm, "White_p"])),
grupo("Normalidad", 4),
fila("Asimetría", fx(tabla_resid[nm, "asimetria"], 2)), fila("Curtosis", fx(tabla_resid[nm, "curtosis"], 2)),
fila("Jarque-Bera", fx(tabla_resid[nm, "JB"], 1)), fila("\\quad $p$-valor", fp(tabla_resid[nm, "JB_p"])),
grupo("Forma funcional y ajuste", 4),
fila("RESET", fx(tabla_valid[nm, "RESET"], 2)), fila("\\quad $p$-valor", fp(tabla_valid[nm, "RESET_p"])),
fila("F global", fx(tabla_valid[nm, "F_global"], 1)), fila("$R^2$ ajustado", fx(tabla_valid[nm, "R2_ajust"])),
grupo("Pertinencia de no ponderar", 4),
fila("DuMouchel-Duncan ($F$)", ifelse(is.na(tabla_dd[nm, "F"]), "---", fx(tabla_dd[nm, "F"], 2))),
fila("\\quad $p$-valor", ifelse(is.na(tabla_dd[nm, "p_valor"]), "---", fp(tabla_dd[nm, "p_valor"]))),
grupo("Observaciones influyentes y multicolinealidad", 4),
fila("Atípicos ($|r_i| > 3$)", tabla_resid[nm, "atipicos"]),
fila("Distancia de Cook máxima", fx(tabla_resid[nm, "cook_max"])),
fila("GVIF$^{1/(2gl)}$ máximo", fx(gvif_max, 2)),
"\\bottomrule", "\\end{tabular}", "\\par\\vspace{2pt}",
"\\begin{minipage}{0.92\\textwidth}",
paste("\\footnotesize Modelos estimados por MCO sin ponderar. Breusch-Pagan (Koenker) y White (forma especial):",
"$H_0$ homocedasticidad. Jarque-Bera: $H_0$ normalidad. RESET con cuadrados y cubos de los ajustados:",
"$H_0$ forma funcional correcta. F global con errores robustos. DuMouchel-Duncan: $H_0$ las estimaciones",
"ponderada y sin ponderar coinciden. El GVIF excluye la experiencia, la indicadora de mujer y sus",
"interacciones (colinealidad estructural). Fuente: GEIH 2025-IV, DANE."),
"\\end{minipage}", "\\end{table}"), "tabla_diagnostico.tex")## \begin{table}[H]
## \centering
## \caption{Pruebas de especificación y diagnóstico}
## \label{tab:diagnostico}
## \small
## \begin{tabular}{lccc}
## \toprule
## & (1) & (2) & (3) \\
## \midrule
## \multicolumn{4}{l}{\textit{Heterocedasticidad}} \\
## Breusch-Pagan & 62{,}2 & 65{,}7 & 89{,}6 \\
## \quad $p$-valor & $<${}0{,}001 & $<${}0{,}001 & $<${}0{,}001 \\
## White & 15{,}5 & 19{,}7 & 31{,}1 \\
## \quad $p$-valor & $<${}0{,}001 & $<${}0{,}001 & $<${}0{,}001 \\
## \multicolumn{4}{l}{\textit{Normalidad}} \\
## Asimetría & $-$0{,}62 & $-$0{,}63 & $-$0{,}61 \\
## Curtosis & 4{,}91 & 5{,}08 & 5{,}00 \\
## Jarque-Bera & 575{,}0 & 651{,}0 & 606{,}6 \\
## \quad $p$-valor & $<${}0{,}001 & $<${}0{,}001 & $<${}0{,}001 \\
## \multicolumn{4}{l}{\textit{Forma funcional y ajuste}} \\
## RESET & 29{,}35 & 59{,}23 & 1{,}31 \\
## \quad $p$-valor & $<${}0{,}001 & $<${}0{,}001 & 0{,}269 \\
## F global & 204{,}0 & 182{,}3 & 139{,}0 \\
## $R^2$ ajustado & 0{,}370 & 0{,}378 & 0{,}410 \\
## \multicolumn{4}{l}{\textit{Pertinencia de no ponderar}} \\
## DuMouchel-Duncan ($F$) & 8{,}01 & 6{,}25 & 3{,}31 \\
## \quad $p$-valor & $<${}0{,}001 & $<${}0{,}001 & $<${}0{,}001 \\
## \multicolumn{4}{l}{\textit{Observaciones influyentes y multicolinealidad}} \\
## Atípicos ($|r_i| > 3$) & 22 & 25 & 29 \\
## Distancia de Cook máxima & 0{,}018 & 0{,}029 & 0{,}016 \\
## GVIF$^{1/(2gl)}$ máximo & 1{,}41 & 1{,}72 & 1{,}48 \\
## \bottomrule
## \end{tabular}
## \par\vspace{2pt}
## \begin{minipage}{0.92\textwidth}
## \footnotesize Modelos estimados por MCO sin ponderar. Breusch-Pagan (Koenker) y White (forma especial): $H_0$ homocedasticidad. Jarque-Bera: $H_0$ normalidad. RESET con cuadrados y cubos de los ajustados: $H_0$ forma funcional correcta. F global con errores robustos. DuMouchel-Duncan: $H_0$ las estimaciones ponderada y sin ponderar coinciden. El GVIF excluye la experiencia, la indicadora de mujer y sus interacciones (colinealidad estructural). Fuente: GEIH 2025-IV, DANE.
## \end{minipage}
## \end{table}# ---------- Tabla de hipótesis ----------
guardar(c("\\begin{table}[H]", "\\centering", "\\caption{Contrastes de hipótesis}", "\\label{tab:hipotesis}",
"\\small", "\\begin{tabular}{lccc}", "\\toprule",
fila("Hipótesis nula", c("Estimación", "Estadístico", "$p$-valor")), "\\midrule",
sapply(seq_len(nrow(tabla_hip)), function(i)
fila(tabla_hip$hipotesis[i], c(ifelse(is.na(tabla_hip$estimacion[i]), "---", fx(tabla_hip$estimacion[i])),
fx(tabla_hip$t[i], 2), fp(tabla_hip$p_valor[i])))),
"\\bottomrule", "\\end{tabular}", "\\par\\vspace{2pt}",
"\\footnotesize Estadístico $t$ con errores robustos, salvo Chow ($F$). La elasticidad a las horas se estima agregando el log de las horas al modelo (1).",
"\\end{table}"), "tabla_hipotesis.tex")## \begin{table}[H]
## \centering
## \caption{Contrastes de hipótesis}
## \label{tab:hipotesis}
## \small
## \begin{tabular}{lccc}
## \toprule
## Hipótesis nula & Estimación & Estadístico & $p$-valor \\
## \midrule
## Igual retorno a la educación por sexo & 0{,}038 & 5{,}81 & $<${}0{,}001 \\
## Brecha por sexo entre universitarios = 0 & 0{,}052 & 0{,}85 & 0{,}396 \\
## Prima de la media = prima de la secundaria & 0{,}082 & 1{,}36 & 0{,}174 \\
## Elasticidad del ingreso a las horas = 1 & 0{,}786 & $-$5{,}35 & $<${}0{,}001 \\
## Chow: igual ecuación por sexo & --- & 13{,}12 & $<${}0{,}001 \\
## \bottomrule
## \end{tabular}
## \par\vspace{2pt}
## \footnotesize Estadístico $t$ con errores robustos, salvo Chow ($F$). La elasticidad a las horas se estima agregando el log de las horas al modelo (1).
## \end{table}# ---------- Tabla de robustez ----------
guardar(c("\\begin{table}[H]", "\\centering", "\\caption{Robustez de los resultados del modelo (2)}",
"\\label{tab:robustez}", "\\footnotesize",
paste0("\\begin{tabular}{l", strrep("c", ncol(tabla_rob)), "}"), "\\toprule",
fila("", colnames(tabla_rob)), "\\midrule",
sapply(seq_len(nrow(tabla_rob)), function(i) fila(rownames(tabla_rob)[i], tabla_rob[i, ])),
"\\bottomrule", "\\end{tabular}", "\\par\\vspace{2pt}",
"\\footnotesize Efectos en porcentaje. *** $p<0{,}01$, ** $p<0{,}05$, * $p<0{,}1$, con errores robustos (HC1), salvo la columna agrupada por vivienda. La columna ponderada usa el factor de expansión.",
"\\end{table}"), "tabla_robustez.tex")## \begin{table}[H]
## \centering
## \caption{Robustez de los resultados del modelo (2)}
## \label{tab:robustez}
## \footnotesize
## \begin{tabular}{lcccccccc}
## \toprule
## & Principal & Con log(horas) & Ingreso por hora & MCGF & Ponderado & Winsorizado 1-99 & Sin influyentes & Agrupado vivienda \\
## \midrule
## Retorno hombres (%) & 7{,}5$^{***}$ & 8{,}0$^{***}$ & 8{,}1$^{***}$ & 7{,}2$^{***}$ & 5{,}6$^{***}$ & 7{,}5$^{***}$ & 7{,}7$^{***}$ & 7{,}5$^{***}$ \\
## Retorno mujeres (%) & 11{,}7$^{***}$ & 10{,}6$^{***}$ & 10{,}2$^{***}$ & 12{,}7$^{***}$ & 8{,}6$^{***}$ & 11{,}6$^{***}$ & 13{,}5$^{***}$ & 11{,}7$^{***}$ \\
## Brecha con 11 años (%) & $-$20{,}5$^{***}$ & $-$10{,}9$^{***}$ & $-$7{,}9$^{***}$ & $-$21{,}0$^{***}$ & $-$27{,}2$^{***}$ & $-$20{,}0$^{***}$ & $-$21{,}6$^{***}$ & $-$20{,}5$^{***}$ \\
## Estratos 4-6 vs E1 (%) & 113{,}8$^{***}$ & 101{,}6$^{***}$ & 98{,}1$^{***}$ & 110{,}0$^{***}$ & 143{,}1$^{***}$ & 108{,}7$^{***}$ & 111{,}8$^{***}$ & 113{,}8$^{***}$ \\
## Rural vs urbano E1 (%) & $-$15{,}6$^{***}$ & $-$6{,}4$^{}$ & $-$3{,}5$^{}$ & $-$15{,}0$^{***}$ & $-$24{,}2$^{***}$ & $-$15{,}6$^{***}$ & $-$13{,}2$^{***}$ & $-$15{,}6$^{***}$ \\
## \bottomrule
## \end{tabular}
## \par\vspace{2pt}
## \footnotesize Efectos en porcentaje. *** $p<0{,}01$, ** $p<0{,}05$, * $p<0{,}1$, con errores robustos (HC1), salvo la columna agrupada por vivienda. La columna ponderada usa el factor de expansión.
## \end{table}# ---------- Tabla de flujo ----------
guardar(c("\\begin{table}[H]", "\\centering", "\\caption{Flujo de depuración de la muestra}", "\\label{tab:flujo}",
"\\small", "\\begin{tabular}{lrr}", "\\toprule", fila("Paso", c("Observaciones", "Excluidas")), "\\midrule",
sapply(seq_len(nrow(tabla_flujo)), function(i)
fila(tabla_flujo$paso[i], c(fi(tabla_flujo$observaciones[i]),
ifelse(is.na(tabla_flujo$excluidas[i]), "---", fi(tabla_flujo$excluidas[i]))))),
"\\bottomrule", "\\end{tabular}", "\\par\\vspace{2pt}", "\\footnotesize Fuente: GEIH 2025-IV, DANE.",
"\\end{table}"), "tabla_flujo.tex")## \begin{table}[H]
## \centering
## \caption{Flujo de depuración de la muestra}
## \label{tab:flujo}
## \small
## \begin{tabular}{lrr}
## \toprule
## Paso & Observaciones & Excluidas \\
## \midrule
## Registros del trimestre & 6.847 & --- \\
## Ocupados con ingreso laboral positivo & 2.875 & --- \\
## Con educación, experiencia y nivel educativo & 2.875 & 0 \\
## Edad de 15 a 65 años y experiencia no negativa & 2.682 & 193 \\
## Al menos 5 horas de trabajo semanal & 2.668 & 14 \\
## Con dato de estrato (urbanos) & 2.654 & 14 \\
## \bottomrule
## \end{tabular}
## \par\vspace{2pt}
## \footnotesize Fuente: GEIH 2025-IV, DANE.
## \end{table}# ---------- Tabla por sexo ----------
proporciones <- c("Zona rural (proporción)", "No cotiza a pensión (proporción)")
filas_sexo <- sapply(seq_len(nrow(tabla_sexo)), function(i) {
d <- if (i <= 2) 2 else if (rownames(tabla_sexo)[i] %in% proporciones) 3 else 1
fila(rownames(tabla_sexo)[i], c(fx(tabla_sexo[i, "hombres"], d), fx(tabla_sexo[i, "mujeres"], d),
fx(tabla_sexo[i, "diferencia"], d), fp(tabla_sexo[i, "p_valor"])))
})
guardar(c("\\begin{table}[H]", "\\centering", "\\caption{Características de los ocupados por sexo}",
"\\label{tab:sexo}", "\\small", "\\begin{tabular}{lcccc}", "\\toprule",
fila("", c("Hombres", "Mujeres", "Diferencia", "$p$-valor")), "\\midrule", filas_sexo,
fila("Observaciones", c(fi(sum(df$mujer == 0)), fi(sum(df$mujer == 1)), "", "")),
"\\bottomrule", "\\end{tabular}", "\\par\\vspace{2pt}",
"\\footnotesize Medias ponderadas con el factor de expansión. Diferencia mujeres menos hombres, con errores robustos. Fuente: GEIH 2025-IV, DANE.",
"\\end{table}"), "tabla_sexo.tex")## \begin{table}[H]
## \centering
## \caption{Características de los ocupados por sexo}
## \label{tab:sexo}
## \small
## \begin{tabular}{lcccc}
## \toprule
## & Hombres & Mujeres & Diferencia & $p$-valor \\
## \midrule
## Log del ingreso mensual & 13{,}69 & 13{,}50 & $-$0{,}19 & 0{,}002 \\
## Log del ingreso por hora & 8{,}62 & 8{,}62 & 0{,}00 & 0{,}954 \\
## Horas mensuales & 171{,}1 & 148{,}3 & $-$22{,}8 & $<${}0{,}001 \\
## Años de educación & 8{,}0 & 9{,}7 & 1{,}7 & $<${}0{,}001 \\
## Experiencia potencial & 24{,}9 & 23{,}0 & $-$1{,}9 & 0{,}079 \\
## Zona rural (proporción) & 0{,}644 & 0{,}485 & $-$0{,}159 & $<${}0{,}001 \\
## No cotiza a pensión (proporción) & 0{,}825 & 0{,}772 & $-$0{,}054 & 0{,}022 \\
## Observaciones & 1.419 & 1.235 & & \\
## \bottomrule
## \end{tabular}
## \par\vspace{2pt}
## \footnotesize Medias ponderadas con el factor de expansión. Diferencia mujeres menos hombres, con errores robustos. Fuente: GEIH 2025-IV, DANE.
## \end{table}# ---------- Tabla por posición ocupacional ----------
tp <- tab_pos %>% filter(grupo != "Otro") %>% arrange(match(grupo, c("Asalariados", "Cuenta propia", "Patrones")))
guardar(c("\\begin{table}[H]", "\\centering", "\\caption{Ocupados según posición ocupacional}",
"\\label{tab:posicion}", "\\small", "\\begin{tabular}{lccc}", "\\toprule", fila("", tp$grupo), "\\midrule",
fila("Participación (\\%)", fx(tp$pct_ocup, 1)), fila("No cotiza a pensión (\\%)", fx(tp$no_cotiza, 1)),
fila("Ingreso mensual mediano (\\$)", fi(tp$mediana)), fila("Gana menos del mínimo (\\%)", fx(tp$bajo_min, 1)),
fila("Horas mensuales", fx(tp$horas_mes, 1)), fila("Años de educación", fx(tp$educ, 1)),
fila("Mujeres (\\%)", fx(tp$mujeres, 1)), fila("Observaciones", fi(tp$n)),
"\\bottomrule", "\\end{tabular}", "\\par\\vspace{2pt}",
"\\footnotesize Estadísticas ponderadas. Asalariados: empleados particulares, del gobierno, domésticos y jornaleros. Fuente: GEIH 2025-IV, DANE.",
"\\end{table}"), "tabla_posicion.tex")## \begin{table}[H]
## \centering
## \caption{Ocupados según posición ocupacional}
## \label{tab:posicion}
## \small
## \begin{tabular}{lccc}
## \toprule
## & Asalariados & Cuenta propia & Patrones \\
## \midrule
## Participación (\%) & 30{,}0 & 65{,}4 & 1{,}9 \\
## No cotiza a pensión (\%) & 50{,}0 & 93{,}6 & 92{,}0 \\
## Ingreso mensual mediano (\$) & 1.423.500 & 600.000 & 2.000.000 \\
## Gana menos del mínimo (\%) & 42{,}4 & 86{,}3 & 25{,}4 \\
## Horas mensuales & 180{,}3 & 152{,}2 & 210{,}0 \\
## Años de educación & 11{,}0 & 7{,}6 & 8{,}7 \\
## Mujeres (\%) & 48{,}1 & 31{,}9 & 20{,}2 \\
## Observaciones & 1.229 & 1.356 & 51 \\
## \bottomrule
## \end{tabular}
## \par\vspace{2pt}
## \footnotesize Estadísticas ponderadas. Asalariados: empleados particulares, del gobierno, domésticos y jornaleros. Fuente: GEIH 2025-IV, DANE.
## \end{table}Anexo técnico. Justificación de la especificación por niveles
j_lineal <- lm(lingreso ~ educ + exper + exper2 + mujer + zona_estrato, data = dm)
j_cuad <- update(j_lineal, . ~ . + I(educ^2))
j_niv <- lm(lingreso ~ nivel + exper + exper2 + mujer + zona_estrato, data = dm)
j_ambos <- lm(lingreso ~ educ + nivel + exper + exper2 + mujer + zona_estrato, data = dm)
# 1. Residuo medio del modelo lineal por años de educación (grupos con 30 o más)
data.frame(educ = dm$educ, e = residuals(j_lineal)) %>%
group_by(educ) %>% summarise(n = n(), residuo_medio = round(mean(e), 3)) %>% filter(n >= 30)## # A tibble: 18 × 3
## educ n residuo_medio
## <dbl> <int> <dbl>
## 1 0 36 0.519
## 2 2 50 0.068
## 3 3 60 0.282
## 4 4 54 0.193
## 5 5 294 0.083
## 6 6 60 0.074
## 7 7 56 -0.067
## 8 8 56 -0.095
## 9 9 72 -0.076
## 10 10 37 -0.158
## 11 11 784 -0.15
## 12 12 52 -0.229
## 13 13 287 -0.115
## 14 14 166 -0.095
## 15 15 41 -0.368
## 16 16 356 0.176
## 17 17 117 0.376
## 18 18 48 0.51# 2. Término cuadrático (convexidad solo si el coeficiente es positivo y significativo)
coeftest(j_cuad, vcov. = vcovHC(j_cuad, type = "HC1"))[c("educ", "I(educ^2)"), ]## Estimate Std. Error t value Pr(>|t|)
## educ -0.093249867 0.0154059590 -6.052844 1.624826e-09
## I(educ^2) 0.009206749 0.0007015655 13.123150 3.682131e-38# 3. ¿Los niveles aportan información a los años? (consistente con efectos diploma; no causal)
waldtest(update(j_ambos, . ~ . - nivel), j_ambos, vcov = vcovHC(j_ambos, type = "HC1"), test = "F")## Wald test
##
## Model 1: lingreso ~ educ + exper + exper2 + mujer + zona_estrato
## Model 2: lingreso ~ educ + nivel + exper + exper2 + mujer + zona_estrato
## Res.Df Df F Pr(>F)
## 1 2645
## 2 2640 5 46.745 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 4. RESET por especificación
sapply(list(Lineal = j_lineal, Cuadratico = j_cuad, Niveles = j_niv), function(m) {
r <- resettest(m, power = 2:3, type = "fitted"); c(RESET = unname(r$statistic), p = unname(r$p.value))
})## Lineal Cuadratico Niveles
## RESET 2.934894e+01 0.7128819 5.49886501
## p 2.473968e-13 0.4903236 0.00413843