6.4 Más alla de muestras aleatorias simples

Introdujimos el bootstrap en el contexto de muestras aleatorias, esto es, suponiendo que las observaciones son independientes; en este escenario basta con aproximar la distribución desconocida \(P\) usando la dsitribución empírica \(P_n\), y el cálculo de los estadísticos es inmediato. Hay casos en los que el mecanismo que generó los datos es más complicado, por ejemplo, cuando tenemos dos muestras, en diseños de encuestas complejas o en series de tiempo.

Ejemplo: Dos muestras

En el ejemplo de experimentos clínicos de aspirina y ataques de de corazón, podemos pensar el modelo probabilístico \(P\) como compuesto por dos distribuciones de probabilidad \(G\) y \(Q\) una correspondiente al grupo control y otra al grupo de tratamiento, entonces las observaciones de cada grupo provienen de distribuciones distintas y el método bootstrap debe tomar en cuenta esto al generar las muestras, en este caso implica seleccionar muesreas de manera independiente dentro de cada grupo.

Ejemplo: Bootstrap en muestreo de encuestas

La necesidad de estimaciones confiables junto con el uso eficiente de recursos conllevan a diseños de muestras complejas. Estos diseños típicamente usan las siguientes técnicas: muestreo sin reemplazo de una población finita, muestreo sistemático, estratificación, conglomerados, ajustes a no-respuesta, postestratificación. Como consecuencia, los valores de la muestra suelen no ser independientes.

La complejidad de los diseños de encuestas conlleva a que el cálculo de errores estándar sea muy complicado, para atacar este problema hay dos técnicas básicas: 1) un enfoque analítico usando linearización, 2) métodos de remuestreo como bootstrap. El incremento en el poder de cómputo ha favorecido los métodos de remuestreo pues la linearización requiere del desarrollo de una fórmula para cada estimación y supuestos adicionales para simplificar.

En 1988 Rao and Wu (1988) propusieron un método de bootstrap para diseños estratificados multietápicos con reemplazo de UPMs que describimos a continuación.

ENIGH. Usaremos como ejemplo la Encuesta Nacional de Ingresos y Gastos de los Hogares, ENIGH 2018 (INEGI 2018), esta encuesta usa un diseño de conglomerados estratificado.

Antes de proceder a bootstrap debemos entender como se seleccionaron los datos, esto es, el diseño de la muestra:

Unidad primaria de muestreo (UPM). Las UPMs están constituidas por agrupaciones de viviendas. Se les denomina unidades primarias pues corresponden a la primera etapa de selección, las unidades secundarias (USMs) serían los hogares.
Estratificación. Los estratos se construyen en base a estado, ámbito (urbano, complemento urbano, rural), características sociodemográficas de los habitantes de las viviendas, características físicas y equipamiento. El proceso de estratificación resulta en 888 subestratos en todo el ámbito nacional.
La selección de la muestra es independiente para cada estrato, y una vez que se obtiene la muestra se calculan los factores de expansión que reflejan las distintas probabilidades de selección. Después se llevan a cabo ajustes por no respuesta y por proyección (calibración), esta última busca que distintos dominios de la muestra coincidan con la proyección de población de INEGI.

library(usethis)
use_zip("https://www.inegi.org.mx/contenidos/programas/enigh/nc/2018/datosabiertos/conjunto_de_datos_enigh_2018_ns_csv.zip", "data")

library(here)

concentrado_hogar <- read_csv(here("data", 
    "conjunto_de_datos_enigh_2018_ns_csv", 
    "conjunto_de_datos_concentradohogar_enigh_2018_ns", "conjunto_de_datos",
    "conjunto_de_datos_concentradohogar_enigh_2018_ns.csv"))
