library(haven)
library(wesanderson)
library(tidyverse)
library(srvyr)
library(kableExtra)
library(survey)Análisis de Encuestas (ENSANUT)
Paquetes y datos a utilizar
A lo largo de esta sección usaremos los siguientes paquetes:
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 |
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.
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 muestraSin 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.71Compara 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.76Muestreo 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í: $$ =
$$
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.62Esta 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.72Nota 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.61Muestreo 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.333Si 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.0353Aquí 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.0421La 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.73Muestreo 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.73Para 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 probabil…¹ ponde…²
<chr> <chr> <chr> <int> <dbl> <int> <dbl> <dbl>
1 René Mayor o igual a 30 Fem 2 1.60 1 0.667 1.50
2 Charlie Mayor o igual a 30 Fem 5 1.48 2 0.667 1.50
3 Akira Mayor o igual a 30 Fem 1 1.74 1 0.667 1.50
4 Sol Mayor o igual a 30 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 igual a 30 Fem 3 1.58 2 0.667 1.50
7 Guadalupe Mayor o igual a 30 Fem 3 1.45 2 0.667 1.50
# … with abbreviated variable names ¹probabilidad, ²ponderador#Después
muestra_ponderada <- muestra_ponderada %>%
as_survey_design(1, strata = sexo, weight = ponderador, psu = upm)
muestra_ponderadaStratified 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.
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.
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.00359En 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.122Mientras 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. 519Podemos 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. 59Ejercicios
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)
Attaching package: 'lubridate'The following objects are masked from 'package:base':
date, intersect, setdiff, unionsuma_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] 63Cuestionario 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.9No 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.472Podemos 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-16En este caso, por ejemplo, nuestro valor \(p < \alpha\) nos indica que hay una asociación entre las variables.
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
svyglmpara 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/
Sistema
sessionInfo()R version 4.2.0 (2022-04-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur/Monterey 10.16
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] grid stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] lubridate_1.8.0 survey_4.1-1 survival_3.4-0 Matrix_1.4-1
[5] kableExtra_1.3.4 srvyr_1.1.1 forcats_0.5.1 stringr_1.4.0
[9] dplyr_1.0.9 purrr_0.3.4 readr_2.1.2 tidyr_1.2.0
[13] tibble_3.1.8 ggplot2_3.3.6 tidyverse_1.3.2 wesanderson_0.3.6
[17] haven_2.5.0
loaded via a namespace (and not attached):
[1] httr_1.4.4 bit64_4.0.5 vroom_1.5.7
[4] jsonlite_1.8.0 viridisLite_0.4.0 splines_4.2.0
[7] modelr_0.1.8 assertthat_0.2.1 highr_0.9
[10] googlesheets4_1.0.1 cellranger_1.1.0 yaml_2.3.5
[13] pillar_1.8.0 backports_1.4.1 lattice_0.20-45
[16] glue_1.6.2 digest_0.6.29 rvest_1.0.2
[19] colorspace_2.0-3 htmltools_0.5.3 pkgconfig_2.0.3
[22] broom_1.0.0 scales_1.2.0 webshot_0.5.3
[25] svglite_2.1.0 tzdb_0.3.0 googledrive_2.0.0
[28] generics_0.1.3 ellipsis_0.3.2 withr_2.5.0
[31] cli_3.3.0 magrittr_2.0.3 crayon_1.5.1
[34] readxl_1.4.0 evaluate_0.16 fs_1.5.2
[37] fansi_1.0.3 xml2_1.3.3 tools_4.2.0
[40] hms_1.1.1 mitools_2.4 gargle_1.2.0
[43] lifecycle_1.0.1 munsell_0.5.0 reprex_2.0.2
[46] compiler_4.2.0 systemfonts_1.0.4 rlang_1.0.4
[49] rstudioapi_0.13 htmlwidgets_1.5.4 rmarkdown_2.15
[52] gtable_0.3.0 DBI_1.1.3 R6_2.5.1
[55] knitr_1.39 fastmap_1.1.0 bit_4.0.4
[58] utf8_1.2.2 stringi_1.7.8 parallel_4.2.0
[61] vctrs_0.4.1 dbplyr_2.2.1 tidyselect_1.1.2
[64] xfun_0.32 References
Martı́nez, Martı́n Romero, Tonatiuh Barrientos-Gutiérrez, Lucia Cuevas-Nasu, Sergio Bautista-Arredondo, M Arantxa Colchero, Elsa Berenice Gaona-Pineda, Jesús Martı́nez-Barnetche, 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.
Romero-Martı́nez, Martı́n, Tonatiuh Barrientos-Gutiérrez, Lucı́a Cuevas-Nasu, Sergio Bautista-Arredondo, Arantxa Colchero, Elsa Berenice Gaona-Pineda, Eduardo Lazcano-Ponce, 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.