complexr: Análisis de Datos de Encuestas Complejas

1.1 Población y muestra

Sea \(U = \{1, 2, \ldots, N\}\) la población finita de tamaño \(N\), y \(s \subset U\) la muestra seleccionada bajo un diseño probabilístico \(p(s)\). Para cada unidad \(k \in U\), se define \(y_k\) como el valor de la variable de interés. El total poblacional y la media poblacional se definen respectivamente como:

\[ Y = \sum_{k \in U} y_k, \qquad \bar{Y} = \frac{Y}{N}. \]

1.2 Pesos de diseño y pesos ajustados

La probabilidad de que la unidad \(k\) sea incluida en la muestra se denota \(\pi_k = \Pr(k \in s) > 0\). El peso básico de diseño es \(d_k = 1/\pi_k\). En la práctica, estos pesos se modifican para incorporar ajustes por no respuesta o calibración a totales conocidos, obteniéndose los pesos ajustados \(w_k\). En esta documentación, \(w_k\) denota los pesos finales disponibles en el archivo de microdatos (variable weight).

1.3 Estimador de Horvitz–Thompson

El estimador de Horvitz–Thompson (HT) del total poblacional es (Horvitz and Thompson 1952):

\[ \hat{Y}_{HT} = \sum_{k \in s} d_k\, y_k, \]

y el estimador del tamaño poblacional se define como:

\[ \hat{N}_{HT} = \sum_{k \in s} d_k. \]

Cuando se trabaja con pesos ajustados \(w_k\), el estimador HT ponderado toma la forma \(\hat{Y}_w = \sum_{k \in s} w_k\, y_k\).

1.4 Varianza bajo el enfoque de diseño

La varianza del estimador HT se estima como (Särndal et al. 1992):

\[ \hat{V}_p\!\left(\hat{Y}_{HT}\right) = \sum_{k \in s}\sum_{l \in s} \bigl(d_k d_l - d_{kl}\bigr)\, y_k\, y_l, \]

donde \(d_{kl} = 1/\pi_{kl}\) y \(\pi_{kl} = \Pr(k, l \in s)\) son las probabilidades de inclusión de segundo orden. En la práctica se emplean métodos equivalentes como la linealización de Taylor o la replicación (jackknife, bootstrap), que no requieren el cálculo explícito de \(\pi_{kl}\).

1.5 Diseño estratificado multietápico

Para un diseño con \(H\) estratos, \(\alpha_h\) unidades primarias de muestreo (UPM) en el estrato \(h\) y \(n_{h\alpha}\) observaciones en la UPM \(\alpha\), el estimador del total es:

\[ \hat{Y}_{HT} = \sum_{h=1}^{H}\sum_{\alpha=1}^{\alpha_h} \sum_{k=1}^{n_{h\alpha}} \omega_{h\alpha k}\, y_{h\alpha k}, \]

donde \(\omega_{h\alpha k}\) es el peso ajustado del individuo \(k\) en la UPM \(\alpha\) del estrato \(h\).

1.6 Efecto del diseño (DEFF)

Siguiendo a Kish (1965), el efecto del diseño se define como la razón entre la varianza del estimador bajo el diseño complejo y la varianza del mismo estimador bajo un muestreo aleatorio simple (MAS) de igual tamaño:

\[ \widehat{\text{DEFF}} = \frac{\hat{V}_p(\hat{\theta})}{\hat{V}_{\text{SRS}}(\hat{\theta})}. \]

Un valor \(\widehat{\text{DEFF}} > 1\) indica que el diseño complejo incrementa la varianza con respecto al MAS, mientras que \(\widehat{\text{DEFF}} < 1\) señala una ganancia de eficiencia, típica de diseños estratificados con estratos homogéneos.

3.1 Variables

3.2 Jerarquía y consistencia interna

La simulación garantiza las siguientes consistencias:

• Los pesos son constantes dentro del hogar: \(w_{h\alpha k} = w_{h\alpha}\) para todo \(k\) en el hogar.
• El ingreso per cápita es constante dentro del hogar: \(y_{h\alpha k} = Y_{h\alpha} / N_{h\alpha}\).
• Los individuos menores de 5 años tienen NA en educación.

El ingreso del hogar sigue un modelo gamma jerárquico con efectos aleatorios de UPM y hogar:

\[ Y_{h\alpha} \sim \text{Gamma}\!\left(\alpha_0,\; \beta_0 \cdot \exp(u_{h\alpha} + v_{h\alpha k})\right) \]