glimpse(concentrado_hogar)
#> Observations: 74,647
#> Variables: 126
#> $ folioviv   <dbl> 100013601, 100013602, 100013603, 100013604, 100013606, 100…
#> $ foliohog   <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
#> $ ubica_geo  <dbl> 1001, 1001, 1001, 1001, 1001, 1001, 1001, 1001, 1001, 1001…
#> $ tam_loc    <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
#> $ est_socio  <dbl> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3…
#> $ est_dis    <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
#> $ upm        <dbl> 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5…
#> $ factor     <dbl> 175, 175, 175, 175, 175, 189, 189, 189, 189, 186, 186, 186…
#> $ clase_hog  <dbl> 2, 2, 2, 2, 2, 2, 1, 2, 2, 3, 2, 1, 2, 2, 2, 2, 3, 1, 3, 2…
#> $ sexo_jefe  <dbl> 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1…
#> $ edad_jefe  <dbl> 74, 48, 39, 70, 51, 41, 57, 53, 30, 69, 76, 77, 70, 29, 68…
#> $ educa_jefe <dbl> 4, 11, 10, 8, 4, 11, 9, 11, 6, 4, 3, 4, 6, 6, 9, 7, 6, 8, …
#> $ tot_integ  <dbl> 3, 5, 2, 2, 4, 4, 1, 2, 3, 4, 2, 1, 2, 4, 4, 2, 5, 1, 4, 3…
#> $ hombres    <dbl> 2, 2, 1, 1, 1, 2, 0, 1, 2, 4, 0, 1, 1, 2, 3, 1, 2, 1, 2, 1…
#> $ mujeres    <dbl> 1, 3, 1, 1, 3, 2, 1, 1, 1, 0, 2, 0, 1, 2, 1, 1, 3, 0, 2, 2…
#> $ mayores    <dbl> 3, 5, 2, 2, 3, 4, 1, 2, 2, 3, 2, 1, 2, 2, 4, 2, 5, 1, 3, 3…
#> $ menores    <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 2, 0, 0, 0, 0, 1, 0…
#> $ p12_64     <dbl> 1, 5, 2, 1, 3, 4, 1, 2, 2, 2, 1, 0, 0, 2, 2, 0, 5, 1, 3, 2…
#> $ p65mas     <dbl> 2, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 2, 0, 2, 2, 0, 0, 0, 1…
#> $ ocupados   <dbl> 2, 2, 2, 0, 2, 2, 1, 2, 2, 2, 1, 1, 1, 1, 2, 0, 4, 1, 2, 2…
#> $ percep_ing <dbl> 3, 5, 2, 2, 2, 2, 1, 2, 2, 4, 2, 1, 2, 1, 4, 2, 4, 1, 2, 3…
#> $ perc_ocupa <dbl> 2, 2, 2, 0, 2, 2, 1, 2, 2, 2, 1, 1, 1, 1, 2, 0, 4, 1, 2, 2…
#> $ ing_cor    <dbl> 76403.70, 42987.73, 580697.74, 46252.71, 53837.09, 237742.…
#> $ ingtrab    <dbl> 53114.74, 15235.06, 141885.21, 0.00, 43229.49, 129836.03, …
#> $ trabajo    <dbl> 53114.74, 0.00, 141885.21, 0.00, 8852.45, 129836.03, 23606…
#> $ sueldos    <dbl> 53114.74, 0.00, 133770.48, 0.00, 8852.45, 95901.63, 23606.…
#> $ horas_extr <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00…
#> $ comisiones <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 22131.14, 0.00, 0.00, 0.00, …
#> $ aguinaldo  <dbl> 0.00, 0.00, 3934.42, 0.00, 0.00, 11803.26, 0.00, 22131.14,…
#> $ indemtrab  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ otra_rem   <dbl> 0.00, 0.00, 4180.31, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0…
#> $ remu_espec <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00…
#> $ negocio    <dbl> 0.00, 13759.66, 0.00, 0.00, 34377.04, 0.00, 0.00, 0.00, 0.…
#> $ noagrop    <dbl> 0.00, 13759.66, 0.00, 0.00, 34377.04, 0.00, 0.00, 0.00, 0.…
#> $ industria  <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00…
#> $ comercio   <dbl> 0.00, 0.00, 0.00, 0.00, 34377.04, 0.00, 0.00, 0.00, 0.00, …
#> $ servicios  <dbl> 0.00, 13759.66, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, …
#> $ agrope     <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ agricolas  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ pecuarios  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ reproducc  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ pesca      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ otros_trab <dbl> 0.0, 1475.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, …
#> $ rentas     <dbl> 0.00, 0.00, 29508.19, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, …
#> $ utilidad   <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00…
#> $ arrenda    <dbl> 0.00, 0.00, 29508.19, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, …
#> $ transfer   <dbl> 11288.96, 3752.67, 391304.34, 34252.71, 107.60, 89906.51, …
#> $ jubilacion <dbl> 9147.54, 0.00, 0.00, 23606.55, 0.00, 23606.55, 0.00, 0.00,…
#> $ becas      <dbl> 0.0, 491.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0…
#> $ donativos  <dbl> 0.00, 147.54, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 23…
#> $ remesas    <dbl> 0.00, 98.36, 0.00, 5901.63, 0.00, 0.00, 0.00, 0.00, 0.00, …
#> $ bene_gob   <dbl> 1622.95, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0…
#> $ transf_hog <dbl> 0.00, 3014.97, 0.00, 0.00, 107.60, 61714.26, 0.00, 0.00, 0…
#> $ trans_inst <dbl> 518.47, 0.00, 391304.34, 4744.53, 0.00, 4585.70, 0.00, 0.0…
#> $ estim_alqu <dbl> 12000.00, 24000.00, 18000.00, 12000.00, 10500.00, 18000.00…
#> $ otros_ing  <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00…
#> $ gasto_mon  <dbl> 18551.47, 55470.99, 103106.89, 19340.06, 13605.03, 33627.7…
#> $ alimentos  <dbl> 5618.47, 20930.29, 37594.06, 2892.84, 7367.09, 0.00, 11455…
#> $ ali_dentro <dbl> 4075.63, 8587.46, 25251.25, 2892.84, 4795.67, 0.00, 8344.1…
#> $ cereales   <dbl> 964.25, 2689.65, 3728.53, 385.71, 257.14, 0.00, 437.13, 10…
#> $ carnes     <dbl> 0.00, 1401.41, 2828.56, 2121.42, 2931.41, 0.00, 1787.13, 3…
#> $ pescado    <dbl> 745.71, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.…
#> $ leche      <dbl> 0.00, 443.55, 4345.70, 0.00, 0.00, 0.00, 2841.41, 1491.40,…
#> $ huevo      <dbl> 719.98, 0.00, 411.42, 0.00, 0.00, 0.00, 308.57, 629.99, 30…
#> $ aceites    <dbl> 0.00, 0.00, 1928.57, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0…
#> $ tuberculo  <dbl> 0.00, 257.14, 385.71, 0.00, 128.57, 0.00, 231.42, 411.42, …
#> $ verduras   <dbl> 745.70, 1893.29, 2635.66, 0.00, 835.70, 0.00, 861.38, 1028…
#> $ frutas     <dbl> 0.00, 533.16, 1864.27, 0.00, 0.00, 0.00, 244.27, 809.98, 0…
#> $ azucar     <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 257.14, 0.00, 0.00, 0.…
#> $ cafe       <dbl> 0.00, 462.85, 1414.28, 0.00, 0.00, 0.00, 964.28, 0.00, 0.0…
#> $ especias   <dbl> 0.00, 167.14, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.…
#> $ otros_alim <dbl> 0.00, 392.13, 2545.71, 385.71, 514.28, 0.00, 0.00, 2699.99…
#> $ bebidas    <dbl> 899.99, 347.14, 3162.84, 0.00, 128.57, 0.00, 411.42, 385.7…
#> $ ali_fuera  <dbl> 771.42, 12342.83, 12342.81, 0.00, 2571.42, 0.00, 1928.57, …
#> $ tabaco     <dbl> 771.42, 0.00, 0.00, 0.00, 0.00, 0.00, 1182.84, 0.00, 0.00,…
#> $ vesti_calz <dbl> 0.00, 401.06, 2015.21, 97.82, 0.00, 0.00, 0.00, 1565.20, 5…
#> $ vestido    <dbl> 0.00, 224.98, 2015.21, 97.82, 0.00, 0.00, 0.00, 293.47, 53…
#> $ calzado    <dbl> 0.00, 176.08, 0.00, 0.00, 0.00, 0.00, 0.00, 1271.73, 0.00,…
#> $ vivienda   <dbl> 3912.00, 2495.00, 4475.00, 1458.00, 300.00, 2801.00, 4405.…
#> $ alquiler   <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 3900.00, 1…
#> $ pred_cons  <dbl> 0.00, 1250.00, 1250.00, 0.00, 0.00, 140.00, 250.00, 0.00, …
#> $ agua       <dbl> 312.00, 750.00, 750.00, 600.00, 0.00, 741.00, 630.00, 690.…
#> $ energia    <dbl> 3600.00, 495.00, 2475.00, 858.00, 300.00, 1920.00, 3525.00…
#> $ limpieza   <dbl> 522.00, 412.16, 3318.26, 5514.00, 3300.00, 5682.00, 2496.6…
#> $ cuidados   <dbl> 522.00, 375.00, 2340.00, 5514.00, 3300.00, 5682.00, 2301.0…
#> $ utensilios <dbl> 0.00, 37.16, 978.26, 0.00, 0.00, 0.00, 195.65, 391.30, 0.0…
#> $ enseres    <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5901.63, 0.00, 0…
#> $ salud      <dbl> 0.00, 1348.99, 28858.68, 322.82, 56.73, 0.00, 4695.64, 0.0…
#> $ atenc_ambu <dbl> 0.00, 1007.59, 28858.68, 0.00, 56.73, 0.00, 4695.64, 0.00,…
#> $ hospital   <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 8217…
#> $ medicinas  <dbl> 0.00, 341.40, 0.00, 322.82, 0.00, 0.00, 0.00, 0.00, 0.00, …
#> $ transporte <dbl> 8400.00, 7628.56, 12325.68, 7350.00, 600.00, 18235.70, 580…
#> $ publico    <dbl> 0.00, 578.56, 4255.68, 0.00, 0.00, 1285.70, 0.00, 0.00, 0.…
#> $ foraneo    <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 590.16, 0.00, 0.00, 70…
#> $ adqui_vehi <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00…
#> $ mantenim   <dbl> 7200.00, 3600.00, 4500.00, 6000.00, 0.00, 13200.00, 4253.1…
#> $ refaccion  <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 2213.11, 0.00, 0.00, 0…
#> $ combus     <dbl> 7200.00, 3600.00, 4500.00, 6000.00, 0.00, 13200.00, 2040.0…
#> $ comunica   <dbl> 1200.00, 3450.00, 3570.00, 1350.00, 600.00, 3750.00, 960.0…
#> $ educa_espa <dbl> 0.00, 17567.05, 0.00, 639.34, 0.00, 1800.00, 627.00, 3600.…
#> $ educacion  <dbl> 0.00, 8547.39, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0…
#> $ esparci    <dbl> 0.00, 167.21, 0.00, 639.34, 0.00, 1800.00, 627.00, 3600.00…
#> $ paq_turist <dbl> 0.00, 8852.45, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0…
#> $ personales <dbl> 99.00, 4663.29, 8520.00, 1065.24, 1686.13, 5109.00, 3335.6…
#> $ cuida_pers <dbl> 99.00, 1497.00, 8520.00, 180.00, 1647.00, 4509.00, 1500.00…
#> $ acces_pers <dbl> 0.00, 166.29, 0.00, 0.00, 39.13, 0.00, 0.00, 0.00, 0.00, 0…
#> $ otros_gas  <dbl> 0.00, 3000.00, 0.00, 885.24, 0.00, 600.00, 1835.65, 0.00, …
#> $ transf_gas <dbl> 0.00, 24.59, 6000.00, 0.00, 295.08, 0.00, 491.80, 23606.55…
#> $ percep_tot <dbl> 0.00, 6073.09, 3857.14, 1380.55, 0.00, 1928.57, 489.13, 20…
#> $ retiro_inv <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ prestamos  <dbl> 0.00, 7.37, 0.00, 737.70, 0.00, 0.00, 0.00, 0.00, 491.80, …
#> $ otras_perc <dbl> 0.00, 462.28, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.…
#> $ ero_nm_viv <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ ero_nm_hog <dbl> 0.00, 5603.44, 3857.14, 642.85, 0.00, 1928.57, 489.13, 205…
#> $ erogac_tot <dbl> 0.00, 9009.82, 81147.53, 0.00, 0.00, 14754.09, 0.00, 58229…
#> $ cuota_viv  <dbl> 0, 0, 0, 0, 0, 0, 0, 12000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
#> $ mater_serv <dbl> 0.00, 147.54, 0.00, 0.00, 0.00, 0.00, 0.00, 7868.85, 0.00,…
#> $ material   <dbl> 0.00, 147.54, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.…
#> $ servicio   <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 7868.85, 0.00, 0…
#> $ deposito   <dbl> 0.00, 9.83, 66393.44, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, …
#> $ prest_terc <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00…
#> $ pago_tarje <dbl> 0.00, 8852.45, 0.00, 0.00, 0.00, 14754.09, 0.00, 0.00, 0.0…
#> $ deudas     <dbl> 0.00, 0.00, 14754.09, 0.00, 0.00, 0.00, 0.00, 38360.65, 0.…
#> $ balance    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ otras_erog <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00…
#> $ smg        <dbl> 7952.4, 7952.4, 7952.4, 7952.4, 7952.4, 7952.4, 7952.4, 79…

