8 Análisis de Encuestas (ENSANUT)

8.1 Paquetes y datos a utilizar

A lo largo de esta sección usaremos los siguientes paquetes:

library(haven)
library(wesanderson)
library(tidyverse)
library(srvyr)
library(kableExtra)
library(survey)
library(broom) # para tidy() sobre los modelos

Para los ejercicios usaremos al Encuesta Nacional de Salud y Nutrición 2021 (ENSANUT 2021) disponible en la siguiente liga: https://ensanut.insp.mx/encuestas/ensanutcontinua2021/descargas.php

En particular usaremos el cuestionario de salud de adultos Cuestionario de salud de adultos (20 años o más) en formato csv así como el Cuestionario de antropometría y tensión arterial igual en formato csv. Acá leeré ambos como salud_adultos_svy y antropometria:

salud_adultos <- read_csv("datasets/ensadul2021_entrega_w_15_12_2021.csv")
antropometria <- read_csv("datasets/ensaantro21_entrega_w_17_12_2021.csv")

Finalmente, a modo de práctica y explicación utilizaremos la siguiente base de datos generada por el código de abajo. La base consiste en \(12\) individuos con dos grupos de edad (Menor a 30 y Mayor o igual a 30) que viven en distintos hogares numerados del 1 al 5 en la columna (hogar). Los hogares se encuentran en diferentes localidades (del 1 al 3) y tienen cierta estatura (estatura) y sexo (sexo).

# Código para generar un marco muestral imaginario (sencillo)
# y entender cómo se construye la ENSANUT
set.seed(2457)
datos <- tibble(
  nombre = c(
    "René", "Alexis", "Francis", "Andy", "Charlie",
    "Andrea", "Sasha", "Guadalupe", "Akira", "Sol",
    "Paris", "Li"
  ),
  edad = sample(
    c("Menor a 30", "Mayor o igual a 30"), 12, T,
    c(0.25, 0.75)
  ),
  sexo = sample(c("Masc", "Fem"), 12, T, c(0.49, 0.51)),
  hogar = sample(1:5, 12, T, c(0.1, 0.2, 0.3, 0.1, 0.3))
)

datos <- datos %>%
  group_by(edad) %>%
  mutate(estatura = if_else(edad == "Menor a 30",
    rnorm(n(), mean = 1.8, sd = 0.1),
    rnorm(n(), mean = 1.6, sd = 0.1)
  )) %>%
  group_by(hogar) %>%
  mutate(localidad = sample(1:4, 1, T, c(0.25, 0.5, 0.1, 0.15))) %>%
  ungroup()

La tabla se ve así:

nombre	edad	sexo	hogar	estatura	localidad
René	Mayor o igual a 30	Fem	2	1.600142	1
Alexis	Menor a 30	Fem	3	1.732241	2
Francis	Menor a 30	Masc	3	1.998119	2
Andy	Mayor o igual a 30	Fem	3	1.580372	2
Charlie	Mayor o igual a 30	Fem	5	1.483513	2
Andrea	Mayor o igual a 30	Fem	2	1.677817	1
Sasha	Menor a 30	Masc	2	1.775873	1
Guadalupe	Mayor o igual a 30	Fem	3	1.451104	2
Akira	Mayor o igual a 30	Fem	1	1.741073	1
Sol	Mayor o igual a 30	Masc	5	1.672928	2
Paris	Menor a 30	Masc	5	1.706307	2
Li	Mayor o igual a 30	Masc	3	1.643326	2

8.2 Breve explicación del muestreo polietápico estratificado

En esta sección discutiremos los tipos de muestreo que existen complicando cada vez un poco más las cosas hasta llegar a un diseño similar al de la ENSANUT.

8.2.1 Muestreo Aleatorio Simple (uniforme) sin reemplazo

Comenzamos con el tipo de muestreo aleatorio más sencillo: el simple. La idea del muestreo aleatorio simple sin reemplazo es de una población de tamaño \(N\) seleccionar \(n\) (donde \(n \leq N\)) individuos (no tienen que ser personas pueden ser animales, árboles, bacterias, piedras) elegidos de manera aleatoria cada uno de ellos.

Podemos pensar que la elección es paso a paso en el sentido de que primero se selecciona \(1\) de todos luego otro (dentro de los \(N-1\) que quedan) y luego otro (dentro de los \(N-2\)) y así… Esto se puede repetir en R en múltiples pasos:

nombres <- datos$nombre

# Extraigo el primer nombre
primer_nombre <- sample(nombres, size = 1)
primer_nombre

[1] "Sol"

# Actualizo los nombres que quedan
nombres <- nombres[which(nombres != primer_nombre)]

# Extraigo el segundo nombre para mi muestra:
segundo_nombre <- sample(nombres, size = 1)
segundo_nombre

[1] "Guadalupe"

# El proceso puede continuar según el tamaño de muestra

Sin embargo, como notarás, es un poco engorroso hacer la extracción de uno en uno. No sólo eso sino que tenemos que muestrear un vector de nombres para luego ir a buscarlos en la tabla. ¡Qué poco óptimo! En R podemos muestrear de la tabla directamente y especificar cuántos queremos que extraiga sin reemplazo (es decir de uno en uno reduciendo el tamaño de la muestra sin que vuelvan a participar). Para ello existe el comando sample_n:

muestra <- datos %>%
  sample_n(size = 7, replace = FALSE)

nombre	edad	sexo	hogar	estatura	localidad
René	Mayor o igual a 30	Fem	2	1.600142	1
Sasha	Menor a 30	Masc	2	1.775873	1
Guadalupe	Mayor o igual a 30	Fem	3	1.451104	2
Li	Mayor o igual a 30	Masc	3	1.643326	2
Charlie	Mayor o igual a 30	Fem	5	1.483513	2
Andy	Mayor o igual a 30	Fem	3	1.580372	2
Alexis	Menor a 30	Fem	3	1.732241	2

Si quisiera calcular estadísticas sobre mi muestra puedo usar las funciones clásicas de dplyr como summarise y group_by para por ejemplo determinar la media de estatura, la proporción de cada sexo o la media de estatura por sexo. Recuerda que la fórmula de la media (promedio) consiste en sumar todos y dividirlos entre \(n\):

\[ \text{Media Clásica} = \frac{\cdot \text{Estatura}_1 + \cdot \text{Estatura}_2 + \dots + \text{Estatura}_n}{n} \]

# Media de estatura
muestra %>%
  summarise(
    Estatura_Promedio = mean(estatura)
  )

# A tibble: 1 × 1
  Estatura_Promedio
              <dbl>
1              1.61

# Proporción de sexos
muestra %>%
  group_by(sexo) %>%
  tally() %>%
  mutate(Proporcion_sexo = n / sum(n))

# A tibble: 2 × 3
  sexo      n Proporcion_sexo
  <chr> <int>           <dbl>
1 Fem       5           0.714
2 Masc      2           0.286

# Media de estatura por sexo
muestra %>%
  group_by(sexo) %>%
  summarise(
    Estatura_Promedio = mean(estatura)
  )

# A tibble: 2 × 2
  sexo  Estatura_Promedio
  <chr>             <dbl>
1 Fem                1.57
2 Masc               1.71

Compara los resultados previos con las métricas poblacionales:

# Media de estatura
datos %>%
  summarise(
    Estatura_Promedio = mean(estatura)
  )

# A tibble: 1 × 1
  Estatura_Promedio
              <dbl>
1              1.67

# Proporción de sexos
datos %>%
  group_by(sexo) %>%
  tally() %>%
  mutate(Proporcion_sexo = n / sum(n))

# A tibble: 2 × 3
  sexo      n Proporcion_sexo
  <chr> <int>           <dbl>
1 Fem       7           0.583
2 Masc      5           0.417

# Media de estatura por sexo
datos %>%
  group_by(sexo) %>%
  summarise(
    Estatura_Promedio = mean(estatura)
  )

# A tibble: 2 × 2
  sexo  Estatura_Promedio
  <chr>             <dbl>
1 Fem                1.61
2 Masc               1.76

8.2.2 Muestreo Aleatorio Simple (ponderado) sin reemplazo

Complicaremos un poco el muestreo. Supongamos que como ser mayor a 30 incluye más grupos de edad que ser menor, nos interesa que las personas mayores a 30 años tengan el doble de probabilidad de aparecer en la muestra que las menores. Para esto, es necesario establecer probabilidades (mal llamadas weights) que le indiquen a R que extraiga con mayor probabilidad 0.66 a los mayores de 30 años que a los menores (proba 0.33):

