Ya estuve hablando anteriormente de los Metalearners o como se diga aquí y aquí. Pero ahora vamos a ver si lo utilizamos en unos datos reales.
La encuestas de estructura salarial se ha usado muchas veces para ver la brecha salarial entre hombre y mujeres. No obstante yo me hago la pregunta de si es posible y cómo estimar la brecha salarial entre sector público y sector privado.
¿Cómo podríamos hacerlo? Está claro que son dos sectores muy distintos y que comparar las medias, tal y como hacen (mal) algunos para comparar brecha salarial de género, no es lo más adecuado.
Mi idea aquí es contar un poco como lo haríamos estilo compadre usando un modelo lineal de toda la vida, luego ver como se haría utilizando metalearners y también usando doubly robust estimator .
Vamos a utilizar los microdatos de la encuestas de estructura salarial del INE. A pesar de ser una encuesta anual, los últimos resultados publicados son de 2021 y los últimos microdatos disponibles los de 2018. La verdad es que me gustaría entender por qué el INE publica tan tarde los microdatos 😥. La nota de prensa con los resultados de 2021 es del 20 de junio de 2023. Y si ya tienen resultados de 2021, ¿por qué los últimos microdatos disponibles son los de hace 5 años?
Sea como fuere vamos a lo nuestro.
library(tidyverse)
library(haven) # Para leer los datos de SPSS
library(survey) # para obtener estimadores correctos por ser una muestra
library(sjPlot) # plot de los modelos,
# library(causact) # usar alguna cosa de numpyro desde R , quiero probar este paqueteAunque en el fichero comprimido que te descargas del INE viene script para leer los datos con R, no me gusta ese script que instala XLConnect y no sé qué más. Así que lo que he hecho es leer el fichero en formato de spss con haven::read_spss() y luego limpiarlos un poco con janitor::clean_names().
ess <- read_sav(here::here("data/INE/datos_2018/SPSS/EES_2018.sav"))
ess <- janitor::clean_names(ess)
head(ess)
#> # A tibble: 6 × 56
#> idenccc ordentra nuts1 cnace estrato2 control mercado regulacion sexo
#> <chr> <chr> <chr+l> <chr+lb> <chr+lb> <chr+l> <chr+l> <chr+lbl> <chr+l>
#> 1 00000025 01 1 [NOR… H1 [Tra… 1 [DE 1… 2 [PRI… 3 [UNI… 2 [SECTOR… 1 [HOM…
#> 2 00000025 02 1 [NOR… H1 [Tra… 1 [DE 1… 2 [PRI… 3 [UNI… 2 [SECTOR… 1 [HOM…
#> 3 00000025 03 1 [NOR… H1 [Tra… 1 [DE 1… 2 [PRI… 3 [UNI… 2 [SECTOR… 6 [MUJ…
#> 4 00000025 04 1 [NOR… H1 [Tra… 1 [DE 1… 2 [PRI… 3 [UNI… 2 [SECTOR… 1 [HOM…
#> 5 00000025 05 1 [NOR… H1 [Tra… 1 [DE 1… 2 [PRI… 3 [UNI… 2 [SECTOR… 1 [HOM…
#> 6 00000025 06 1 [NOR… H1 [Tra… 1 [DE 1… 2 [PRI… 3 [UNI… 2 [SECTOR… 1 [HOM…
#> # ℹ 47 more variables: tipopais <chr+lbl>, cno1 <chr+lbl>, responsa <chr+lbl>,
#> # estu <chr+lbl>, anoanti <dbl>, mesanti <dbl>, tipojor <chr+lbl>,
#> # tipocon <chr+lbl>, fijodism <dbl>, fijodisd <dbl>, val <dbl>, van <dbl>,
#> # puentes <dbl>, jap <dbl>, jsp1 <dbl>, jsp2 <dbl>, hextra <dbl>,
#> # drelabm <dbl>, siespm1 <chr+lbl>, dsiespm1 <dbl>, siespm2 <chr+lbl>,
#> # dsiespm2 <dbl>, salbase <dbl>, extraorm <dbl>, phextra <dbl>, comsal <dbl>,
#> # comsaltt <dbl>, irpfmes <dbl>, cotiza <dbl>, base <dbl>, drelabam <dbl>, …Dejo enlace al diseño del registro para ver qué es cada variable en los microdatos.
Lo que quiero comparar es el salario neto mensual, ¿por qué? porque me da la gana y porque en la docu del INE explican como se calcula el salario neto partiendo de los microdatos.
ess <- ess |>
mutate(
diasmes = drelabm - dsiespm2,
diasrelaba = drelabam * 30.42 + drelabad,
diasrelaba = ifelse(diasrelaba > 365, 365, diasrelaba),
diasano = diasrelaba - dsiespa2 - dsiespa4,
salbase = ifelse(siespm1 == "6", (31 / diasmes) * salbase, salbase),
comsal = ifelse(siespm1 == "6", (31 / diasmes) * comsal, comsal),
comsaltt = ifelse(siespm1 == "6", (31 / diasmes) * comsaltt, comsaltt),
salmes = salbase + comsal + extraorm + phextra,
salmor = salbase + comsal + phextra,
salneto = salmes - cotiza - irpfmes,
salanual = (365 / diasano) * (retrinoin + retriin + vespnoin + vespin),
salaor = (365 / diasano) * ((retrinoin + retriin) - gextra),
vespnoin = (365 / diasano) * vespnoin,
jmp1 = (jsp1 + jsp2 / 60) * 4.35 + hextra,
salhora = salmes / jmp1
)La variable dónde se consigna si el sector es público o privado es control
ess |>
group_by(control) |>
count()
#> # A tibble: 2 × 2
#> # Groups: control [2]
#> control n
#> <chr+lbl> <int>
#> 1 1 [PUBLICO] 35553
#> 2 2 [PRIVADO] 181173Voy a crearme variable treatment que valga 1 cuando sea sector público y 0 para el sector privado
ess$treatment = ess$control
ess$treatment = ifelse(ess$control == "1", 1, 0)
# también llamo outcome al salario neto
ess$outcome = ess$salnetoLo más simple , hacemos un group by y calculamos medias
ess |>
group_by(treatment) |>
summarise(
mean = mean(outcome),
n = n()
)
#> # A tibble: 2 × 3
#> treatment mean n
#> <dbl> <dbl> <int>
#> 1 0 1561. 181173
#> 2 1 1864. 35553Así de primeras, pues parece que se gana más en el sector público que en el privado, pero ¡ojo! que la encuesta tiene una variable de ponderación, que el INE ha calculado para que los resultados sean representativos de la población. En la nota metodológica el INE dice lo siguiente sobre el plan de muestreo
El procedimiento de selección aleatoria de unidades corresponde a un muestreo bietápico estratificado, donde las unidades de primera etapa son las cuentas de cotización a la Seguridad Social (CC), mientras que las de segunda etapa son los trabajadores (asalariados). En la primera etapa tanto el diseño muestral como la muestra obtenida de CC coincide con la ETCL (para una mayor información consultar la metodología de la ETCL). Las unidades de primera etapa se estratifican según las siguientes variables:
En los microdatos tenemos la variable factotal que es la ponderación que el INE dice que hay que usar a la hora de hacer estimaciones.
ess |>
group_by(treatment) |>
summarise(
media_ponderada = weighted.mean(outcome, w = factotal)
)
#> # A tibble: 2 × 2
#> treatment media_ponderada
#> <dbl> <dbl>
#> 1 0 1351.
#> 2 1 1799.Pero sabemos que la media tal cual puede no ser un buen indicador, lo suyo sería controlar (condicionar) por otras variables, tales como el sexo, nivel de estudio, edad, años de antigüedad , tipo de jornada laboral, y cosas así.
Hagámoslo, pero usando que tenemos pesos en la encuesta.
disenno <- svydesign(id = ~1, weight = ~factotal, data = ess)Modelo simple dónde uso variables como edad, tipo contrato, área nuts, antigüedad, nivel de estudios, etc..
mod_simple <- svyglm(outcome ~ treatment + sexo + anos2 + estu + cnace + cno1 + estrato2 + tipojor + anoanti + mesanti + tipocon + nuts1, design = disenno)summary(mod_simple)
#>
#> Call:
#> svyglm(formula = outcome ~ treatment + sexo + anos2 + estu +
#> cnace + cno1 + estrato2 + tipojor + anoanti + mesanti + tipocon +
#> nuts1, design = disenno)
#>
#> Survey design:
#> svydesign(id = ~1, weight = ~factotal, data = ess)
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2267.4873 58.5408 38.733 < 2e-16 ***
#> treatment 112.7778 8.1697 13.804 < 2e-16 ***
#> sexo6 -158.0960 5.4276 -29.128 < 2e-16 ***
#> anos202 21.0286 25.4780 0.825 0.409167
#> anos203 101.1315 25.4521 3.973 7.09e-05 ***
#> anos204 149.8967 25.3798 5.906 3.51e-09 ***
#> anos205 176.1458 25.6102 6.878 6.09e-12 ***
#> anos206 52.5466 27.0543 1.942 0.052106 .
#> estu2 29.6478 14.4793 2.048 0.040601 *
#> estu3 38.1183 14.3801 2.651 0.008031 **
#> estu4 142.3435 14.8437 9.590 < 2e-16 ***
#> estu5 166.5131 15.9282 10.454 < 2e-16 ***
#> estu6 265.4148 17.3323 15.313 < 2e-16 ***
#> estu7 510.3787 18.1827 28.070 < 2e-16 ***
#> cnaceC1 -324.1132 25.8960 -12.516 < 2e-16 ***
#> cnaceC2 -278.7694 31.4639 -8.860 < 2e-16 ***
#> cnaceC3 -366.8802 28.5581 -12.847 < 2e-16 ***
#> cnaceC4 -107.0397 26.9400 -3.973 7.09e-05 ***
#> cnaceC5 -252.2397 27.0669 -9.319 < 2e-16 ***
#> cnaceC6 -167.9342 26.4826 -6.341 2.28e-10 ***
#> cnaceC7 -180.7201 28.5478 -6.330 2.45e-10 ***
#> cnaceC8 -171.5305 25.6146 -6.697 2.14e-11 ***
#> cnaceD0 323.7104 36.5705 8.852 < 2e-16 ***
#> cnaceE0 -282.7253 26.5709 -10.640 < 2e-16 ***
#> cnaceF0 -222.3453 25.6119 -8.681 < 2e-16 ***
#> cnaceG1 -266.4662 27.1602 -9.811 < 2e-16 ***
#> cnaceG2 -409.9509 28.0075 -14.637 < 2e-16 ***
#> cnaceH1 -231.9251 27.2267 -8.518 < 2e-16 ***
#> cnaceH2 -288.5758 26.5848 -10.855 < 2e-16 ***
#> cnaceI0 -314.6540 26.9575 -11.672 < 2e-16 ***
#> cnaceJ0 -305.2619 27.9572 -10.919 < 2e-16 ***
#> cnaceK0 -26.9416 29.4084 -0.916 0.359606
#> cnaceL0 -426.2992 31.8821 -13.371 < 2e-16 ***
#> cnaceM0 -382.3510 26.5815 -14.384 < 2e-16 ***
#> cnaceN0 -426.7819 26.7070 -15.980 < 2e-16 ***
#> cnaceO0 -489.0672 26.5529 -18.419 < 2e-16 ***
#> cnaceP0 -677.0378 28.3344 -23.895 < 2e-16 ***
#> cnaceQ0 -417.2582 26.2085 -15.921 < 2e-16 ***
#> cnaceR0 -446.4843 27.0788 -16.488 < 2e-16 ***
#> cnaceS0 -441.9862 25.9167 -17.054 < 2e-16 ***
#> cno1B0 -564.2450 40.6680 -13.874 < 2e-16 ***
#> cno1C0 -581.5582 40.3270 -14.421 < 2e-16 ***
#> cno1D0 -772.4413 38.8858 -19.864 < 2e-16 ***
#> cno1E0 -964.0718 38.7102 -24.905 < 2e-16 ***
#> cno1F0 -963.2193 39.0883 -24.642 < 2e-16 ***
#> cno1G0 -944.6769 40.6941 -23.214 < 2e-16 ***
#> cno1H0 -1007.1839 39.2405 -25.667 < 2e-16 ***
#> cno1I0 -805.5218 41.4801 -19.419 < 2e-16 ***
#> cno1J0 -833.4489 50.5645 -16.483 < 2e-16 ***
#> cno1K0 -971.3301 39.7779 -24.419 < 2e-16 ***
#> cno1L0 -948.3725 38.9086 -24.374 < 2e-16 ***
#> cno1M0 -977.4409 38.9573 -25.090 < 2e-16 ***
#> cno1N0 -976.4017 40.1926 -24.293 < 2e-16 ***
#> cno1O0 -1009.1702 40.1926 -25.108 < 2e-16 ***
#> cno1P0 -1033.0883 39.0592 -26.449 < 2e-16 ***
#> cno1Q0 -981.9190 65.8301 -14.916 < 2e-16 ***
#> estrato21 25.7381 28.0871 0.916 0.359476
#> estrato22 130.3521 28.3970 4.590 4.43e-06 ***
#> estrato23 208.3260 28.3907 7.338 2.18e-13 ***
#> estrato24 130.3223 29.7704 4.378 1.20e-05 ***
#> tipojor2 -593.8166 5.3784 -110.408 < 2e-16 ***
#> anoanti 12.3393 0.3434 35.929 < 2e-16 ***
#> mesanti 2.0881 0.6091 3.428 0.000608 ***
#> tipocon2 -119.7352 5.0449 -23.734 < 2e-16 ***
#> nuts12 139.1728 8.4611 16.449 < 2e-16 ***
#> nuts13 97.8708 9.3378 10.481 < 2e-16 ***
#> nuts14 -4.8016 7.8750 -0.610 0.542040
#> nuts15 104.8487 7.5335 13.918 < 2e-16 ***
#> nuts16 48.8032 8.2790 5.895 3.76e-09 ***
#> nuts17 17.5482 10.3777 1.691 0.090848 .
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for gaussian family taken to be 373331.4)
#>
#> Number of Fisher Scoring iterations: 2Y el coeficiente asociaso al sector público indica que se gana en media unos 112 euros más que en el sector privado, según este modelo.
¿Cuánto sería para alguien que trabaja a jornada completa, nivel de estudios superior o igual a licenciado?
Para eso podemos hacer lo que se conoce como una “intervención”, que es crear dos conjuntos de datos copias del original, con la diferencia de que en uno todo el mundo es del sector privado y en el otro todos del sector público y comparamos las medias estimadas de salario neto que nos da el modelo para el subgrupo de población que queramos.
A esto se le conoce por los modernos como un S-learner
ess_fake_publico <- ess
ess_fake_publico$treatment <- 1
estim_publico <- predict(mod_simple, newdata = ess_fake_publico |> filter(tipojor == "1", estu == "7", sexo == "1" ) )
ess_fake_privado <- ess
ess_fake_privado$treatment <- 0
estim_privado <- predict(mod_simple, newdata = ess_fake_privado |> filter(tipojor == "1", estu == "7", sexo == "1") )
mean(estim_publico)
#> [1] 2516.276
mean(estim_privado)
#> [1] 2403.498
(s_learner_with_pond <- mean(estim_publico) - mean(estim_privado))
#> [1] 112.7778Y coincide con el coeficiente que daba el modelo. Y eso es así porque no he metido interacciones en el modelo. Si metemos una simple interacción entre ser del sector público y privado con la zona Nuts1.
ess |>
group_by(nuts1) |>
count()
#> # A tibble: 7 × 2
#> # Groups: nuts1 [7]
#> nuts1 n
#> <chr+lbl> <int>
#> 1 1 [NOROESTE] 24806
#> 2 2 [NORESTE] 33624
#> 3 3 [COMUNIDAD DE MADRID] 34269
#> 4 4 [CENTRO] 26428
#> 5 5 [ESTE] 58852
#> 6 6 [SUR] 29413
#> 7 7 [CANARIAS] 9334mod_inter_con_nuts <- svyglm(outcome ~ treatment* nuts1 + sexo + anos2 + estu + cnace + cno1 + estrato2 + tipojor + anoanti + mesanti + tipocon , design = disenno)
summary(mod_inter_con_nuts)
#>
#> Call:
#> svyglm(formula = outcome ~ treatment * nuts1 + sexo + anos2 +
#> estu + cnace + cno1 + estrato2 + tipojor + anoanti + mesanti +
#> tipocon, design = disenno)
#>
#> Survey design:
#> svydesign(id = ~1, weight = ~factotal, data = ess)
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2267.0375 58.7580 38.583 < 2e-16 ***
#> treatment 139.5808 16.4879 8.466 < 2e-16 ***
#> nuts12 136.3734 9.6420 14.144 < 2e-16 ***
#> nuts13 109.8222 10.5875 10.373 < 2e-16 ***
#> nuts14 -13.0969 8.7817 -1.491 0.135862
#> nuts15 115.8699 8.4975 13.636 < 2e-16 ***
#> nuts16 50.2008 9.4214 5.328 9.92e-08 ***
#> nuts17 11.3094 11.6282 0.973 0.330765
#> sexo6 -157.8517 5.4260 -29.092 < 2e-16 ***
#> anos202 22.9645 25.5775 0.898 0.369272
#> anos203 103.1094 25.5532 4.035 5.46e-05 ***
#> anos204 151.6792 25.4825 5.952 2.65e-09 ***
#> anos205 177.3203 25.7096 6.897 5.32e-12 ***
#> anos206 53.6148 27.1482 1.975 0.048282 *
#> estu2 29.8615 14.4816 2.062 0.039207 *
#> estu3 37.3601 14.3817 2.598 0.009384 **
#> estu4 142.6627 14.8458 9.610 < 2e-16 ***
#> estu5 167.2457 15.9296 10.499 < 2e-16 ***
#> estu6 264.4340 17.3311 15.258 < 2e-16 ***
#> estu7 510.9210 18.1818 28.101 < 2e-16 ***
#> cnaceC1 -322.4719 25.8565 -12.472 < 2e-16 ***
#> cnaceC2 -277.8307 31.4386 -8.837 < 2e-16 ***
#> cnaceC3 -368.1116 28.5083 -12.912 < 2e-16 ***
#> cnaceC4 -107.5481 26.8992 -3.998 6.38e-05 ***
#> cnaceC5 -251.5735 27.0238 -9.309 < 2e-16 ***
#> cnaceC6 -165.5994 26.4505 -6.261 3.84e-10 ***
#> cnaceC7 -179.9980 28.5076 -6.314 2.72e-10 ***
#> cnaceC8 -170.7310 25.5703 -6.677 2.45e-11 ***
#> cnaceD0 323.6055 36.5367 8.857 < 2e-16 ***
#> cnaceE0 -281.1503 26.5217 -10.601 < 2e-16 ***
#> cnaceF0 -221.8244 25.5641 -8.677 < 2e-16 ***
#> cnaceG1 -266.6874 27.1071 -9.838 < 2e-16 ***
#> cnaceG2 -410.0799 27.9517 -14.671 < 2e-16 ***
#> cnaceH1 -228.7886 27.1892 -8.415 < 2e-16 ***
#> cnaceH2 -287.2581 26.5307 -10.827 < 2e-16 ***
#> cnaceI0 -314.4745 26.9273 -11.679 < 2e-16 ***
#> cnaceJ0 -307.8448 27.9093 -11.030 < 2e-16 ***
#> cnaceK0 -25.9796 29.3654 -0.885 0.376318
#> cnaceL0 -428.7358 31.8459 -13.463 < 2e-16 ***
#> cnaceM0 -382.9944 26.5201 -14.442 < 2e-16 ***
#> cnaceN0 -427.5663 26.6470 -16.046 < 2e-16 ***
#> cnaceO0 -491.1736 26.4855 -18.545 < 2e-16 ***
#> cnaceP0 -679.4070 28.3151 -23.995 < 2e-16 ***
#> cnaceQ0 -418.4000 26.1469 -16.002 < 2e-16 ***
#> cnaceR0 -446.9767 27.0305 -16.536 < 2e-16 ***
#> cnaceS0 -441.7225 25.8682 -17.076 < 2e-16 ***
#> cno1B0 -563.9700 40.6603 -13.870 < 2e-16 ***
#> cno1C0 -580.4722 40.3078 -14.401 < 2e-16 ***
#> cno1D0 -771.5835 38.8634 -19.854 < 2e-16 ***
#> cno1E0 -963.6523 38.6908 -24.907 < 2e-16 ***
#> cno1F0 -962.4620 39.0754 -24.631 < 2e-16 ***
#> cno1G0 -943.1722 40.6791 -23.186 < 2e-16 ***
#> cno1H0 -1006.9339 39.2172 -25.676 < 2e-16 ***
#> cno1I0 -799.1017 41.4442 -19.281 < 2e-16 ***
#> cno1J0 -833.4729 50.5531 -16.487 < 2e-16 ***
#> cno1K0 -970.9199 39.7730 -24.412 < 2e-16 ***
#> cno1L0 -946.7735 38.8851 -24.348 < 2e-16 ***
#> cno1M0 -975.6476 38.9286 -25.062 < 2e-16 ***
#> cno1N0 -974.9346 40.1732 -24.268 < 2e-16 ***
#> cno1O0 -1008.8982 40.1716 -25.115 < 2e-16 ***
#> cno1P0 -1032.8826 39.0379 -26.458 < 2e-16 ***
#> cno1Q0 -941.5737 66.0825 -14.248 < 2e-16 ***
#> estrato21 18.7235 28.1838 0.664 0.506476
#> estrato22 122.1463 28.5034 4.285 1.83e-05 ***
#> estrato23 201.3305 28.4886 7.067 1.59e-12 ***
#> estrato24 126.1192 29.8559 4.224 2.40e-05 ***
#> tipojor2 -593.1639 5.3768 -110.318 < 2e-16 ***
#> anoanti 12.3334 0.3434 35.921 < 2e-16 ***
#> mesanti 2.0871 0.6094 3.425 0.000615 ***
#> tipocon2 -119.0692 5.0451 -23.601 < 2e-16 ***
#> treatment:nuts12 15.2013 19.3060 0.787 0.431056
#> treatment:nuts13 -71.6951 20.4396 -3.508 0.000452 ***
#> treatment:nuts14 27.4429 19.1259 1.435 0.151329
#> treatment:nuts15 -66.6916 17.9164 -3.722 0.000197 ***
#> treatment:nuts16 -8.9706 19.3186 -0.464 0.642398
#> treatment:nuts17 39.0070 25.1314 1.552 0.120634
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for gaussian family taken to be 373097.8)
#>
#> Number of Fisher Scoring iterations: 2Estimamos diferencias entre sector público y privado para Madrid y Andalucía, para un hombre a jornada completa y con estudios de licenciatura o superior.
estim_publico_madrid <- predict(mod_inter_con_nuts, newdata = ess_fake_publico |> filter(tipojor == "1", estu == "7", sexo == "1" , nuts1 == "3" ) )
estim_privado_madrid <- predict(mod_inter_con_nuts, newdata = ess_fake_privado |> filter(tipojor == "1", estu == "7", sexo == "1", nuts1 == "3") )
mean(estim_publico_madrid)
#> [1] 2533.962
mean(estim_privado_madrid)
#> [1] 2466.076
(s_learner_with_pond_madrid <- mean(estim_publico_madrid) - mean(estim_privado_madrid))
#> [1] 67.88572
estim_publico_sur <- predict(mod_inter_con_nuts, newdata = ess_fake_publico |> filter(tipojor == "1", estu == "7", sexo == "1" , nuts1 == "1" ) )
estim_privado_sur <- predict(mod_inter_con_nuts, newdata = ess_fake_privado |> filter(tipojor == "1", estu == "7", sexo == "1", nuts1 == "1") )
mean(estim_publico_sur)
#> [1] 2449.94
mean(estim_privado_sur)
#> [1] 2310.359
(s_learner_with_pond_sur <- mean(estim_publico_sur) - mean(estim_privado_sur))
#> [1] 139.5808Bueno, pues según esto, parece que para ese perfil, dónde se ha tenido en cuenta edad, años de antigüedad y demás, se gana un poco más en el sector público que en el privado, aunque esa diferencia es mayor en el Sur que en Madrid.
Otro de los metalearners empleados es el T-learner, ya explicado en post anteriores. Aquí vamos a usarlo sin tener en cuenta la ponderación de la encuesta.
En el T-learner se ajusta un modelo para cuando sea sector público y otro para cuando sea sector privado y se ve la diferencia de las medias de sus estimaciones.
modpublico <- lm(outcome ~ sexo + anos2 + estu + cno1 + estrato2 + tipojor + anoanti + mesanti + tipocon + nuts1, data = ess[ess$treatment==1, ])
modprivado <- lm(outcome ~ sexo + anos2 + estu + cno1 + estrato2 + tipojor + anoanti + mesanti + tipocon + nuts1, data = ess[ess$treatment==0, ])
ess_sub <- ess %>%
filter(tipojor == "1", estu == "7", sexo == "1")
# t-learner
(t_learner <- mean(predict(modpublico, ess_sub)) -
mean(predict(modprivado, ess_sub)) )
#> [1] -175.951
ess_sub_madrid <- ess %>%
filter(tipojor == "1", estu == "7", sexo == "1", nuts1 == "3")
# t-learner
(t_learner_madrid <- mean(predict(modpublico, ess_sub_madrid)) -
mean(predict(modprivado, ess_sub_madrid)) )
#> [1] -296.4373
ess_sub_sur <- ess %>%
filter(tipojor == "1", estu == "7", sexo == "1", nuts1 == "1")
# t-learner
(t_learner_sur <- mean(predict(modpublico, ess_sub_sur)) -
mean(predict(modprivado, ess_sub_sur)) )
#> [1] -159.0581En este caso, nos sale que se ganaría más en el sector público que en el privado. ¿Con qué nos quedamos?
Ya expliqué en su día en que consiste un X-learner (../2020/12/30/y-si-parte-ii/#x-learner)
Básicamente, usas el modelo ajustado con treatment = 1 para predecir las observaciones con treatment = 0 y al revés en un intento de estimar el potential outcome. Luego haces dos modelos para modelar las diferencias entre el outcome y las predicciones anteriores y otro modelo de propensity score que se usará para ponderar esas dos predicciones.
m1 <- lm(outcome ~ sexo + anos2 + estu + estrato2 + tipojor + anoanti + mesanti + tipocon + nuts1, data = ess[ess$treatment==1, ])
m2 <- lm(outcome ~ sexo + anos2 + estu + estrato2 + tipojor + anoanti + mesanti + tipocon + nuts1, data = ess[ess$treatment==0, ])
# Usamos modelo 1 para estimar cuando W=0 y el modelo 2 para estimar cuando W = 1
# Con el viejo R-base sería
ess$Difer[ess$treatment==0] <- ess$outcome[ess$treatment==0] - predict(m1, ess[ess$treatment==0, ])
head(ess[ess$treatment==0, c("outcome", "Difer")])
#> # A tibble: 6 × 2
#> outcome Difer
#> <dbl> <dbl>
#> 1 302. -267.
#> 2 1181. -263.
#> 3 1314. -784.
#> 4 1210. -285.
#> 5 1152. -201.
#> 6 1184. -237.
ess$Difer[ess$treatment==1] <- ess$outcome[ess$treatment==1] - predict(m2, ess[ess$treatment==1, ])
head(ess[ess$treatment==1, c("outcome", "Difer")])
#> # A tibble: 6 × 2
#> outcome Difer
#> <dbl> <dbl>
#> 1 2112. 288.
#> 2 1989. -405.
#> 3 2183. -36.2
#> 4 1954. -212.
#> 5 1779. -483.
#> 6 1738. -447.
# Modelamos las diferencias
m3 <- lm(Difer ~ sexo + anos2 + estu + estrato2 + tipojor + anoanti + mesanti + tipocon + nuts1, data = ess[ess$treatment==1, ])
m4 <- lm(Difer ~ sexo + anos2 + estu + estrato2 + tipojor + anoanti + mesanti + tipocon + nuts1, data = ess[ess$treatment==0, ])
# Combinamos
glm1 <- glm(treatment ~ sexo + anos2 + estu + estrato2 + tipojor + anoanti + mesanti + tipocon + nuts1, data = ess, family= binomial)
ess$pesos <- predict(glm1, ess, type = "response")
ess$combinado <- ess$pesos * predict(m4, ess) + (1-ess$pesos) * predict(m3, ess)
head(ess[, c("outcome", "treatment","Difer", "pesos", "combinado")])
#> # A tibble: 6 × 5
#> outcome treatment Difer pesos combinado
#> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 302. 0 -267. 0.111 -254.
#> 2 1181. 0 -263. 0.0238 125.
#> 3 1314. 0 -784. 0.0591 93.2
#> 4 1210. 0 -285. 0.0332 74.0
#> 5 1152. 0 -201. 0.0107 -19.7
#> 6 1184. 0 -237. 0.0201 -62.1
(x_learner <- ess %>%
filter(tipojor == "1", estu == "7", sexo == "1") |>
group_by(treatment) %>%
summarise(mean = mean(outcome)) |>
pivot_wider(names_from = treatment, values_from = mean, names_prefix = "mean_") |>
mutate(
estim_xlearner = mean_1 - mean_0) |>
pull(estim_xlearner))
#> [1] -228.745
(x_learner_madrid <- ess %>%
filter(tipojor == "1", estu == "7", sexo == "1", nuts1=="3") |>
group_by(treatment) %>%
summarise(mean = mean(outcome)) |>
pivot_wider(names_from = treatment, values_from = mean, names_prefix = "mean_") |>
mutate(
estim_xlearner = mean_1 - mean_0) |>
pull(estim_xlearner))
#> [1] -501.2355
(x_learner_sur <- ess %>%
filter(tipojor == "1", estu == "7", sexo == "1", nuts1=="1") |>
group_by(treatment) %>%
summarise(mean = mean(outcome)) |>
pivot_wider(names_from = treatment, values_from = mean, names_prefix = "mean_") |>
mutate(
estim_xlearner = mean_1 - mean_0) |>
pull(estim_xlearner))
#> [1] -55.67715Con idea parecida al X-learner , en el sentido de mezclar las estrategias de usar Inverse probability weighting y el de hacer un modelo de la respuesta condicionando por los counfounders.
De nuevo, al igual que con el T-Learner o el X-Learner no vamos a tener en cuenta la variable de ponderación de casos.
El estimador sería algo así como
\[\dfrac{1}{n} \sum_{i=1}^n \left[ \dfrac{Y_i \cdot A_i - \color{red}{ \left(A_i -\pi(X_i)\right) \mu(X_i, A_i)})} {\pi(X_i)} - \dfrac{Y_i \cdot (1-A_i) - \color{red}{ \left(A_i -\pi(X_i)\right) \mu(X_i,A_i)})} {1-\pi(X_i)} \right ] \tag{1}\]
Dónde \(\mu\) hace referencia al modelo para estimar el outcome y \(\pi\) al modelo de propensity score.
Este Doubly robust estimator es una combinación entre usar inverse probability weighting y el modelo de la media del outcome. Este estimador suele ser consistnete si al menos uno de los dos modelos es correcto. A la expresión coloreda en rojo se le denomina augmented IPW estimator
En código es bastante sencillo.
dr_estimator <- function(data, prop_model, mean_model){
data %>%
mutate(
prob = predict(prop_model, newdata = data, type = "response"),
pred = predict(mean_model, newdata = data, type = "response"),
augm = (treatment - prob) * pred
) %>%
summarise(
EYpublico = mean((outcome * treatment -augm) / prob),
EYprivado= mean((outcome * (1 - treatment) - augm) / (1 - prob))
) %>%
mutate(dre = EYpublico - EYprivado)
}Y si usamos ese estimador tenemos
prop_model <- glm(treatment ~ sexo + anos2+ cnace + cno1 + estrato2 + estu + tipojor + anoanti + mesanti + tipocon + nuts1, data = ess, family = binomial)
mean_model <- glm(outcome ~ treatment + sexo + anos2+ cnace + cno1 + estrato2 + estu + tipojor + anoanti + mesanti + tipocon + nuts1 , data = ess, family = gaussian)
summary(prop_model)
#>
#> Call:
#> glm(formula = treatment ~ sexo + anos2 + cnace + cno1 + estrato2 +
#> estu + tipojor + anoanti + mesanti + tipocon + nuts1, family = binomial,
#> data = ess)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -8.090e+00 5.934e-01 -13.634 < 2e-16 ***
#> sexo6 2.674e-01 2.116e-02 12.637 < 2e-16 ***
#> anos202 8.244e-01 5.243e-01 1.572 0.115866
#> anos203 1.488e+00 5.233e-01 2.844 0.004459 **
#> anos204 1.817e+00 5.233e-01 3.473 0.000515 ***
#> anos205 1.981e+00 5.235e-01 3.785 0.000154 ***
#> anos206 2.064e+00 5.242e-01 3.937 8.24e-05 ***
#> cnaceC1 -2.069e+00 1.909e-01 -10.835 < 2e-16 ***
#> cnaceC2 -2.340e+00 2.854e-01 -8.200 2.41e-16 ***
#> cnaceC3 4.823e-01 1.773e-01 2.720 0.006533 **
#> cnaceC4 -1.680e+01 1.053e+02 -0.160 0.873244
#> cnaceC5 -1.634e+01 1.897e+02 -0.086 0.931349
#> cnaceC6 -2.705e+00 2.597e-01 -10.416 < 2e-16 ***
#> cnaceC7 -1.649e+01 1.280e+02 -0.129 0.897495
#> cnaceC8 -6.740e-01 1.545e-01 -4.363 1.28e-05 ***
#> cnaceD0 -5.031e+00 7.223e-01 -6.964 3.30e-12 ***
#> cnaceE0 1.948e+00 1.483e-01 13.139 < 2e-16 ***
#> cnaceF0 2.271e-01 1.511e-01 1.503 0.132789
#> cnaceG1 -1.622e+01 1.124e+02 -0.144 0.885252
#> cnaceG2 -1.600e+01 9.413e+01 -0.170 0.865012
#> cnaceH1 1.437e+00 1.507e-01 9.534 < 2e-16 ***
#> cnaceH2 2.936e+00 1.482e-01 19.808 < 2e-16 ***
#> cnaceI0 -7.434e-01 1.783e-01 -4.170 3.04e-05 ***
#> cnaceJ0 3.013e-01 1.472e-01 2.047 0.040683 *
#> cnaceK0 1.752e-02 1.484e-01 0.118 0.905995
#> cnaceL0 6.883e-01 1.787e-01 3.852 0.000117 ***
#> cnaceM0 1.066e+00 1.455e-01 7.326 2.38e-13 ***
#> cnaceN0 -7.697e-01 1.550e-01 -4.966 6.85e-07 ***
#> cnaceO0 2.246e+01 9.488e+01 0.237 0.812845
#> cnaceP0 2.748e+00 1.487e-01 18.482 < 2e-16 ***
#> cnaceQ0 2.013e+00 1.464e-01 13.750 < 2e-16 ***
#> cnaceR0 1.759e+00 1.487e-01 11.834 < 2e-16 ***
#> cnaceS0 5.370e-01 1.554e-01 3.455 0.000550 ***
#> cno1B0 1.487e+00 6.226e-02 23.887 < 2e-16 ***
#> cno1C0 8.397e-01 5.779e-02 14.529 < 2e-16 ***
#> cno1D0 4.886e-01 5.826e-02 8.386 < 2e-16 ***
#> cno1E0 3.641e-01 6.314e-02 5.767 8.08e-09 ***
#> cno1F0 2.401e-02 6.783e-02 0.354 0.723379
#> cno1G0 -2.950e-01 1.191e-01 -2.476 0.013269 *
#> cno1H0 4.008e-01 6.743e-02 5.944 2.79e-09 ***
#> cno1I0 5.683e-01 1.077e-01 5.277 1.31e-07 ***
#> cno1J0 1.599e+00 1.437e-01 11.121 < 2e-16 ***
#> cno1K0 4.359e-01 9.292e-02 4.691 2.72e-06 ***
#> cno1L0 2.333e-01 7.777e-02 2.999 0.002705 **
#> cno1M0 -8.548e-01 1.435e-01 -5.956 2.59e-09 ***
#> cno1N0 2.546e-01 7.800e-02 3.264 0.001098 **
#> cno1O0 3.076e-01 7.608e-02 4.043 5.27e-05 ***
#> cno1P0 -8.458e-01 9.808e-02 -8.624 < 2e-16 ***
#> cno1Q0 3.002e+00 1.482e+03 0.002 0.998384
#> estrato21 -9.996e-01 2.206e-01 -4.531 5.87e-06 ***
#> estrato22 9.719e-02 2.202e-01 0.441 0.658908
#> estrato23 9.132e-01 2.202e-01 4.147 3.37e-05 ***
#> estrato24 5.801e-01 2.244e-01 2.585 0.009724 **
#> estu2 2.584e-01 1.589e-01 1.626 0.103944
#> estu3 1.217e+00 1.575e-01 7.727 1.10e-14 ***
#> estu4 1.253e+00 1.582e-01 7.917 2.43e-15 ***
#> estu5 1.644e+00 1.600e-01 10.275 < 2e-16 ***
#> estu6 1.470e+00 1.607e-01 9.145 < 2e-16 ***
#> estu7 1.858e+00 1.603e-01 11.593 < 2e-16 ***
#> tipojor2 -1.526e+00 3.257e-02 -46.857 < 2e-16 ***
#> anoanti 7.184e-02 1.209e-03 59.437 < 2e-16 ***
#> mesanti 3.078e-02 2.627e-03 11.718 < 2e-16 ***
#> tipocon2 1.800e+00 2.660e-02 67.697 < 2e-16 ***
#> nuts12 -7.224e-02 3.819e-02 -1.891 0.058560 .
#> nuts13 -5.637e-01 3.770e-02 -14.952 < 2e-16 ***
#> nuts14 5.716e-01 3.905e-02 14.636 < 2e-16 ***
#> nuts15 -3.778e-01 3.500e-02 -10.795 < 2e-16 ***
#> nuts16 3.272e-01 3.778e-02 8.660 < 2e-16 ***
#> nuts17 1.661e-02 5.427e-02 0.306 0.759563
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#> Null deviance: 193458 on 216725 degrees of freedom
#> Residual deviance: 76006 on 216657 degrees of freedom
#> AIC: 76144
#>
#> Number of Fisher Scoring iterations: 18
summary(mean_model)
#>
#> Call:
#> glm(formula = outcome ~ treatment + sexo + anos2 + cnace + cno1 +
#> estrato2 + estu + tipojor + anoanti + mesanti + tipocon +
#> nuts1, family = gaussian, data = ess)
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 2877.5791 61.1896 47.027 < 2e-16 ***
#> treatment 17.9283 8.3590 2.145 0.03197 *
#> sexo6 -240.4080 4.9123 -48.940 < 2e-16 ***
#> anos202 22.1784 43.0544 0.515 0.60647
#> anos203 89.1187 42.8391 2.080 0.03750 *
#> anos204 178.3491 42.8356 4.164 3.13e-05 ***
#> anos205 240.2115 42.9841 5.588 2.29e-08 ***
#> anos206 99.7613 43.5939 2.288 0.02211 *
#> cnaceC1 -362.2154 26.7350 -13.548 < 2e-16 ***
#> cnaceC2 -336.6893 29.3593 -11.468 < 2e-16 ***
#> cnaceC3 -407.7262 33.1492 -12.300 < 2e-16 ***
#> cnaceC4 -142.3327 26.9322 -5.285 1.26e-07 ***
#> cnaceC5 -320.6788 30.9230 -10.370 < 2e-16 ***
#> cnaceC6 -214.7569 28.1272 -7.635 2.26e-14 ***
#> cnaceC7 -295.1496 27.9070 -10.576 < 2e-16 ***
#> cnaceC8 -304.6001 26.3723 -11.550 < 2e-16 ***
#> cnaceD0 315.2232 33.1870 9.498 < 2e-16 ***
#> cnaceE0 -368.5855 27.9188 -13.202 < 2e-16 ***
#> cnaceF0 -251.8171 26.4790 -9.510 < 2e-16 ***
#> cnaceG1 -317.1349 27.2848 -11.623 < 2e-16 ***
#> cnaceG2 -443.6205 28.1035 -15.785 < 2e-16 ***
#> cnaceH1 -24.2486 28.2536 -0.858 0.39076
#> cnaceH2 -252.9940 28.3037 -8.939 < 2e-16 ***
#> cnaceI0 -364.2829 27.7660 -13.120 < 2e-16 ***
#> cnaceJ0 -399.0452 26.5212 -15.046 < 2e-16 ***
#> cnaceK0 -184.2903 27.1146 -6.797 1.07e-11 ***
#> cnaceL0 -376.1283 34.1175 -11.025 < 2e-16 ***
#> cnaceM0 -481.0621 26.1363 -18.406 < 2e-16 ***
#> cnaceN0 -465.0980 26.1408 -17.792 < 2e-16 ***
#> cnaceO0 -498.8197 27.9612 -17.840 < 2e-16 ***
#> cnaceP0 -785.3254 28.6574 -27.404 < 2e-16 ***
#> cnaceQ0 -448.4801 26.8000 -16.734 < 2e-16 ***
#> cnaceR0 -286.1541 27.6853 -10.336 < 2e-16 ***
#> cnaceS0 -523.3793 28.0319 -18.671 < 2e-16 ***
#> cno1B0 -1061.3503 16.4387 -64.564 < 2e-16 ***
#> cno1C0 -1133.4003 13.6831 -82.832 < 2e-16 ***
#> cno1D0 -1275.7747 13.3916 -95.267 < 2e-16 ***
#> cno1E0 -1533.8695 14.7719 -103.837 < 2e-16 ***
#> cno1F0 -1503.8142 15.1869 -99.020 < 2e-16 ***
#> cno1G0 -1478.8524 17.2092 -85.934 < 2e-16 ***
#> cno1H0 -1541.8639 16.5003 -93.444 < 2e-16 ***
#> cno1I0 -1498.1010 20.0039 -74.890 < 2e-16 ***
#> cno1J0 -1517.7995 36.9359 -41.093 < 2e-16 ***
#> cno1K0 -1527.8727 18.7324 -81.563 < 2e-16 ***
#> cno1L0 -1494.9651 15.1867 -98.439 < 2e-16 ***
#> cno1M0 -1549.5391 16.7938 -92.269 < 2e-16 ***
#> cno1N0 -1589.9336 17.9925 -88.367 < 2e-16 ***
#> cno1O0 -1570.6894 16.5873 -94.692 < 2e-16 ***
#> cno1P0 -1595.2569 16.1226 -98.945 < 2e-16 ***
#> cno1Q0 -1530.0462 138.3092 -11.063 < 2e-16 ***
#> estrato21 -26.9047 30.6473 -0.878 0.38001
#> estrato22 93.7567 30.8162 3.042 0.00235 **
#> estrato23 189.6530 30.7961 6.158 7.36e-10 ***
#> estrato24 99.0138 33.0210 2.999 0.00271 **
#> estu2 17.6977 22.9002 0.773 0.43963
#> estu3 31.8206 22.8239 1.394 0.16326
#> estu4 146.0304 23.0680 6.330 2.45e-10 ***
#> estu5 188.1811 23.7531 7.922 2.34e-15 ***
#> estu6 340.0061 24.1380 14.086 < 2e-16 ***
#> estu7 629.4307 23.9967 26.230 < 2e-16 ***
#> tipojor2 -632.4750 6.1453 -102.921 < 2e-16 ***
#> anoanti 13.7806 0.2783 49.511 < 2e-16 ***
#> mesanti 1.9103 0.5911 3.232 0.00123 **
#> tipocon2 -105.5518 5.9315 -17.795 < 2e-16 ***
#> nuts12 112.1766 8.2021 13.677 < 2e-16 ***
#> nuts13 159.6686 8.4096 18.986 < 2e-16 ***
#> nuts14 -0.7868 8.6804 -0.091 0.92778
#> nuts15 115.6430 7.5034 15.412 < 2e-16 ***
#> nuts16 70.6160 8.4813 8.326 < 2e-16 ***
#> nuts17 13.0092 12.1986 1.066 0.28622
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for gaussian family taken to be 955907.9)
#>
#> Null deviance: 3.0241e+11 on 216725 degrees of freedom
#> Residual deviance: 2.0710e+11 on 216656 degrees of freedom
#> AIC: 3599521
#>
#> Number of Fisher Scoring iterations: 2
(dre_estimator <- ess %>%
filter(tipojor == "1", estu == "7", sexo == "1") %>%
dr_estimator(prop_model, mean_model) |>
pull(dre))
#> [1] 153.843
(dre_estimator_madrid <- ess %>%
filter(tipojor == "1", estu == "7", sexo == "1", nuts1 == "3") %>%
dr_estimator(prop_model, mean_model) |>
pull(dre))
#> [1] 168.0452
(dre_estimator_sur <- ess %>%
filter(tipojor == "1", estu == "7", sexo == "1", nuts1 == "1") %>%
dr_estimator(prop_model, mean_model) |>
pull(dre))
#> [1] 183.358El S-learner usando ponderación de observaciones y el doubly robust estimator (sin usar ponderaciones) nos dan estimaciones diciendo que se gana más en el sector público que en el privado, mientras que el t-learner y el x-learner nos dicen lo contrario.
Así que, no me queda claro la respueta a la pregunta inicial.
res <- data.frame(s_learner_with_pond_madrid = s_learner_with_pond_madrid,
s_learner_with_pond_sur = s_learner_with_pond_sur,
t_learner_madrid= t_learner_madrid, t_learner_sur = t_learner_sur,
x_learner_madrid = x_learner_madrid, x_learner_sur = x_learner_sur,
dre_estimator_madrid = dre_estimator_madrid, dre_estimator_sur = dre_estimator_sur )
res |>
pivot_longer(everything(), names_to = "estimador", values_to = "valor")
#> # A tibble: 8 × 2
#> estimador valor
#> <chr> <dbl>
#> 1 s_learner_with_pond_madrid 67.9
#> 2 s_learner_with_pond_sur 140.
#> 3 t_learner_madrid -296.
#> 4 t_learner_sur -159.
#> 5 x_learner_madrid -501.
#> 6 x_learner_sur -55.7
#> 7 dre_estimator_madrid 168.
#> 8 dre_estimator_sur 183.Me acabo de acordar de otra forma de estimar esto. Consiste en:
Hago un modelo para estimar el salario neto pero sólo usando la población que trabaja en el sector privado.
Aplico ese modelo para obtener estimaciones sobre la población que trabaja en el sector público.
Comparo la estimación obtenida con el salario neto de esa subpoblación.
Es parecido al X-learner, pero sin tanta complicación. Es como decir ¿cuánto ganarían los que están en el sector público si estuvieran en el privado?
Al hacerlo así hay que obviar en el modelo en la subpoblación para el sector privado, las variables de cnace y de cno1 puesto que tienen niveles en el sector público que no están en el privado y el modelo daría error por niveles nuevos. Un modelo mixto si podría hacer eso.
mod_sector_privado <- svyglm(outcome ~ sexo + anos2 + estu + estrato2 + tipojor + anoanti + mesanti + tipocon + nuts1, design = disenno, subset = treatment == 0)
ess_sub_publico_madrid <- ess |>
filter(treatment == 1,tipojor == "1", estu == "7", sexo == "1", nuts1 == "3" )
estim_publico_con_mod_privado_madrid <- predict(mod_sector_privado,
ess_sub_publico_madrid)
(media_estimada_sector_publico_madrid <- weighted.mean(estim_publico_con_mod_privado_madrid, ess_sub_publico_madrid$factotal) )
#> [1] 2412.159
(media_observada_sector_publico <- weighted.mean(ess_sub_publico_madrid$outcome ,
ess_sub_publico_madrid$factotal) )
#> [1] 2535.845
(diferencia <- media_observada_sector_publico - media_estimada_sector_publico_madrid)
#> [1] 123.6857Y haciéndolo así se tendría que si los trabajadores del sector público hombres a jornada completa y con estudios de licenciados o superiores lse cambiaran al privado, manteniendo el resto igual ganarían unos 124 euros menos al mes de media.
Así que tengo varias preguntas
¿Qué piensan mis escasos lectores a la vista de estas estimaciones? ¿Hay brecha salarial?
Sea cual sea la respuesta, ¿no os parece que podría utilizar un método u otro según lo que me interese contar? Dan ganas de escribir un manual sobre como “engañar con estadística de forma avanzada”, pero ya conozco a quien tiene esa idea en mente
¿Cuál sería la metodología correcta si es que existe? ¿Quizá adentrándonos en el mundo bayesiano? ¿o es todo un artificio “técnico”?
Una de las hipótesis es que en el sector privado la distribución es más dispersa.
Todo esto debería haberse hecho antes que todos los modelos. EL EDA es lo primero.
# En general parece que no
ess |>
mutate(nuts1 = as_factor(nuts1)) |>
ggplot(aes(x =outcome,fill = as_factor(treatment))) +
geom_density(alpha = 0.5) + scale_x_continuous(limits = c(0,8000)) +
facet_wrap( ~ nuts1)
Pero y si la vemos, para los licenciados, hombres y a jornada completa.
En la comunidad de Madrid y en Canarias si se aprecia que en la cola de la derecha es superior la función de densidad en el sector privado.
ess |>
filter(tipojor == "1", estu == "7", sexo == "1") |>
mutate(nuts1 = as_factor(nuts1)) |>
ggplot(aes(x =outcome,fill = as_factor(treatment))) +
geom_density(alpha = 0.4) + scale_x_continuous(limits = c(0,8000)) +
facet_wrap( ~ nuts1)
Y si lo vemos para un par de ocupaciones, tales como Directores y gerentes (A0) y para técnicos y profesionales científicos e intelectuales de la salud y la enseñanza(B0).
Gerentes (A0) mejor en el privado, curritos de la enseñanza y la salud pues..
ess |>
filter(tipojor == "1", estu == "7", sexo == "1", cno1 %in% c("A0", "B0")) |>
mutate(nuts1 = as_factor(nuts1)) |>
ggplot(aes(x = outcome, fill = as_factor(treatment))) +
geom_density(alpha = 0.4) + scale_x_continuous(limits = c(0, 10000)) +
facet_grid( nuts1 ~ cno1) 
Y si usamos el primer modelo que vimos
ess |>
mutate(estim = predict(mod_sector_privado, ess)) |>
filter(tipojor == "1", estu == "7", sexo == "1", cno1 %in% c("A0", "B0")) |>
mutate(nuts1 = as_factor(nuts1)) |>
ggplot(aes(x = estim, fill = as_factor(treatment))) +
geom_density(alpha = 0.4) + scale_x_continuous(limits = c(0, 10000)) +
facet_grid( nuts1 ~ cno1) 
Pero claramente como lo que predice es la media condicionada queda todo muy “centrado”. Lo suyo sería un modelo bayesiano o hacer boostraping y tener la posterior predictive, para incorporar correctamente la variabilidad. A ver si lo hago en otro post, pero con numpyro, que stan no puede con estos datos