# seleccionar variable de ingreso corriente
hogar <- concentrado_hogar %>% 
    mutate(
        upm = as.integer(upm),
        jefe_hombre = sexo_jefe == 1, 
        edo = str_sub(ubica_geo, 1, 2), 
        jefa_50 = (sexo_jefe == 2) & (edad_jefe > 50)
        ) %>% 
    select(folioviv, foliohog, est_dis, upm, factor, ing_cor, sexo_jefe, 
       edad_jefe, edo, jefa_50) %>% 
    group_by(est_dis) %>% 
    mutate(n = n_distinct(upm)) %>% # número de upms por estrato
    ungroup()
hogar
#> # A tibble: 74,647 x 11
#>    folioviv foliohog est_dis   upm factor ing_cor sexo_jefe edad_jefe edo  
#>       <dbl>    <dbl>   <dbl> <int>  <dbl>   <dbl>     <dbl>     <dbl> <chr>
#>  1   1.00e8        1       2     1    175  76404.         1        74 10   
#>  2   1.00e8        1       2     1    175  42988.         1        48 10   
#>  3   1.00e8        1       2     1    175 580698.         1        39 10   
#>  4   1.00e8        1       2     1    175  46253.         2        70 10   
#>  5   1.00e8        1       2     1    175  53837.         2        51 10   
#>  6   1.00e8        1       2     2    189 237743.         2        41 10   
#>  7   1.00e8        1       2     2    189  32607.         2        57 10   
#>  8   1.00e8        1       2     2    189 169918.         1        53 10   
#>  9   1.00e8        1       2     2    189  17311.         1        30 10   
#> 10   1.00e8        1       2     3    186 120488.         1        69 10   
#> # … with 74,637 more rows, and 2 more variables: jefa_50 <lgl>, n <int>