donde \(u_{h\alpha} \sim N(0, 0{,}09)\) es el efecto de UPM y \(v_{h\alpha k} \sim N(0, 0{,}04)\) es el efecto de hogar.

6.1 Construcción del objeto de diseño

as_survey_design_tbl() envuelve survey::svydesign() y retorna un objeto tbl_svy compatible con el ecosistema srvyr / survey (Lumley 2010).

design <- as_survey_design_tbl(
  data    = data,
  weight  = "weight",
  strata  = "strata",
  cluster = "upm",
  nest    = TRUE
)

class(design)
#> [1] "tbl_svy"        "survey.design2" "survey.design"

Configuraciones soportadas:

La función valida que los pesos \(w_k\) sean estrictamente positivos (\(w_k > 0\)) y sin valores faltantes, y opcionalmente verifica que las UPMs no se compartan entre estratos (check_psu = TRUE).

Nota: Cuando un estrato \(h\) contiene una única UPM (\(\alpha_h = 1\)), la estimación de \(\hat{V}_p\) por linealización de Taylor no está definida. La función establece automáticamente options(survey.lonely.psu = "adjust") para usar la aproximación conservadora centrada en la media del estrato (Cochran 1977).

6.2 Diagnóstico del diseño

describe_survey_design(design)
#> # A tibble: 1 × 7
#>   n_obs n_strata n_clusters weight_min weight_max weight_mean weight_cv
#>   <int>    <int>      <int>      <dbl>      <dbl>       <dbl>     <dbl>
#> 1  3438        5        100       205.      1500.        848.     0.444

La tabla de diagnóstico reporta:

7.1 Media poblacional

El total poblacional ponderado se estima mediante el estimador HT:

\[ \hat{Y}_w = \sum_{k \in s} w_k\, y_k, \qquad \hat{N}_w = \sum_{k \in s} w_k. \]

La media ponderada (estimador de razón de Horvitz–Thompson) es:

\[ \bar{y}_w = \frac{\hat{Y}_w}{\hat{N}_w} = \frac{\displaystyle\sum_{k \in s} w_k\, y_k}{\displaystyle\sum_{k \in s} w_k}. \]

La varianza de \(\bar{y}_w\) se estima por linealización de Taylor (Särndal et al. 1992):

\[ \hat{V}_p\!\left(\bar{y}_w\right) = \frac{1}{\hat{N}_w^2}\, \hat{V}_p\!\left(\hat{Y}_w\right). \]

r_media <- estimate_survey(
  design    = design,
  variable  = "ingreso_pc",
  estimator = "mean"
)
r_media
#> # A tibble: 1 × 9
#>   variable   estimator estimate    se     cv  deff   lci   uci quality          
#>   <chr>      <chr>        <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <chr>            
#> 1 ingreso_pc mean         1579.  71.6 0.0453  8.05 1439. 1720. Very high precis…

7.2 Total poblacional

El estimador HT del total, bajo un diseño estratificado con UPMs, es:

\[ \hat{Y}_w = \sum_{h=1}^{H}\sum_{\alpha=1}^{\alpha_h}\sum_{k=1}^{n_{h\alpha}} \omega_{h\alpha k}\, y_{h\alpha k}. \]

Su varianza se estima por estratos y conglomerados:

\[ \hat{V}_p\!\left(\hat{Y}_w\right) = \sum_{h=1}^{H} \hat{V}_{p,h}\!\left(\hat{Y}_{w,h}\right), \]

donde \(\hat{V}_{p,h}\) se calcula dentro de cada estrato \(h\) usando los residuos de los totales de UPM respecto a su media en el estrato.

r_total <- estimate_survey(
  design    = design,
  variable  = "ingreso_pc",
  estimator = "total"
)
r_total
#> # A tibble: 1 × 9
#>   variable   estimator    estimate         se     cv  deff    lci    uci quality
#>   <chr>      <chr>           <dbl>      <dbl>  <dbl> <dbl>  <dbl>  <dbl> <chr>  
#> 1 ingreso_pc total     4607014476. 261136568. 0.0567  12.6 4.10e9 5.12e9 High p…

7.3 Proporciones

r_pobre <- estimate_survey(
  design    = design,
  variable  = "pobre",
  estimator = "prop"
)
r_pobre
#> # A tibble: 1 × 9
#>   variable estimator estimate      se     cv  deff   lci   uci quality          
#>   <chr>    <chr>        <dbl>   <dbl>  <dbl> <dbl> <dbl> <dbl> <chr>            
#> 1 pobre    prop         0.259 0.00936 0.0361  1.57 0.241 0.277 Very high precis…