# Agregamos la columna del ponderador a los datos
# nota que está amañanda para sumar 1
datos <- datos %>%
  mutate(
    probabilidad =
      if_else(edad == "Mayor o igual a 30", 0.667, 0.333)
  )

muestra_ponderada <- datos %>%
  sample_n(size = 7, weight = probabilidad)

nombre	edad	sexo	hogar	estatura	localidad	probabilidad
René	Mayor o igual a 30	Fem	2	1.600142	1	0.667
Charlie	Mayor o igual a 30	Fem	5	1.483513	2	0.667
Akira	Mayor o igual a 30	Fem	1	1.741073	1	0.667
Sol	Mayor o igual a 30	Masc	5	1.672928	2	0.667
Paris	Menor a 30	Masc	5	1.706307	2	0.333
Andy	Mayor o igual a 30	Fem	3	1.580372	2	0.667
Guadalupe	Mayor o igual a 30	Fem	3	1.451104	2	0.667

En este caso ya vale la pena distinguir la media muestral del estimador de la media poblacional. Por un lado la media muestral corresponde al promedio clásico (sumar todos los de la muestra y dividirlos entre el total \(n\)); sin embargo, ésta no es la mejor representación de cómo está la población pues como hubo gente que por el ponderador tenía más probabilidad de salir, éstas personas sesgan la media. Para corregir y dar una mejor foto de cómo se ve la población está el estimador de la media poblacional el cual corresponde al promedio (ponderado) de las estaturas. La fórmula es más o menos así: \[ \text{Media Ponderada} = \frac{w_1\cdot \text{Estatura}_1 + w_2\cdot \text{Estatura}_2 + \dots + w_n\cdot \text{Estatura}_n}{n} \]

En R una forma es usar la función weighted.mean con el ponderador. Usualmente (checa acá) el ponderador para calcular la media es inversamente proporcional a la probabilidad (i.e. \(\text{ponderador} = 1/\text{probabilidad}\)). Esto viene de teoría de encuestas y no entraremos mucho en ella. Sin embargo, brevemente podemos explicar que la lógica es que para aquellos con probabilidad demasiado alta de salir el ponderador \(1/\text{probabilidad}\) da un número más pequeño que para aquellos con probabilidad muy baja de salir (experimenta haciendo \(1/0.99\) contra \(1/0.001\) y ve cómo ponderan cada uno). De alguna forma el ponderador “corrige” que el muestreo tuviera distintas probabilidades.

En R podemos construir el ponderador como sigue:

muestra_ponderada <- muestra_ponderada %>%
  mutate(ponderador = 1 / probabilidad)

Y una forma en la que podemos calcular la media ponderada es con weighted.mean

muestra_ponderada %>%
  summarise(`Media ponderada` = weighted.mean(estatura, ponderador))

# A tibble: 1 × 1
  `Media ponderada`
              <dbl>
1              1.62

Esta es la opción más sencilla pero se queda corta para otros resultados. La opción que nosotros usaremos está dentro del paquete srvyr y es un poco más complicada pero valdrá la pena:

muestra_ponderada %>%
  as_survey_design(1, weight = ponderador) %>%
  summarize(`Media ponderada` = survey_mean(estatura, vartype = "ci"))

# A tibble: 1 × 3
  `Media ponderada` `Media ponderada_low` `Media ponderada_upp`
              <dbl>                 <dbl>                 <dbl>
1              1.62                  1.51                  1.72

Nota que ésta nos da incluso un intervalo de confianza (default 95%). Ésta, por cierto, da un resultado distinto a que si sólo tomamos la media de la muestra:

muestra_ponderada %>%
  summarize(`Media` = mean(estatura))

# A tibble: 1 × 1
  Media
  <dbl>
1  1.61

8.2.3 Muestreo Aleatorio Estratificado sin reemplazo

En nuestra encuesta podría interesarnos garantizar que ciertos grupos estén en la encuesta a fuerza. Como es aleatorio el muestreo siempre puede pasar que tengamos mala suerte y nadie de uno de los grupos esté La solución: pensar nuestra muestra como dos (o más) muestras distintas de distintos grupos (llamados estratos). Por ejemplo podríamos garantizar que a fuerza haya hombres y mujeres en nuestra muestra agrupando por sexo y obteniendo muestras de tamaño 3 ahí:

datos %>%
  group_by(sexo) %>%
  sample_n(size = 3)

# A tibble: 6 × 7
# Groups:   sexo [2]
  nombre    edad               sexo  hogar estatura localidad probabilidad
  <chr>     <chr>              <chr> <int>    <dbl>     <int>        <dbl>
1 Guadalupe Mayor o igual a 30 Fem       3     1.45         2        0.667
2 Andy      Mayor o igual a 30 Fem       3     1.58         2        0.667
3 Akira     Mayor o igual a 30 Fem       1     1.74         1        0.667
4 Sasha     Menor a 30         Masc      2     1.78         1        0.333
5 Sol       Mayor o igual a 30 Masc      5     1.67         2        0.667
6 Paris     Menor a 30         Masc      5     1.71         2        0.333

Si se quisiera un tamaño distinto (digamos \(4\) Masc y \(3\) Fem) puede agregarse una columna de tamaño de muestra y muestrear de acuerdo a esa columna:

# Fuente:
# https://stackoverflow.com/questions/51671856/dplyr-sample-n-by-group-with-unique-size-argument-per-group
muestra_estratificada <- datos %>%
  mutate(tamaño = if_else(sexo == "Fem", 3, 4)) %>%
  group_by(sexo) %>%
  sample_n(size = tamaño[1], weight = probabilidad)

nombre	edad	sexo	hogar	estatura	localidad	probabilidad	ponderador
René	Mayor o igual a 30	Fem	2	1.600142	1	0.667	1.499250
Charlie	Mayor o igual a 30	Fem	5	1.483513	2	0.667	1.499250
Akira	Mayor o igual a 30	Fem	1	1.741073	1	0.667	1.499250
Sol	Mayor o igual a 30	Masc	5	1.672928	2	0.667	1.499250
Paris	Menor a 30	Masc	5	1.706307	2	0.333	3.003003
Andy	Mayor o igual a 30	Fem	3	1.580372	2	0.667	1.499250
Guadalupe	Mayor o igual a 30	Fem	3	1.451104	2	0.667	1.499250

En este caso podemos especificar también cuáles fueron los estratos dentro del diseño:

muestra_ponderada %>%
  as_survey_design(1, strata = sexo, weight = ponderador) %>%
  summarize(`Media ponderada` = survey_mean(estatura, vartype = "se"))

# A tibble: 1 × 2
  `Media ponderada` `Media ponderada_se`
              <dbl>                <dbl>
1              1.62               0.0353

Aquí la opción se representa el error estándar. Nota que el estrato influye en el error estándar (y la varianza) pues de no especificarlo obtenemos un número distinto:

muestra_ponderada %>%
  as_survey_design(1, weight = ponderador) %>%
  summarize(`Media ponderada` = survey_mean(estatura, vartype = "se"))

# A tibble: 1 × 2
  `Media ponderada` `Media ponderada_se`
              <dbl>                <dbl>
1              1.62               0.0421

La idea de una muestra estratificada no sólo es que los estratos representen grupos poblacionales sino que si los estratos se arman “adecuadamente” (como aquí) podemos reducir el error estándar de nuestro estimador y por tanto la longitud de los intervalos de confianza.

Podemos también usar la variable de agrupación para ver las medias por sexo:

muestra_ponderada %>%
  as_survey_design(1, strata = sexo, weight = ponderador) %>%
  group_by(sexo) %>%
  summarize(`Media ponderada` = survey_mean(estatura, vartype = "ci"))

# A tibble: 2 × 4
  sexo  `Media ponderada` `Media ponderada_low` `Media ponderada_upp`
  <chr>             <dbl>                 <dbl>                 <dbl>
1 Fem                1.57                  1.44                  1.70
2 Masc               1.70                  1.66                  1.73

8.2.4 Muestreo Aleatorio Bietápico

Usualmente en las encuestas existe más de un nivel de aleatoriedad. Por ejemplo, se selecciona aleatoriamente una vivienda (unidad primaria) y luego se selecciona aleatoriamente una persona (unidad secundaria) dentro de la vivienda. Esto porque de inicio se desconocen cuántas personas hay por vivienda. Para muestrear de esta forma en nuestro código en el caso donde haya que muestrear \(3\) viviendas y luego máximo dos personas de cada estrato (sexo) por vivienda:

# Primero muestreo n viviendas
viviendas <- datos %>%
  distinct(hogar) %>%
  sample_n(3)

# Luego muestreo personas por vivienda
muestra_hogar <- datos %>%
  filter(hogar %in% !!viviendas$hogar) %>%
  group_by(hogar, sexo) %>%
  sample_n(min(2, n()), weight = probabilidad)

Nota que en el primer hogar sólo se muestreo una persona pues no vivía nadie más ahí. Esto va a requerir un ajuste dado que el individuo está solitario (no hay uno de los estratos en la unidad). Hay varias opciones en R para manejar esto. Para la ENSANUT nosotros pediremos que se ajuste la varianza de dicho estrato (para los intervalos) usando las varianzas de los demás con la siguiente opción:

options(survey.lonely.psu = "adjust")

Una vez establecida la opción podemos indicar a R el diseño de la encuesta estableciendo cuáles fueron las unidades primarias de muestreo:

# Nota la opción de psu no es necesaria en R ponerla pero
# si se requieren opciones más avanzadas de Bootstrap/Jacknife sí
muestra_ponderada %>%
  as_survey_design(1, strata = sexo, weight = ponderador, psu = upm) %>%
  group_by(sexo) %>%
  summarize(`Media ponderada` = survey_mean(estatura, vartype = "ci"))

# A tibble: 2 × 4
  sexo  `Media ponderada` `Media ponderada_low` `Media ponderada_upp`
  <chr>             <dbl>                 <dbl>                 <dbl>
1 Fem                1.57                  1.44                  1.70
2 Masc               1.70                  1.66                  1.73

Para no tener que estar llamando el diseño de la encuesta cada que operamos podemos guardar el objeto como una tabla tipo encuesta:

# Antes de volverla encuesta
muestra_ponderada

# A tibble: 7 × 8
  nombre    edad          sexo  hogar estatura localidad probabilidad ponderador
  <chr>     <chr>         <chr> <int>    <dbl>     <int>        <dbl>      <dbl>
1 René      Mayor o igua… Fem       2     1.60         1        0.667       1.50
2 Charlie   Mayor o igua… Fem       5     1.48         2        0.667       1.50
3 Akira     Mayor o igua… Fem       1     1.74         1        0.667       1.50
4 Sol       Mayor o igua… Masc      5     1.67         2        0.667       1.50
5 Paris     Menor a 30    Masc      5     1.71         2        0.333       3.00
6 Andy      Mayor o igua… Fem       3     1.58         2        0.667       1.50
7 Guadalupe Mayor o igua… Fem       3     1.45         2        0.667       1.50

# Después
muestra_ponderada <- muestra_ponderada %>%
  as_survey_design(1, strata = sexo, weight = ponderador, psu = upm)
muestra_ponderada

Stratified Independent Sampling design (with replacement)
Called via srvyr
Sampling variables:
  - ids: `1` 
  - strata: sexo 
  - weights: ponderador 
Data variables: 
  - nombre (chr), edad (chr), sexo (chr), hogar (int), estatura (dbl),
    localidad (int), probabilidad (dbl), ponderador (dbl)

En general los comandos de dplyr tipo mutate, summarise, rename están disponibles para las tablas de encuesta. Sin embargo no todas las funciones que existen en R para tibbles estàn. La sugerencia es volver tu tabla tipo encuesta una vez hayas limpiado todo y ya sólo falte analizar.

8.3 Análisis de la ENSANUT

La ENSANUT 2021 (ver Romero-Martı́nez et al. (2021) y Martı́nez et al. (2021)) es una encuesta probabilística polietápica estratificada. La forma de elaborarla fue obteniendo una muestra aleatoria de Áreas Geoestadísticas Básicas (AGEB). Éstas fueron las Unidades Primarias de Muestreo (UPM) que se seleccionaron con ponderadores de acuerdo a su población. Una vez seleccionada la AGEB, las Unidades Secundarias de Muestreo (USM) fueron las manzanas seleccionadas aleatoriamente. Finalmente, dentro de cada manzana se seleccionaron varias viviendas (\(n = 6\)). Una vez seleccionadas las viviendas, dentro de cada vivienda se seleccionaron algunas personas de acuerdo a la edad para obtener al menos una por cada grupo etario de interés. Finalmente, para mediciones extras (antropometría / sangre / etc) se obtuvieron submuestras de personas a las que se les aplicó un cuestionario extra.

Romero-Martı́nez, Martı́n, Tonatiuh Barrientos-Gutiérrez, Lucı́a Cuevas-Nasu, et al. 2021. «Metodologı́a de la Encuesta Nacional de Salud y Nutrición 2020 sobre Covid-19». salud pública de méxico 63 (3 May-Jun): 444-51.

Martı́nez, Martı́n Romero, Tonatiuh Barrientos-Gutiérrez, Lucia Cuevas-Nasu, et al. 2021. «Metodologı́a de la Encuesta Nacional de Salud y Nutrición 2021». Salud Pública de México 63 (6, Nov-Dic): 813-18.

8.3.1 Cuestionario de adultos

En nuestro caso para el de adultos se establece así la estructura de encuesta:

# Armamos nuestra tabla con estructura de encuesta
# Nota la opción de psu no es necesaria en R ponerla pero
# si se requieren opciones más avanzadas de Bootstrap/Jacknife sí
options(survey.lonely.psu = "adjust") # singleunit(centered) en Stata
salud_adultos_svy <- salud_adultos %>%
  filter(!is.na(ponde_f)) %>%
  as_survey_design(
    id = FOLIO_INT, strata = est_sel,
    psu = upm, weights = ponde_f, nest = TRUE
  )

Una vez armada la estructura podemos empezar a hacerle preguntas a la base. Por ejemplo a qué porcentaje de personas les dijeron que tenían diabetes o niveles de azucar altos:

salud_adultos_svy %>%
  summarise(
    Diabetes = survey_mean(a0301 == 1, na.rm = T)
  )

# A tibble: 1 × 2
  Diabetes Diabetes_se
     <dbl>       <dbl>
1    0.102     0.00359

En hombres y mujeres así se vieron:

salud_adultos_svy %>%
  mutate(sexo_name = if_else(sexo == 1, "Hombre", "Mujer")) %>%
  group_by(sexo_name) %>%
  summarise(
    Diabetes = survey_mean(a0301 == 1, na.rm = T, vartype = "ci")
  )

# A tibble: 2 × 4
  sexo_name Diabetes Diabetes_low Diabetes_upp
  <chr>        <dbl>        <dbl>        <dbl>
1 Hombre      0.0903       0.0793        0.101
2 Mujer       0.113        0.104         0.122

Mientras que una tabla por edad y sexo es así:

salud_adultos_svy %>%
  mutate(sexo_name = if_else(sexo == 1, "Hombre", "Mujer")) %>%
  mutate(edad_name = case_when(
    edad < 40 ~ "20 a 39",
    edad >= 40 & edad < 60 ~ "49 a 59",
    edad >= 60 ~ "60 y más",
  )) %>%
  group_by(sexo_name, edad_name) %>%
  summarise(
    Diabetes = survey_mean(a0301 == 1, na.rm = T, vartype = "ci")
  )

# A tibble: 6 × 5
# Groups:   sexo_name [2]
  sexo_name edad_name Diabetes Diabetes_low Diabetes_upp
  <chr>     <chr>        <dbl>        <dbl>        <dbl>
1 Hombre    20 a 39     0.0218       0.0103       0.0333
2 Hombre    49 a 59     0.103        0.0851       0.121 
3 Hombre    60 y más    0.229        0.195        0.264 
4 Mujer     20 a 39     0.0196       0.0139       0.0253
5 Mujer     49 a 59     0.149        0.132        0.167 
6 Mujer     60 y más    0.281        0.252        0.310

Si quisiéramos conteos en lugar de utilizar survey_mean podemos usar el survey_total:

salud_adultos_svy %>%
  mutate(sexo_name = if_else(sexo == 1, "Hombre", "Mujer")) %>%
  mutate(edad_name = case_when(
    edad < 40 ~ "20 a 39",
    edad >= 40 & edad < 60 ~ "49 a 59",
    edad >= 60 ~ "60 y más",
  )) %>%
  group_by(sexo_name, edad_name) %>%
  summarise(
    N = survey_total(a0301 == 1, na.rm = T, vartype = "ci")
  )