Para el cálculo de estadísticos debemos usar los factores de expansión, por ejemplo el ingreso trimestral total sería:

sum(hogar$factor * hogar$ing_cor / 1000)
#> [1] 1723700566

y ingreso trimestral medio (miles pesos)

sum(hogar$factor * hogar$ing_cor / 1000) / sum(hogar$factor)
#> [1] 49.61029

La estimación del error estándar, por otro lado, no es sencilla y requiere usar aproximaciones, en la metodología de INEGI proponen una aproximación con series de Taylor.

Figure 6.1: Extracto de estimación de errores de muestreo, ENIGH 2018.

Veamos ahora como calcular el error estándar siguiendo el bootstrap de Rao y Wu:

En cada estrato se seleccionan con reemplazo \(m_h\) UPMs de las \(n_h\) de la muestra original. Denotamos por \(m_{hi}^*\) el número de veces que se seleccionó la UPM \(i\) en el estrato \(h\) (de tal manera que \(\sum m_{hi}^*=m_h\)). Creamos una replicación del ponderador correspondiente a la \(k\)-ésima unidad (USM) como:

\[d_k^*=d_k \bigg[\bigg(1-\sqrt{\frac{m_h}{n_h - 1}}\bigg) + \bigg(\sqrt{\frac{m_h}{n_h - 1}}\frac{n_h}{m_h}m_{h}^*\bigg)\bigg]\]