r_empleo <- estimate_survey(
  design    = design,
  variable  = "empleo",
  estimator = "prop"
)
r_empleo
#> # A tibble: 3 × 10
#>   variable estimator estimate      se     cv  deff empleo      lci   uci quality
#>   <chr>    <chr>        <dbl>   <dbl>  <dbl> <dbl> <chr>     <dbl> <dbl> <chr>  
#> 1 empleo   prop         0.331 0.00815 0.0246  1.03 Formal    0.315 0.347 Very h…
#> 2 empleo   prop         0.341 0.00903 0.0265  1.25 Informal  0.323 0.359 Very h…
#> 3 empleo   prop         0.328 0.00866 0.0264  1.17 Unemploy… 0.311 0.345 Very h…

7.4 Estimador de razón

El estimador de razón de dos totales poblacionales es (Cochran 1977):

\[ \hat{R} = \frac{\hat{Y}_w}{\hat{X}_w} = \frac{\displaystyle\sum_{k \in s} w_k\, y_k} {\displaystyle\sum_{k \in s} w_k\, x_k}. \]

La varianza se estima por linealización de Taylor de primer orden:

\[ \hat{V}_p(\hat{R}) \approx \frac{1}{\hat{X}_w^2}\, \hat{V}_p\!\left(\hat{Y}_w - \hat{R}\,\hat{X}_w\right). \]

Numérico / Numérico:

r_razon <- estimate_survey(
  design      = design,
  estimator   = "ratio",
  numerator   = "ingreso_pc",
  denominator = "gasto_pc"
)
r_razon
#> # A tibble: 1 × 9
#>   variable          estimator estimate      se      cv  deff   lci   uci quality
#>   <chr>             <chr>        <dbl>   <dbl>   <dbl> <dbl> <dbl> <dbl> <chr>  
#> 1 ingreso_pc_over_… ratio         1.34 0.00579 0.00432    NA  1.33  1.35 Very h…

Categórico / Categórico — razón de trabajadores Formales respecto a Informales (\(\hat{N}_{\text{Formal}} / \hat{N}_{\text{Informal}}\)):

r_razon_cat <- estimate_survey(
  design          = design,
  estimator       = "ratio",
  numerator       = "empleo",
  denominator     = "empleo",
  ratio_num_level = "Formal",
  ratio_den_level = "Informal"
)
r_razon_cat
#> # A tibble: 1 × 9
#>   variable           estimator estimate     se     cv  deff   lci   uci quality 
#>   <chr>              <chr>        <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl> <chr>   
#> 1 empleo_over_empleo ratio        0.970 0.0429 0.0442    NA 0.886  1.05 Very hi…

Numérico / Categórico — ingreso promedio entre trabajadores formales, equivalente a \(\hat{Y}_{\text{ingreso}} / \hat{N}_{\text{Formal}}\):

r_razon_mix <- estimate_survey(
  design          = design,
  estimator       = "ratio",
  numerator       = "ingreso_pc",
  denominator     = "empleo",
  ratio_den_level = "Formal"
)
r_razon_mix
#> # A tibble: 1 × 9
#>   variable             estimator estimate    se     cv  deff   lci   uci quality
#>   <chr>                <chr>        <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <chr>  
#> 1 ingreso_pc_over_emp… ratio        4773.  242. 0.0507    NA 4299. 5247. High p…

7.5 Cuantiles

Los cuantiles se derivan de la función de distribución acumulada empírica ponderada (Woodruff 1952):

\[ \hat{F}_w(t) = \frac{\displaystyle\sum_{k \in s} w_k\, I(y_k \le t)} {\displaystyle\sum_{k \in s} w_k} = \frac{\hat{N}(y \le t)}{\hat{N}_w}. \]

El cuantil de orden \(p \in (0,1)\) se define como:

\[ \hat{q}_p = \inf\bigl\{t : \hat{F}_w(t) \ge p\bigr\}. \]

Los intervalos de confianza se calculan con el método de linealización de Woodruff, que transforma el problema a la escala de la proporción acumulada:

\[ IC_p[\hat{q}_p] = \left\{t : \hat{F}_w(t) \in \left[p \pm t_{1-\alpha/2,\,gl}\; ee(\hat{F}_w(t))\right]\right\}. \]