# A tibble: 6 × 5
# Groups:   sexo_name [2]
  sexo_name edad_name        N    N_low    N_upp
  <chr>     <chr>        <dbl>    <dbl>    <dbl>
1 Hombre    20 a 39    412720.  191303.  634136.
2 Hombre    49 a 59   1400260. 1145975. 1654544.
3 Hombre    60 y más  1850151. 1530672. 2169630.
4 Mujer     20 a 39    399334.  282878.  515790.
5 Mujer     49 a 59   2374376. 2078941. 2669812.
6 Mujer     60 y más  2233392. 1971571. 2495213.

Y si queremos resultados de la muestra hay que usar unweighted. ¡Es posible en summarise combinar múltiples objetos:

salud_adultos_svy %>%
  mutate(sexo_name = if_else(sexo == 1, "Hombre", "Mujer")) %>%
  mutate(edad_name = case_when(
    edad < 40 ~ "20 a 39",
    edad >= 40 & edad < 60 ~ "49 a 59",
    edad >= 60 ~ "60 y más",
  )) %>%
  group_by(sexo_name, edad_name) %>%
  summarise(
    N  = survey_total(a0301 == 1, na.rm = T, vartype = "ci"),
    n  = unweighted(sum(a0301 == 1, na.rm = T))
  )

# A tibble: 6 × 6
# Groups:   sexo_name [2]
  sexo_name edad_name        N    N_low    N_upp     n
  <chr>     <chr>        <dbl>    <dbl>    <dbl> <int>
1 Hombre    20 a 39    412720.  191303.  634136.    41
2 Hombre    49 a 59   1400260. 1145975. 1654544.   212
3 Hombre    60 y más  1850151. 1530672. 2169630.   261
4 Mujer     20 a 39    399334.  282878.  515790.    72
5 Mujer     49 a 59   2374376. 2078941. 2669812.   485
6 Mujer     60 y más  2233392. 1971571. 2495213.   519

Podemos filtrar para obtener la tabla para una entidad en particular. Por ejemplo para la Ciudad de México (entidad == "09")

salud_adultos_svy %>%
  filter(entidad == "09") %>%
  mutate(sexo_name = if_else(sexo == 1, "Hombre", "Mujer")) %>%
  mutate(edad_name = case_when(
    edad < 40 ~ "20 a 39",
    edad >= 40 & edad < 60 ~ "49 a 59",
    edad >= 60 ~ "60 y más",
  )) %>%
  group_by(sexo_name, edad_name) %>%
  summarise(
    N  = survey_total(a0301 == 1, na.rm = T, vartype = "ci"),
    n  = unweighted(sum(a0301 == 1, na.rm = T))
  )

# A tibble: 6 × 6
# Groups:   sexo_name [2]
  sexo_name edad_name       N   N_low   N_upp     n
  <chr>     <chr>       <dbl>   <dbl>   <dbl> <int>
1 Hombre    20 a 39    12860.  -3471.  29191.     3
2 Hombre    49 a 59   132681.  76513. 188849.    26
3 Hombre    60 y más  189190. 119159. 259221.    38
4 Mujer     20 a 39    20803.   3596.  38010.     6
5 Mujer     49 a 59   169581. 109481. 229682.    41
6 Mujer     60 y más  220766. 154535. 286998.    59

8.4 Ejercicios

Responde las siguientes preguntas:

¿Qué proporción de personas tienen diagnóstico previo de enfermedad renal? ¿Varía entre hombres y mujeres?
Replica los resultados de la siguiente tabla del Porcentaje de adultos que reportan medición de colesterol en la sangre, y tenían niveles altos (A0604):

Grupo de edad	Proporción de Hombres	Total (muestra) Hombres	Proporción de Mujeres	Total (muestra) Mujeres
20 a 39	6% [6%, 8%]	1,316,391 (149)	8% [6%, 10%]	1,603,452 (261)
49 a 59	18% [14%, 20%]	2,324,265 (308)	22% [20%, 24%]	3,503,055 (682)
60 y más	16% [14%, 20%]	1,349,788 (194)	24% [22%, 28%]	1,951,682 (431)

¿En qué entidad federativa más gente tiene una vacuna contra el tétanos (A0906)?
De las personas que se accidentaron en vehículos de cuatro o más ruedas (A1102), ¿cuántas llevaban puesto el cinturón de seguridad (A1103)?
Los pacientes que toman pastillas para controlar su insulina, en promedio, ¿cuánto tiempo tienen tomándolas?

Nota Apóyate del paquete lubridate para sumar meses más años. Por ejemplo para sumar \(3\) meses más \(5\) años se haría así:

library(lubridate)
suma_periodo <- months(3) + years(5) # En meses: 5*12 + 3 = 63

# Para convertir cambia unit a:
# meses (months) / días (days) / años (years) / semanas (weeks)
time_length(suma_periodo, unit = "months")

[1] 63

8.4.1 Cuestionario de antropometría

Armamos el diseño de cuestionaro

antropometria_cuest <- antropometria %>%
  left_join(salud_adultos,
    by = c("FOLIO_INT", "upm", "est_sel")
  ) %>%
  rename(ponde_f = ponde_f.x) # Nos quedamos con el ponderador adecuado

# Armamos nuestra tabla con estructura de encuesta
options(survey.lonely.psu = "adjust") # singleunit(centered) en Stata
antropometria_svy <- antropometria_cuest %>%
  filter(!is.na(ponde_f)) %>%
  as_survey_design(
    id = FOLIO_INT, strata = est_sel,
    psu = upm, weights = ponde_f, nest = TRUE
  )

El promedio de peso en menores de 60 (quitando a las embarazadas) puede calcularse como siempre con svy_mean:

antropometria_svy %>%
  filter(an06 != 1 & an06 != 3) %>% # Embarazadas
  summarise(Promedio_peso = survey_mean(an01_1, na.rm = T, "ci"))

# A tibble: 1 × 3
  Promedio_peso Promedio_peso_low Promedio_peso_upp
          <dbl>             <dbl>             <dbl>
1          67.1              66.4              67.9

No sólo podemos calcular la media sino también la mediana y los cuantiles (para por ejemplo calcular dónde está el 25% y el 75%)

antropometria_svy %>%
  filter(an06 != 1 & an06 != 3) %>% # Embarazadas
  summarise(
    Mediana_peso = survey_median(an01_1, na.rm = T, "se"),
    Q_peso = survey_quantile(an01_1,
      quantiles = c(0.25, 0.75),
      na.rm = T, "se"
    )
  )

# A tibble: 1 × 6
  Mediana_peso Mediana_peso_se Q_peso_q25 Q_peso_q75 Q_peso_q25_se Q_peso_q75_se
         <dbl>           <dbl>      <dbl>      <dbl>         <dbl>         <dbl>
1         65.0           0.395       55.2       75.8         0.434         0.472

Podemos dentro del mismo paquete hacer algunas pruebas estadísticas (por ahora las haremos, luego las explicamos) como una prueba t de S tudent para diferencia de medias:

antropometria_svy %>%
  svyttest(
    design = ., # Notación para decirle que el diseño ya está
    formula = an01_1 ~ sexo, # Peso por sexo
    na.rm = TRUE
  )


    Design-based t-test

data:  an01_1 ~ sexo
t = -13.244, df = 6180, p-value < 2.2e-16
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
 -12.561163  -9.321981
sample estimates:
difference in mean 
         -10.94157

Podemos crear nuevas variables como en tidyverse por ejemplo una del índice cintura altura que combina la medida de la cintura (an08_1) con la estatura (an04_1) en personas de 20 a 60 años:

antropometria_svy <- antropometria_svy %>%
  mutate(ica = ifelse(edad > 20 & edad < 60, an08_1 / an04_1, NA)) %>%
  mutate(
    interpreta_ica = case_when(
      ica < 0.35 ~ "Extremadamente delgado",
      ica >= 0.35 & ica < 0.43 & sexo == 1 ~ "Delgado sano",
      ica >= 0.35 & ica < 0.42 & sexo == 2 ~ "Delgado sano",
      ica >= 0.43 & ica < 0.53 & sexo == 1 ~ "Sano",
      ica >= 0.42 & ica < 0.49 & sexo == 2 ~ "Sano",
      ica >= 0.53 & ica < 0.58 & sexo == 1 ~ "Sobrepeso",
      ica >= 0.49 & ica < 0.54 & sexo == 2 ~ "Sobrepeso",
      ica >= 0.58 & ica < 0.63 & sexo == 1 ~ "Sobrepeso elevado",
      ica >= 0.54 & ica < 0.58 & sexo == 2 ~ "Sobrepeso elevado",
      ica >= 0.63 & sexo == 1 ~ "Obesidad mórbida",
      ica >= 0.58 & sexo == 2 ~ "Obesidad mórbida",
      TRUE ~ NA_character_
    )
  )