donde \(d_k\) es el inverso de la probabilidad de selección. Si \(m_h<(n_h -1)\) todos los pesos definidos de esta manera serán no negativos. Calculamos el peso final \(w_k^*\) aplicando a \(d_k^*\) los mismos ajustes que se hicieron a los ponderadores originales.

Calculamos el estadístico de interés \(\hat{\theta}\) usando los ponderadores \(w_k^*\) en lugar de los originales \(w_k\).
Repetimos los pasos 1 y 2 \(B\) veces para obtener \(\hat{\theta}^{*1},\hat{\theta}^{*2},...,\hat{\theta}^{*B}\).
Calculamos el error estándar como:

\[\hat{se}_B = \bigg\{\frac{\sum_{b=1}^B[\hat{\theta}^*(b)-\hat{\theta}^*(\cdot)]^2 }{B}\bigg\}^{1/2}\]

Podemos elegir cualquier valor de \(m_h \geq 1\), el más sencillo es elegir \(m_h=n_h-1\), en este caso: \[d_k^*=d_k \frac{n_h}{n_h-1}m_{hi}^*\] en este escenario las unidades que no se incluyen en la muestra tienen un valor de cero como ponderador. Si elegimos \(n_h \ne n_h-1\) las unidades que no están en la muestra tienen ponderador distinto a cero, si \(m_h=n_h\) el ponderador podría tomar valores negativos.