r_cuant <- estimate_survey(
  design    = design,
  variable  = "ingreso_pc",
  estimator = "quantile",
  probs     = c(0.10, 0.25, 0.50, 0.75, 0.90)
)
r_cuant
#> # A tibble: 5 × 10
#>   variable   estimator quantile estimate    se     cv  deff   lci   uci quality 
#>   <chr>      <chr>        <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <chr>   
#> 1 ingreso_pc quantile      0.1      499.  28.0 0.0562    NA  444.  554. High pr…
#> 2 ingreso_pc quantile      0.25     736.  32.3 0.0438    NA  673.  800. Very hi…
#> 3 ingreso_pc quantile      0.5     1127.  42.0 0.0373    NA 1045. 1209. Very hi…
#> 4 ingreso_pc quantile      0.75    1895. 106.  0.0561    NA 1686. 2103. High pr…
#> 5 ingreso_pc quantile      0.9     3196. 221.  0.0691    NA 2763. 3629. High pr…

8.1 Efecto de diseño de dominio

El efecto de diseño para el dominio \(d\) se define de forma análoga al DEFF global (Kish 1965):

\[ \widehat{\text{DEFF}}_d = \frac{\hat{V}_p(\hat{\theta}_d)} {\hat{V}_{\text{SRS}}(\hat{\theta}_d)}, \]

donde \(\hat{V}_{\text{SRS}}(\hat{\theta}_d)\) es la varianza que se obtendría bajo un MAS restringido al dominio \(d\). Para la media de dominio, esta expresión se simplifica a:

\[ \hat{V}_{\text{SRS}}(\bar{y}_{w,d}) = \left(1 - \frac{n_d}{N_d}\right)\frac{S_{y,d}^2}{n_d}, \]

con \(n_d = \sum_{k \in s} I(k \in U_d)\) el tamaño muestral en el dominio, \(N_d \approx \hat{N}_{w,d}\) el tamaño poblacional estimado del dominio, y \(S_{y,d}^2\) la varianza muestral no ponderada dentro del dominio. Valores \(\widehat{\text{DEFF}}_d > 1\) indican que la conglomeración o la ponderación desigual inflan la varianza incluso dentro del dominio.

El argumento by controla los dominios de estimación:

r_region <- estimate_survey(
  design    = design,
  variable  = "ingreso_pc",
  estimator = "mean",
  by        = "region"
)
r_region
#> # A tibble: 3 × 10
#>   variable   estimator region estimate    se     cv  deff   lci   uci quality   
#>   <chr>      <chr>     <chr>     <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <chr>     
#> 1 ingreso_pc mean      Center    1597.  88.6 0.0555  4.03 1423. 1770. High prec…
#> 2 ingreso_pc mean      North     1616.  85.2 0.0527  3.62 1449. 1783. High prec…
#> 3 ingreso_pc mean      South     1523.  71.1 0.0467  2.89 1384. 1662. Very high…

Se admiten múltiples variables de dominio (dominios cruzados \(U_{d_1} \cap U_{d_2}\)):

r_region_area <- estimate_survey(
  design    = design,
  variable  = "ingreso_pc",
  estimator = "mean",
  by        = c("region", "area")
)
r_region_area
#> # A tibble: 6 × 11
#> # Groups:   region [3]
#>   variable   estimator region area  estimate    se     cv  deff   lci   uci
#>   <chr>      <chr>     <chr>  <chr>    <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
#> 1 ingreso_pc mean      Center Rural    1642. 109.  0.0663  2.91 1429. 1856.
#> 2 ingreso_pc mean      Center Urban    1541. 110.  0.0713  3.44 1325. 1756.
#> 3 ingreso_pc mean      North  Rural    1687. 115.  0.0679  2.76 1462. 1911.
#> 4 ingreso_pc mean      North  Urban    1545.  90.9 0.0588  2.57 1367. 1723.
#> 5 ingreso_pc mean      South  Rural    1512.  91.6 0.0606  2.39 1332. 1691.
#> 6 ingreso_pc mean      South  Urban    1536.  89.0 0.0579  2.29 1362. 1711.
#> # ℹ 1 more variable: quality <chr>

Proporciones por dominio:

r_pobre_region <- estimate_survey(
  design    = design,
  variable  = "pobre",
  estimator = "prop",
  by        = "region"
)
r_pobre_region
#> # A tibble: 3 × 10
#>   variable estimator region estimate     se     cv  deff   lci   uci quality    
#>   <chr>    <chr>     <chr>     <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl> <chr>      
#> 1 pobre    prop      Center    0.250 0.0155 0.0621  1.42 0.220 0.281 High preci…
#> 2 pobre    prop      North     0.252 0.0132 0.0523  1.11 0.227 0.278 High preci…
#> 3 pobre    prop      South     0.275 0.0148 0.0537  1.24 0.246 0.304 High preci…