antropometria_svy %>%
  svychisq(
    design = ., # Notación para decirle que el diseño ya está
    formula = ~ interpreta_ica + sexo, # Peso por sexo
    na.rm = TRUE,
    statistic = "Wald"
  )


    Design-based Wald test of association

data:  NextMethod()
F = 81.34, ndf = 5, ddf = 15835, p-value < 2.2e-16

En este caso, por ejemplo, nuestro valor \(p < \alpha\) nos indica que hay una asociación entre las variables.

8.4.2 Ejercicios

Determine la media, mediana y cuantiles de índice de masa corporal en el cuestionario de antropometría. Recuerda que índice de masa corporal se calcula como: \[ \text{IMC} = \frac{Peso (kg)}{(Altura (m))^2} \]

No olvides convertir la altura en metros pues está en centímetros.
Clasifique a las personas en normal, sobrepeso, peso bajo y obesidad y determine la proporción de individuos (hombres y mujeres) en cada categoría.
Determine si hay una asociación entre IMC y sexo mediante una prueba t.
Determine si hay una asociación entre IMC y entidad mediante una prueba ji cuadrada \(\chi^2\).
Averigüe como utilizar svyglm para generar una regresión de peso contra altura con las covariables de sexo y edad. Como recomendación cheque el siguiente curso: https://tidy-survey-r.github.io/tidy-survey-short-course/

8.5 Modelos de regresión con diseño de encuesta

Hasta aquí hemos calculado promedios, totales y medianas. Ahora vamos a ajustar modelos, que es lo que contesta preguntas del tipo “¿cuánto sube la presión por cada año de edad, si mantenemos el IMC constante?”.

La regla de oro: svyglm, nunca glm

Si usaras glm con estos datos estarías tratando a la ENSANUT como si fuera una muestra aleatoria simple. Los coeficientes te saldrían parecidos, pero los errores estándar estarían mal (casi siempre demasiado chicos), y con ellos los valores p y los intervalos de confianza.

svyglm usa el mismo formula = y ~ x de siempre; lo único que cambia es que en vez de data = recibe design =. Todo lo que aprendiste en el capítulo de regresiones sigue aplicando.

8.5.1 Preparar la base analítica

Vamos a juntar el cuestionario de adultos con el de antropometría y a construir las variables que necesitamos. La parte más importante de este bloque no es la regresión: es la limpieza.

library(svydiags) # diagnósticos de regresión que respetan el diseño
library(WeightedROC) # curva ROC con ponderadores

base_modelos <- antropometria %>%
  left_join(salud_adultos, by = c("FOLIO_INT", "upm", "est_sel")) %>%
  rename(ponde_f = ponde_f.x) %>%
  mutate(
    # Presión arterial: 999 es código de NO RESPUESTA, no una presión de 999
    across(
      c(an27_01s, an27_02s, an27_03s, an27_01d, an27_02d, an27_03d),
      ~ if_else(.x >= 900, NA_real_, as.numeric(.x))
    ),

    # Peso y talla: 222.2 también es código de no respuesta
    peso = if_else(an01_1 > 200, NA_real_, as.numeric(an01_1)),
    talla = if_else(an04_1 > 200 | an04_1 < 100, NA_real_, as.numeric(an04_1)),
    imc = peso / (talla / 100)^2,

    # Se descarta la 1a toma de presión y se promedian la 2a y la 3a,
    # que es la práctica estándar: la primera suele salir alta por el nervio
    sistolica = rowMeans(cbind(an27_02s, an27_03s), na.rm = TRUE),
    diastolica = rowMeans(cbind(an27_02d, an27_03d), na.rm = TRUE),

    # Conteo: ¿en cuántas de las 3 tomas salió la presión elevada?
    tomas_altas = (an27_01s >= 140 | an27_01d >= 90) +
      (an27_02s >= 140 | an27_02d >= 90) +
      (an27_03s >= 140 | an27_03d >= 90),
    diabetes = if_else(a0301 == 1, 1, 0),
    sexo_f = factor(if_else(sexo == 1, "Hombre", "Mujer")),
    edad = as.numeric(edad)
  ) %>%
  filter(edad >= 20, !is.na(ponde_f))

Los códigos de no respuesta arruinan silenciosamente tu modelo

En la ENSANUT el 999 de presión arterial y el 222.2 de peso y talla no son mediciones: son la forma de decir “no se pudo medir”. Si no los quitas, R los va a tratar como números reales y vas a estimar la presión de gente con 999 mmHg y el IMC de personas de 2.22 metros.

Nadie te avisa. El modelo corre, da coeficientes, y están mal. Siempre grafica o saca el rango (range) de tus variables antes de modelar: un máximo sospechosamente redondo o repetido es la señal.

Armamos el diseño como siempre:

options(survey.lonely.psu = "adjust")

diseno_modelos <- base_modelos %>%
  as_survey_design(
    id = FOLIO_INT, strata = est_sel,
    psu = upm, weights = ponde_f, nest = TRUE
  )

8.5.2 Regresión lineal: presión sistólica

La pregunta: ¿cómo cambia la presión sistólica con la edad y el IMC?

modelo_lineal <- svyglm(sistolica ~ edad + imc + sexo_f,
  design = diseno_modelos
)

tidy(modelo_lineal, conf.int = TRUE)

# A tibble: 4 × 7
  term        estimate std.error statistic  p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 (Intercept)   93.1      1.48        62.7 0          90.1      96.0  
2 edad           0.392    0.0257      15.3 9.65e-52    0.342     0.443
3 imc            0.507    0.0494      10.3 1.40e-24    0.411     0.604
4 sexo_fMujer  -10.5      0.526      -19.9 2.77e-85  -11.5      -9.43

Cómo se leen estos coeficientes

edad: por cada año más, la presión sistólica sube en promedio 0.39 mmHg, manteniendo IMC y sexo constantes. Sobre 40 años de diferencia, eso son unos 16 mmHg: muchísimo.
imc: por cada unidad de IMC, la presión sube 0.51 mmHg.
sexo_fMujer: las mujeres tienen en promedio 10.5 mmHg menos que los hombres de la misma edad e IMC.

Nota que Hombre no aparece: es la categoría de referencia (R toma la primera en orden alfabético). Todos los demás coeficientes se leen contra ella.

Visualizar el modelo

Un modelo con tres variables no se ve en una sola gráfica de puntos. Lo que se hace es predecir para valores que a uno le interesen y graficar eso:

# Una malla de valores: edades de 20 a 90, con IMC fijo en tres niveles
malla <- expand_grid(
  edad   = seq(20, 90, by = 1),
  imc    = c(22, 27, 32), # normal, sobrepeso, obesidad
  sexo_f = factor(c("Hombre", "Mujer"))
)

pred <- predict(modelo_lineal, newdata = malla, se.fit = TRUE)

malla <- malla %>%
  mutate(
    ajuste = as.numeric(pred),
    ee = sqrt(attr(pred, "var")),
    bajo = ajuste - 1.96 * ee,
    alto = ajuste + 1.96 * ee,
    imc_et = factor(imc,
      levels = c(22, 27, 32),
      labels = c(
        "IMC 22 (normal)", "IMC 27 (sobrepeso)",
        "IMC 32 (obesidad)"
      )
    )
  )

ggplot(malla, aes(x = edad, y = ajuste, color = imc_et, fill = imc_et)) +
  geom_ribbon(aes(ymin = bajo, ymax = alto), alpha = 0.2, color = NA) +
  geom_line(linewidth = 1) +
  geom_hline(yintercept = 140, linetype = "dashed", color = "gray40") +
  annotate("text",
    x = 26, y = 143, label = "Umbral de hipertensión",
    color = "gray40", size = 3
  ) +
  facet_wrap(~sexo_f) +
  labs(
    x = "Edad (años)", y = "Presión sistólica predicha (mmHg)",
    color = "", fill = "",
    title = "Presión sistólica según edad, IMC y sexo",
    caption = "ENSANUT 2021. Estimación con diseño de encuesta."
  ) +
  theme_bw() +
  theme(legend.position = "bottom")