Implementemos el bootstrap de Rao y Wu a la ENIGH, usaremos \(m_h=n_h-1\)

# creamos una tabla con los estratos y upms
est_upm <- hogar %>% 
    distinct(est_dis, upm, n)

hogar_factor <- est_upm %>% 
    split(.$est_dis) %>% # dentro de cada estrato tomamos muestra (n_h-1)
    map_df(~sample_n(., size = first(.$n) - 1, replace = TRUE)) %>% 
    add_count(upm, name = "m_hi") %>% # calculamos m_hi*
    left_join(hogar, by = c("est_dis", "upm", "n")) %>% 
    mutate(factor_b = factor * m_hi * n / (n - 1))

# unimos los pasos anteriores en una función para replicar en cada muestra bootstrap
svy_boot <- function(est_upm, hogar){
    m_hi <- est_upm %>% 
        split(.$est_dis) %>% 
        map(~sample(.$upm, size = first(.$n) - 1, replace = TRUE)) %>% 
        flatten_int() %>% 
        plyr::count() %>% 
        select(upm = x, m_h = freq)
    m_hi %>% 
        left_join(hogar, by = c("upm")) %>% 
        mutate(factor_b = factor * m_h * n / (n - 1))
}
set.seed(1038984)
boot_rep <- rerun(500, svy_boot(est_upm, hogar))

# Aplicación a ingreso medio
wtd_mean <- function(w, x, na.rm = FALSE) {
    sum(w * x, na.rm = na.rm) / sum(w, na.rm = na.rm)
} 

# La media es:
hogar %>% 
    summarise(media = wtd_mean(factor, ing_cor))
#> # A tibble: 1 x 1
#>    media
#>    <dbl>
#> 1 49610.

Y el error estándar:

map_dbl(boot_rep, ~wtd_mean(w = .$factor_b, x = .$ing_cor)) %>% sd()
#> [1] 441.0439

El método bootstrap está implementado en el paquete survey y más recientemente en srvyr que es una versión tidy que utiliza las funciones en survey.

Podemos comparar nuestros resultados con la implementación en survey.

# 1. Definimos el diseño de la encuesta
library(survey)
library(srvyr)

enigh_design <- hogar %>% 
    as_survey_design(ids = upm, weights = factor, strata = est_dis)

# 2. Elegimos bootstrap como el método para el cálculo de errores estándar
set.seed(7398731)
enigh_boot <- enigh_design %>% 
    as_survey_rep(type = "subbootstrap", replicates = 500)

# 3. Así calculamos la media
enigh_boot %>% 
    srvyr::summarise(mean_ingcor = survey_mean(ing_cor))
#> # A tibble: 1 x 2
#>   mean_ingcor mean_ingcor_se
#>         <dbl>          <dbl>
#> 1      49610.           459.

enigh_boot %>% 
    group_by(edo) %>% 
    srvyr::summarise(mean_ingcor = survey_mean(ing_cor)) 