Presión sistólica predicha por el modelo, con su intervalo de confianza al 95%.

Nota

Esta gráfica es mucho más informativa que la tabla de coeficientes, y es la que va en un informe. Se ve de inmediato que la edad pesa más que el IMC (las líneas suben mucho más de lo que se separan entre sí) y que los hombres arrancan más alto.

8.5.3 Regresión de Poisson: un conteo

La regresión de Poisson se usa cuando el desenlace es un conteo: número de casos, de visitas, de eventos. Aquí usaremos cuántas de las tres tomas de presión salieron elevadas (tomas_altas, que va de 0 a 3).

modelo_poisson <- svyglm(tomas_altas ~ edad + imc + sexo_f,
  design = diseno_modelos,
  family = quasipoisson()
)

tidy(modelo_poisson, conf.int = TRUE)

# A tibble: 4 × 7
  term        estimate std.error statistic  p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 (Intercept)  -3.78     0.213      -17.7  1.61e-68  -4.19     -3.36  
2 edad          0.0407   0.00378     10.8  8.23e-27   0.0333    0.0481
3 imc           0.0513   0.00490     10.5  1.91e-25   0.0417    0.0609
4 sexo_fMujer  -0.558    0.0783      -7.14 1.08e-12  -0.712    -0.405

En Poisson los coeficientes están en escala logarítmica, así que hay que exponenciarlos para interpretarlos. Al resultado se le llama razón de tasas (o IRR):

exp(cbind(IRR = coef(modelo_poisson), confint(modelo_poisson))) %>%
  round(3)

              IRR 2.5 % 97.5 %
(Intercept) 0.023 0.015  0.035
edad        1.042 1.034  1.049
imc         1.053 1.043  1.063
sexo_fMujer 0.572 0.491  0.667

Cómo se lee una razón de tasas

Un IRR de 1.042 para la edad significa que por cada año más, el número esperado de tomas altas se multiplica por ese factor (o sea, sube alrededor de 4.2%).

A diferencia de la regresión lineal, aquí el efecto es multiplicativo: no se suman casos, se multiplican.

Verificar el supuesto de Poisson

La distribución de Poisson supone algo muy fuerte: que la varianza es igual a la media. Casi nunca es cierto en datos reales, y hay que revisarlo.

summary(modelo_poisson)$dispersion

     variance     SE
[1,]   1.9664 0.1015

Aquí hay sobredispersión

El parámetro de dispersión salió alrededor de 1.97, no de 1. Eso significa que los datos varían casi el doble de lo que la Poisson supone.

Mira la distribución del conteo y vas a entender por qué:


   0    1    2    3 
6134  635  364 1073

No es una Poisson: es casi todo ceros con un montón de treses. La gente con presión alta la tiene alta en las tres tomas, no en una. Es un conteo “pegajoso”.

Por eso usamos quasipoisson() y no poisson(): la quasi estima ese parámetro de dispersión y corrige los errores estándar. Si hubiéramos usado poisson(), los intervalos habrían salido artificialmente angostos.

Regla práctica: en datos de salud pública, empieza siempre con quasipoisson(). Si la dispersión sale cerca de 1, no perdiste nada.

8.5.4 Regresión logística: diabetes

Cuando el desenlace es sí/no se usa la logística. Aquí: ¿tiene diagnóstico previo de diabetes?

modelo_logistico <- svyglm(diabetes ~ edad + imc + sexo_f,
  design = diseno_modelos,
  family = quasibinomial()
)

exp(cbind(OR = coef(modelo_logistico), confint(modelo_logistico))) %>%
  round(3)

               OR 2.5 % 97.5 %
(Intercept) 0.000 0.000  0.000
edad        1.114 1.097  1.131
imc         1.039 1.018  1.061
sexo_fMujer 1.432 1.042  1.966

Tip

Un OR (razón de momios) de 1.114 para la edad significa que por cada año más, los momios de tener diabetes se multiplican por ese número. Como es mayor que 1, la edad aumenta el riesgo. Si el intervalo de confianza no incluye al 1, el efecto es distinguible de “nada”.

8.5.5 El detalle que casi nadie enseña: `OR` no es lo mismo que riesgo

Aquí viene algo que vale oro en salud pública. Podemos ajustar el mismo modelo con Poisson en lugar de logística. Cuando el desenlace es 0/1, eso estima la razón de prevalencias (RP) en vez de la razón de momios:

modelo_rp <- svyglm(diabetes ~ edad + imc + sexo_f,
  design = diseno_modelos,
  family = quasipoisson()
)

exp(cbind(RP = coef(modelo_rp), confint(modelo_rp))) %>%
  round(3)

               RP 2.5 % 97.5 %
(Intercept) 0.000 0.000  0.001
edad        1.101 1.087  1.116
imc         1.033 1.016  1.051
sexo_fMujer 1.354 1.030  1.780

Compara los dos resultados para sexo

	Estimación
`OR` (logística)	1.432
`RP` (Poisson)	1.354

El OR es más grande que la RP. No es un error: es una propiedad matemática. El OR exagera el efecto cuando el desenlace es común, y aquí la diabetes ronda el 13% de la población, que ya es bastante común.

El problema es que casi todo el mundo lee un OR de 1.43 como “tienen 43% más riesgo”, y eso es falso. El número correcto para esa frase es la RP: 35%.

Regla: si el desenlace es raro (menos del 10%), OR y RP se parecen y da casi igual. Si es común, reporta razón de prevalencias, que es la que la gente va a interpretar bien de todos modos.

A esto se le llama regresión de Poisson modificada y es de uso estándar en epidemiología. Con diseño de encuesta sale gratis, porque svyglm ya calcula los errores estándar de forma robusta.

8.6 Cómo verificar tu regresión (más allá de los residuales)

En el capítulo de regresiones revisamos el modelo con gráficas de residuales. Eso sigue sirviendo, pero tiene dos límites: es subjetivo (“¿se ve con patrón o no?”) y las funciones normales de diagnóstico ignoran el diseño de muestreo.

Aquí van cinco herramientas que sí respetan el diseño.

8.6.1 1. Comparar modelos con el `AIC`

Igual que en regresión normal, pero AIC sobre un svyglm usa la versión corregida por diseño:

modelo_solo_edad <- svyglm(sistolica ~ edad, design = diseno_modelos)

rbind(
  "Sólo edad"          = AIC(modelo_solo_edad),
  "Edad + IMC + sexo"  = AIC(modelo_lineal)
) %>% round(1)

                  eff.p     AIC deltabar
Sólo edad           2.8 69167.0      1.4
Edad + IMC + sexo   7.9 52772.1      2.0

Entre más chico, mejor. La columna eff.p es el número efectivo de parámetros: con diseño complejo no es simplemente cuántas variables metiste.

8.6.2 2. Prueba de Wald para un término (`regTermTest`)

¿La variable imc aporta algo, más allá de lo que ya explica la edad?

regTermTest(modelo_lineal, ~imc)

Wald test for imc
 in svyglm(formula = sistolica ~ edad + imc + sexo_f, design = diseno_modelos)
F =  105.6474  on  1  and  6090  df: p= < 2.22e-16

Nota

regTermTest es especialmente útil para variables categóricas con varios niveles, donde no basta mirar el valor p de cada categoría por separado: te da una sola prueba para “¿esta variable, en conjunto, aporta?”.

8.6.3 3. Comparar modelos anidados con `anova`

anova(modelo_solo_edad, modelo_lineal)

Working (Rao-Scott+F) LRT for imc sexo_f
 in svyglm(formula = sistolica ~ edad + imc + sexo_f, design = diseno_modelos)
Working 2logLR =  441.837 p= < 2.22e-16 
(scale factors:  1.1 0.87 );  denominator df= 6090

Ésta es la versión Rao-Scott de la prueba de razón de verosimilitudes, que es la que corresponde cuando hay diseño de encuesta. La versión normal de anova daría un resultado equivocado.

8.6.4 4. Multicolinealidad con diseño: `svyvif`

El VIF que vimos en el capítulo de regresiones no toma en cuenta el diseño. El paquete svydiags tiene la versión que sí:

X_modelo <- model.matrix(modelo_lineal)[, -1, drop = FALSE]
filas_usadas <- as.integer(rownames(modelo_lineal$model))
pesos <- weights(diseno_modelos, "sampling")[filas_usadas]

vif_diseno <- svyvif(mobj = modelo_lineal, X = X_modelo, w = pesos)