#> # A tibble: 30 x 3
#>    edo   mean_ingcor mean_ingcor_se
#>    <chr>       <dbl>          <dbl>
#>  1 10         50161.           942.
#>  2 11         46142.          1252.
#>  3 12         29334.          1067.
#>  4 13         38783.           933.
#>  5 14         60541.          1873.
#>  6 15         48013.          1245.
#>  7 16         42653.          1239.
#>  8 17         42973.          1675.
#>  9 18         48148.          1822.
#> 10 19         68959.          3625.
#> # … with 20 more rows

# cuantiles
svyquantile(~ing_cor, enigh_boot, quantiles = seq(0.1, 1, 0.1), 
    interval.type = "quantile")
#> Statistic:
#>         ing_cor
#> q0.1   13155.75
#> q0.2   18895.37
#> q0.3   24041.89
#> q0.4   29358.29
#> q0.5   35505.47
#> q0.6   42695.44
#> q0.7   52426.32
#> q0.8   66594.08
#> q0.9   94613.04
#> q1   4501830.28
#> SE:
#>          ing_cor
#> q0.1    114.2707
#> q0.2    110.1885
#> q0.3    130.8151
#> q0.4    152.8712
#> q0.5    199.3702
#> q0.6    241.1244
#> q0.7    339.4501
#> q0.8    479.4980
#> q0.9    908.6814
#> q1   384477.9727

Supongamos que queremos calcular la media para los hogares con jefe de familia mujer mayor a 50 años.

# Creamos datos con filter y repetimos lo de arriba
hogar_mujer <- filter(hogar, jefa_50)
est_upm_mujer <- hogar_mujer %>% 
    distinct(est_dis, upm, n)
# bootstrap
boot_rep_mujer <- rerun(500, svy_boot(est_upm_mujer, hogar_mujer))
# media y error estándar
hogar_mujer %>% 
    summarise(media = wtd_mean(factor, ing_cor))
#> # A tibble: 1 x 1
#>    media
#>    <dbl>
#> 1 44356.
# usamos bootstrap para calcular los errores estándar
map_dbl(boot_rep_mujer, ~wtd_mean(w = .$factor_b, x = .$ing_cor, na.rm = TRUE)) %>% 
    sd()
#> [1] 546.8034

Comparemos con los resultados de srvyr. ¿qué pasa?

library(srvyr)
enigh_boot %>% 
    srvyr::group_by(jefa_50) %>% 
    srvyr::summarise(mean_ingcor = survey_mean(ing_cor))
#> # A tibble: 2 x 3
#>   jefa_50 mean_ingcor mean_ingcor_se
#>   <lgl>         <dbl>          <dbl>
#> 1 FALSE        50574.           502.
#> 2 TRUE         44356.           726.

Sub-poblaciones como “jefas de familia mayores a 50” se conocen como un dominio, esto es un subgrupo cuyo tamaño de muestra es aleatorio, este ejemplo nos recalca la importancia de considerar el proceso en que se generó la muestra para calcular los errores estándar bootstrap.

map_dbl(boot_rep, 
    function(x){hm <- filter(x, jefa_50); 
    wtd_mean(w = hm$factor_b, x = hm$ing_cor)}) %>% 
    sd()
#> [1] 715.9535

Resumiendo:

El bootstrap de Rao y Wu genera un estimador consistente y aproximadamente insesgado de la varianza de estadísticos no lineales y para la varianza de un cuantil.
Este método supone que la seleccion de UPMs es con reemplazo; hay variaciones del estimador bootstrap de Rao y Wu que extienden el método que acabamos de estudiar; sin embargo, es común ignorar este aspecto, por ejemplo Mach et al estudian las propiedades del estimador de varianza bootstrap de Rao y Wu cuando la muestra se seleccionó sin reemplazo.

Referencias

INEGI. 2018. “Encuesta Nacional de Ingresos Y Gastos de Los Hogares (Enigh-2018). Diseño Muestral.” http://internet.contenidos.inegi.org.mx/contenidos/Productos/prod_serv/contenidos/espanol/bvinegi/productos/nueva_estruc/702825070359.pdf.

Rao, J. N. K., and C. F. J. Wu. 1988. “Resampling Inference with Complex Survey Data.” Journal of the American Statistical Association 83 (401). Taylor & Francis:231–41. https://doi.org/10.1080/01621459.1988.10478591.