# La función devuelve muchas cosas; nos quedamos con el VIF de diseño
setNames(
  round(unlist(vif_diseno$`Intercept adjusted`$svy.vif.m), 3),
  colnames(X_modelo)
)

       edad         imc sexo_fMujer 
      1.073       1.082       1.004

Con valores tan cerca de 1 no hay problema de colinealidad: edad, IMC y sexo aportan información distinta.

8.6.5 5. Observaciones influyentes: distancia de Cook con diseño

cook <- svyCooksD(modelo_lineal)

c(
  observaciones = length(cook),
  maximo = round(max(cook, na.rm = TRUE), 3),
  por_encima_del_umbral = sum(cook > 4 / length(cook), na.rm = TRUE)
)

        observaciones                maximo por_encima_del_umbral 
             6184.000                 2.491              6184.000

Advertencia

Que haya observaciones influyentes no significa que haya que borrarlas. En una encuesta nacional, una persona con IMC de 60 es una persona real que representa a mucha gente. Lo que hay que hacer es:

Revisar que no sea un error de captura (¿un peso de 222 kg? ¿una talla de 50 cm en un adulto?).
Si el dato es válido, correr el modelo con y sin ella y reportar si cambian las conclusiones.

Borrar datos porque estorban al modelo es de las peores prácticas que existen.

8.6.6 6. Para la logística: la curva ROC ponderada

En un modelo logístico lo que importa no son los residuales sino qué tan bien separa a quienes tienen el desenlace de quienes no. Eso se mide con el área bajo la curva ROC (AUC), y hay que calcularla con los ponderadores:

prob <- as.numeric(predict(modelo_logistico, type = "response"))
obs <- as.numeric(modelo_logistico$y)
filas <- as.integer(rownames(modelo_logistico$model))
w_roc <- weights(diseno_modelos, "sampling")[filas]

roc <- WeightedROC(prob, obs, weight = w_roc)
auc <- WeightedAUC(roc)

ggplot(roc, aes(x = FPR, y = TPR)) +
  geom_line(color = "#003f5c", linewidth = 1) +
  geom_abline(linetype = "dashed", color = "gray50") +
  annotate("text",
    x = 0.6, y = 0.2,
    label = paste0("AUC = ", round(auc, 3)), size = 5
  ) +
  labs(
    x = "1 - Especificidad", y = "Sensibilidad",
    title = "¿Qué tan bien distingue el modelo quién tiene diabetes?"
  ) +
  theme_bw()

Curva ROC ponderada del modelo logístico de diabetes.

Cómo se lee el AUC

Es la probabilidad de que, tomando al azar una persona con diabetes y una sin ella, el modelo le dé mayor probabilidad a la que sí la tiene.

`AUC`	Interpretación
0.5	El modelo no sirve (es una moneda)
0.7 – 0.8	Aceptable
0.8 – 0.9	Bueno
> 0.9	Excelente (o hiciste trampa)

El nuestro da 0.804, que está bien para un modelo con sólo tres variables. La línea punteada es el modelo que adivina al azar: entre más lejos esté tu curva de esa línea, mejor.

8.7 Bootstrap: cuando no existe la función `survey` que necesitas

srvyr tiene survey_mean, survey_total, survey_median y survey_quantile. ¿Pero qué haces cuando necesitas un indicador que no tiene función?

La respuesta es el bootstrap con réplicas. La idea es sencilla:

R genera muchos juegos de ponderadores alternativos (las réplicas), cada uno simulando “qué tal si la muestra hubiera salido un poco distinta”.
Tú calculas tu indicador una vez con cada juego de ponderadores, con un for.
La desviación estándar de todos esos valores es tu error estándar, y los percentiles 2.5 y 97.5 son tu intervalo de confianza.

Lo bonito es que funciona para cualquier indicador, por raro que sea.

8.7.1 Crear las réplicas

Usaremos srs.bootstrap.sample del paquete surveybootstrap. La función es directa: le das tu base y cuántas réplicas quieres.

library(surveybootstrap)

set.seed(2026)

base_imc <- base_modelos %>% filter(!is.na(imc))

replicas <- srs.bootstrap.sample(base_imc, num.reps = 300)

length(replicas)

[1] 300

¿Qué hay dentro de cada réplica? No son los datos, sino instrucciones para rearmarlos:

head(replicas[[1]])

  index weight.scale
1  4829            1
2  3705            1
3   993            1
4  2342            1
5  3629            1
6  1647            1

index: qué renglones de tu base entran en esta muestra bootstrap. Se eligen con reemplazo, así que un mismo renglón puede salir varias veces (o ninguna). Ésa es toda la magia del bootstrap.
weight.scale: un factor por el que hay que multiplicar los ponderadores originales. En el bootstrap srs vale 1.

c(
  renglones = nrow(replicas[[1]]),
  distintos = length(unique(replicas[[1]]$index))
)

renglones distintos 
     6218      3938

Lo que srs supone (y por qué hay que decirlo)

El srs de srs.bootstrap.sample significa simple random sample: la función remuestrea renglones como si la ENSANUT fuera una muestra aleatoria simple, sin tomar en cuenta los estratos ni las UPM.

Es la versión más fácil de entender, y por eso la usamos para aprender el mecanismo. Pero tenlo presente: cuando hay conglomerados, la gente de una misma manzana se parece entre sí y aporta menos información de la que aparenta, así que el error estándar “verdadero” suele ser mayor.

El mismo paquete trae rescaled.bootstrap.sample, que sí respeta los conglomerados y se usa igual (también devuelve index y weight.scale):

replicas_diseno <- rescaled.bootstrap.sample(
  base_imc,
  survey.design = ~upm, # la unidad primaria de muestreo
  num.reps = 300
)

Al final del capítulo comparamos los tres resultados.

8.7.2 Definir los indicadores a mano

Necesitamos versiones ponderadas de lo que queremos calcular:

# Mediana ponderada: el valor donde el peso acumulado llega al 50%
mediana_ponderada <- function(x, w) {
  orden <- order(x)
  x <- x[orden]
  w <- w[orden]
  x[which(cumsum(w) / sum(w) >= 0.5)[1]]
}

# Coeficiente de Gini: qué tan desigual está repartida una variable.
# NO existe survey_gini(), por eso lo hacemos a mano.
gini_ponderado <- function(x, w) {
  orden <- order(x)
  x <- x[orden]
  w <- w[orden]
  p <- cumsum(w) / sum(w)
  cw <- cumsum(w * x) / sum(w * x)
  sum(cw[-1] * p[-length(p)]) - sum(cw[-length(cw)] * p[-1])
}

8.7.3 El `for` sobre las réplicas

n_rep <- length(replicas)

medianas <- numeric(n_rep)
ginis <- numeric(n_rep)

for (i in 1:n_rep) {
  # Qué renglones entran en esta réplica
  filas <- replicas[[i]]$index

  # Sus ponderadores, ajustados por el factor de escala
  w_i <- base_imc$ponde_f[filas] * replicas[[i]]$weight.scale

  # Y el indicador calculado con esa muestra
  medianas[i] <- mediana_ponderada(base_imc$imc[filas], w_i)
  ginis[i] <- gini_ponderado(base_imc$imc[filas], w_i)
}

Nota

Fíjate en que el for no tiene nada de especial: es el mismo ciclo que viste en el capítulo de ciclos. Lo único que cambia en cada vuelta son los renglones que entran y sus ponderadores.

Y ya con eso tenemos todo:

resultados_boot <- tibble(
  Indicador = c("Mediana de IMC", "Gini del IMC"),
  Estimacion = c(
    mediana_ponderada(base_imc$imc, base_imc$ponde_f),
    gini_ponderado(base_imc$imc, base_imc$ponde_f)
  ),
  EE = c(sd(medianas), sd(ginis)),
  IC_bajo = c(quantile(medianas, 0.025), quantile(ginis, 0.025)),
  IC_alto = c(quantile(medianas, 0.975), quantile(ginis, 0.975))
)

resultados_boot

# A tibble: 2 × 5
  Indicador      Estimacion      EE IC_bajo IC_alto
  <chr>               <dbl>   <dbl>   <dbl>   <dbl>
1 Mediana de IMC     28.4   0.123    28.1    28.6  
2 Gini del IMC        0.114 0.00169   0.111   0.117

8.7.4 ¿Cómo sabemos que esto está bien hecho?

Ésta es la parte importante y la que casi siempre se salta. Antes de confiar en el bootstrap para el Gini —que no tiene función directa— hay que comprobarlo en algo que sí tiene función directa. Por eso incluimos la mediana:

diseno_imc <- base_imc %>%
  as_survey_design(
    id = FOLIO_INT, strata = est_sel,
    psu = upm, weights = ponde_f, nest = TRUE
  )

diseno_imc %>%
  summarise(mediana = survey_median(imc, na.rm = TRUE, vartype = "se"))

# A tibble: 1 × 2
  mediana mediana_se
    <dbl>      <dbl>
1    28.4      0.117

La comprobación

Compara los dos resultados de la mediana: la estimación puntual y el error estándar del bootstrap salen prácticamente idénticos a los que da survey_median por su cuenta.

Eso es lo que nos autoriza a creerle al Gini. No lo verificamos porque el número se vea razonable, sino porque el método reprodujo un resultado que ya conocíamos.

Es la misma lógica del paso 6 del flujo de trabajo: antes de usar un procedimiento para lo que no sabes, comprueba que recupera lo que sí sabes.

Otros indicadores que puedes calcular igual

El mismo for sirve, cambiando nada más la función, para: la diferencia de medianas entre dos grupos, la razón P90/P10, el índice de concentración, la asimetría, o cualquier indicador que se te ocurra. Ninguno tiene función directa en srvyr y todos se resuelven con estas quince líneas.

Nota sobre el número de réplicas: 200 alcanza para el error estándar, pero si vas a reportar intervalos de confianza en una publicación, sube a 500 o 1,000. Lo único que cuesta es tiempo de cómputo.

8.8 Ejercicios de modelos y bootstrap

8.8.1 Modelos

Ajusta la regresión lineal de presión diastólica (diastolica) contra edad, IMC y sexo. ¿El efecto de la edad es igual de fuerte que en la sistólica? Grafica las predicciones como lo hicimos arriba.
Al modelo logístico de diabetes agrégale la entidad (desc_ent) como variable. Usa regTermTest para decidir si aporta en conjunto, en lugar de mirar los 31 valores p por separado.
Ajusta el modelo lineal de presión sistólica sin el diseño de encuesta (con glm y data = base_modelos) y compara los errores estándar contra los de svyglm. ¿Cuáles son más chicos? ¿Qué implicaría eso para tus valores p?
Calcula el AUC del modelo logístico quitando el IMC. ¿Cuánto se pierde? ¿Vale la pena la variable?

8.8.2 A mano

Con los coeficientes del modelo lineal, calcula con calculadora la presión sistólica predicha para una mujer de 55 años con IMC de 30. Después compruébalo con predict.
Si el OR de una exposición es 2.0 y el desenlace ocurre en el 40% de la población, ¿la razón de prevalencias será mayor, menor o igual a 2.0? Explica por qué sin hacer cuentas.

8.8.3 Bootstrap

Modifica el for para calcular la diferencia entre la mediana de IMC de mujeres y la de hombres, con su intervalo de confianza. Pista: tu función tiene que recibir también el vector de sexo.
Repite el bootstrap del Gini con replicates = 50 y con replicates = 1000. ¿Cuánto cambia la estimación puntual? ¿Y el intervalo? ¿Qué te dice eso sobre cuántas réplicas necesitas?

8.8.4 Con apoyo de una IA

Pídele a una IA que te explique la diferencia entre OR y RP. Verifica si menciona que depende de qué tan común sea el desenlace. Si no lo menciona, es una omisión grave: pregúntale directamente por ese punto.
Pásale el código del for del bootstrap y pídele que lo reescriba sin ciclo (vectorizado con apply). Comprueba que da exactamente el mismo resultado antes de quedarte con su versión.

8.9 Sistema

sessionInfo()

R version 4.5.3 (2026-03-11)
Platform: x86_64-apple-darwin20
Running under: macOS Sequoia 15.7.3

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] es_ES.UTF-8/es_ES.UTF-8/es_ES.UTF-8/C/es_ES.UTF-8/es_ES.UTF-8

time zone: America/Mexico_City
tzcode source: internal

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] surveybootstrap_0.0.3 WeightedROC_2020.1.31 svydiags_0.7         
 [4] MASS_7.3-65           broom_1.0.12          survey_4.5           
 [7] survival_3.8-6        Matrix_1.7-5          kableExtra_1.4.0     
[10] srvyr_1.3.1           lubridate_1.9.5       forcats_1.0.1        
[13] stringr_1.6.0         dplyr_1.2.1           purrr_1.2.2          
[16] readr_2.2.0           tidyr_1.3.2           tibble_3.3.1         
[19] ggplot2_4.0.3         tidyverse_2.0.0       wesanderson_0.3.7    
[22] haven_2.5.5          

loaded via a namespace (and not attached):
 [1] gtable_0.3.6       xfun_0.60          htmlwidgets_1.6.4  lattice_0.22-9    
 [5] tzdb_0.5.0         vctrs_0.7.3        tools_4.5.3        generics_0.1.4    
 [9] parallel_4.5.3     pkgconfig_2.0.3    RColorBrewer_1.1-3 S7_0.2.2          
[13] lifecycle_1.0.5    compiler_4.5.3     farver_2.1.2       textshaping_1.0.5 
[17] codetools_0.2-20   mitools_2.4        htmltools_0.5.9    yaml_2.3.12       
[21] pillar_1.11.1      crayon_1.5.3       tidyselect_1.2.1   digest_0.6.39     
[25] stringi_1.8.7      labeling_0.4.3     splines_4.5.3      fastmap_1.2.0     
[29] cli_3.6.6          magrittr_2.0.5     utf8_1.2.6         withr_3.0.3       
[33] scales_1.4.0       backports_1.5.1    bit64_4.8.2        timechange_0.4.0  
[37] rmarkdown_2.31     bit_4.6.0          otel_0.2.0         hms_1.1.4         
[41] evaluate_1.0.5     knitr_1.51         viridisLite_0.4.3  rlang_1.3.0       
[45] functional_0.6     Rcpp_1.1.2         glue_1.8.1         DBI_1.3.0         
[49] xml2_1.6.0         svglite_2.2.2     
 [ reached 'max' / getOption("max.print") -- omitted 6 entries ]

8.1 Paquetes y datos a utilizar

8.2 Breve explicación del muestreo polietápico estratificado

8.2.1 Muestreo Aleatorio Simple (uniforme) sin reemplazo

8.2.2 Muestreo Aleatorio Simple (ponderado) sin reemplazo

8.2.3 Muestreo Aleatorio Estratificado sin reemplazo

8.2.4 Muestreo Aleatorio Bietápico

8.3 Análisis de la ENSANUT

8.3.1 Cuestionario de adultos

8.4 Ejercicios

8.4.1 Cuestionario de antropometría

8.4.2 Ejercicios

8.5 Modelos de regresión con diseño de encuesta

8.5.1 Preparar la base analítica

8.5.2 Regresión lineal: presión sistólica

Visualizar el modelo

8.5.3 Regresión de Poisson: un conteo

Verificar el supuesto de Poisson

8.5.4 Regresión logística: diabetes

8.5.5 El detalle que casi nadie enseña: OR no es lo mismo que riesgo

8.6 Cómo verificar tu regresión (más allá de los residuales)

8.6.1 1. Comparar modelos con el AIC

8.6.2 2. Prueba de Wald para un término (regTermTest)

8.6.3 3. Comparar modelos anidados con anova

8.6.4 4. Multicolinealidad con diseño: svyvif

8.6.5 5. Observaciones influyentes: distancia de Cook con diseño

8.6.6 6. Para la logística: la curva ROC ponderada

8.7 Bootstrap: cuando no existe la función survey que necesitas

8.7.1 Crear las réplicas

8.7.2 Definir los indicadores a mano

8.7.3 El for sobre las réplicas

8.7.4 ¿Cómo sabemos que esto está bien hecho?

8.8 Ejercicios de modelos y bootstrap

8.8.1 Modelos

8.8.2 A mano

8.8.3 Bootstrap

8.8.4 Con apoyo de una IA

8.9 Sistema

8.5.5 El detalle que casi nadie enseña: `OR` no es lo mismo que riesgo

8.6.1 1. Comparar modelos con el `AIC`

8.6.2 2. Prueba de Wald para un término (`regTermTest`)

8.6.3 3. Comparar modelos anidados con `anova`

8.6.4 4. Multicolinealidad con diseño: `svyvif`

8.7 Bootstrap: cuando no existe la función `survey` que necesitas

8.7.3 El `for` sobre las réplicas