Modelización avanzada con R: Análisis de Encuestas.
Doctorado en Economía. Sesión 2ª de 13-XI-2020
José Antonio Ortega
13 de noviembre de 2020
Análisis de encuestas complejas con survey
Introducción a las encuestas complejas
En la sesión de hoy vamos a trabajar con datos de una encuesta compleja, la Encuesta Demográfica y de Salud Familiar de Perú, 2019, que forma parte del programa Demographic and Health Surveys. Podemos verlo como un ejemplo de encuesta compleja, y la analizaremos utilizando el paquete survey.
La filosofía del análisis de encuestas con survey es diferente a la que vimos la semana pasada con las funciones básicas. En vez de tener que especificar en cada estimación las variables de peso o, en su caso, las correspondientes a cada uno de los conglomerados o PSUs, se basa en, ligado a un objeto de datos, declarar inicialmente las características técnicas de la encuesta. Una vez el paquete reconoce a la encuesta, se pueden utilizar cualquiera de las técnicas de análisis y modelización que nos proporciona el paquete sin necesidad de reespecificar la parte técnica.
Antes de continuar vamos a cargar los paquetes que necesitaremos. En caso de que alguno no esté instalado, habrá que instalarlo con install.packages.
library(tidyverse)
library(haven) # Lectura de datos
library(sjmisc) # Utilidades para encuestas con etiquetas
library(sjPlot) # Gráficos y tablas "bonitos" con etiquetas
library(broom) # Convierte en tibbles los resultados de modelos
library(survey) # Análisis de encuestas complejas
library(srvyr) # Funciones survey a la tidyverse
library(knitr) # Tablas con kable
library(sandwich); library(lmtest) # Varianzas-covarianzas
library(car) # Inferencia de modelos¿Por qué es importante hacer las cosas de esta manera? Las encuestas complejas distan mucho de ser un muestreo aleatorio simple, la hipótesis subyacente en la mayor parte de los modelos microeconométricos. Por un lado, la probabilidad de seleccionar cada unidad depende de la estructura de la encuesta. Esto puede controlarse especificando el vector de pesos correcto. Por otro, la estructura de la encuesta influye sobre la variabilidad estimada: nuestra muestra no es una muestra aleatoria, sino que se muestrea típicamente dentro de diferentes Unidades Primarias de Muestreo (PSU), en las cuáles es posible que los individuos sean más cercanos entre sí que en la población en general. Si estimamos la variabilidad como si fuera una muestra aleatoria, en general, estaremos subvalorando la dispersión: los errores estándar correctos son más grandes que los que se darían en una muestra aleatoria del mismo tamaño. A la diferencia entre ambos se le denomina efecto de diseño.
Estructura de encuestas complejas:
Las encuestas complejas suelen usar muestreo aleatorio pero dentro de unidades relativamente pequeñas. Para llegar a estas unidades pueden existir:
- Estratos (strata): Los estratos son subpoblaciones naturales en una encuesta que, a priori, son homogéneos en su interior pero heterogéneos entre sí. El diseño de la encuesta se hace de modo que se garantiza cubrir adecuadamente todos los estratos. Ejemplos típicos: CCAA, urbano/rural, sexo, edad (a veces estos dos últimos se tratan como estratos implícitos o pseudo estratos ajustándose los pesos como describimos más adelante)
- Conglomerados (cluster): Los conglomerados son unidades definidas dentro de cada estrato (si es muestreo polietápico), cuyo tamaño es conocido, pero que no tenemos interés en cubrir por separado. La práctica es seleccionar en una primera etapa una muestra de PSUs con probabilidad proporcional a su tamaño, y extraer una muestra aleatoria dentro de cada PSU.
Los motivos para hacer un muestreo de este tipo son claros: los estratos sirven para garantizar que tengamos resultados representativos de los subgrupos que nos interesen. Los conglomerados, sobre todo, para reducir los costes de elaboración de la encuesta.
Conectado con el diseño muestral están los pesos. Tenemos los siguientes conceptos relacionados:
Probabilidad de inclusión: Es la probabilidad que, a priori, tiene de estar incluida en la muestra cada unidad. Se obtiene, en general, multiplicando las probabilidades asociadas al estrato, a cada unidad primaria de muestreo (PSU) dentro del estrato, y a la unidad de segunda etapa dentro cada PSU.
Peso de diseño: Se calcula como la inversa de la probabilidad de inclusión. Se corresponde con a cuántos individuos representa esa unidad antes de hacer la encuesta.
Pesos post-estratificación: Tras realizar la encuesta la tasa de respuesta puede variar dentro de cada estrato y conglomerado, o incluso de acuerdo a estratos implicitos como pueda ser sexos, grupos de edad. Los pesos post-estratificación se definen de modo que la suma de los pesos en cada estrato se corresponda con la población en el estrato.
Pesos replicados: Algunas encuestas prefieren no revelar cuáles son las unidades de muestreo, pero esta información es en principio necesaria para calcular la variabilidad correctamente. Una alternativa inspirada en el método bootstrap pero adaptado al diseño de la encuesta es el uso de los pesos replicados (replicate weights). La idea es proporcionar una serie de pesos que permitan reestimar los modelos o las tablas y calibrar la incertidumbre en la estimación. Ejemplos: American Community Survey, PISA, …
Se han incluido referencias que pueden servir para profundizar en los temas de diseño y análisis de encuestas, que contienen los detalles técnicos y aplicaciones en R. Jiménez (2014) y Gutiérrez (2017) presentan sendos resumenes de todas las fases de elaboración y análisis de encuestas con R. Lumley (2010) y Mayor (2015) presentan ejemplos de análisis para encuestas de diferentes tipos en R. Cabrera (2016) desarrolla la teoría y ejemplos de las técnicas más recientes como el análisis replicado o la calibración. UCLA (2020) proporciona un ejemplo de análisis completo con el paquete survey similar en espíritu al de esta sesión. Lumley y Scott (2017) presenta una panorámica sobre regresión con datos de encuesta con un ejemplo en R. Por último, Damico(2020) proporciona tutoriales de análisis con el paquete survey para gran número de programas de encuestas, aparte de incluir funciones auxiliares para bajar y analizar los datos en el paquete lodown. Para comparar cómo se analizan encuestas con survey en R con otros paquetes estadísticos, ver CDC(2016). Finalmente, UN Statistical Division (2005) presenta la teoría sobre el diseño de encuestas de hogares con ejemplos de programas de países en desarrollo, incluidas las DHS.
Encuesta Demográfica y de Salud (ENDES-DHS) del Perú (2019): Lectura de datos y preparación para el análisis con survey
Bajada de datos al ordenador
Los microdatos de la ENDES están disponible en la página web del INEI peruano, además de estarlas bajo registro en las del programa DHS. Debemos decir que los ficheros no son iguales: por un lado, los del INEI están en español, por otro, en vez de proporcionar todos los items de cada cuestionario en un fichero, se han dividido de acuerdo a los distintos módulos del cuestionario. Para nuestros efectos, sin embargo, la gran ventaja de utilizar los datos del INEI es que nos permiten hacer replicable el análisis sin necesidad de registrarnos en las páginas del DHS.
Bajamos los datos con una estrategia parecida a la de la sesión anterior: definiendo un objeto con las direcciones url que después bajamos a la carpeta data/ que tenemos creada previamente. A diferencia del otro caso, no respetamos los nombres de los ficheros originales sino que utilizamos nombres más descriptivos. Esto nos obliga a utilizar la función pwalk (o, alternativamente, walk2) de purrr en vez de walk al tener dos argumentos.
urlbase="http://iinei.inei.gob.pe/iinei/srienaho/descarga/SPSS/"
# Generamos una tiabla con las extensiones y contenidos de los ficheros
ficheros=
tribble(~contenido,~fichero,
"Hogar","691-Modulo64.zip",
"Vivienda","691-Modulo65.zip",
"DatosBasicos","691-Modulo66.zip",
"Nacimientos","691-Modulo67.zip",
"Embarazos","691-Modulo69.zip",
"Nupcialidad","691-Modulo71.zip"
)
# Bajamos ahora a la carpeta data aprovechando las
# funciones del paquete purrr del tidyverse
ficheros %>% pwalk(~download.file(paste0(urlbase,.y), destfile=paste0("data/",.x,".zip"),mode="wb"))Estrategia de análisis y lectura de datos
En la aplicación de ejemplo, vamos estudiar los factores asociados al uso de anticonceptivos modernos mediante regresión logística. La variable principal es si en la actualidad está utilizando métodos anticonceptivos modernos, que se corresponde con que la variable de uso actual (V313) tome el valor \(3\) (V313==3). Nos interesan las siguientes variables:
## Generalmente esto lo ocultaríamos con echo=FALSE,
# lo mostramos para garantizar la replicabilidad
varsel=tribble(~Variable,~Significado,~Módulo,~Fichero,~Tipo,
"V313","Uso actual por tipo de método","RE223132","Nacimientos","Explicada",
"V201","Total de niños nacidos","RE223132","Nacimientos","Explicativa",
"V218","Número de niños vivos","RE223132","Nacimientos","Explicativa",
"V106","Nivel educativo más alto","REC0111","DatosBasicos","Explicativa",
"V190","Índice de riqueza","REC0111","DatosBasicos","Explicativa",
"V012","Edad actual - Entrevistada","REC0111","DatosBasicos","Explicativa",
"V013","Edad actual por grupos de 5 años","REC0111","DatosBasicos","Explicativa",
"V025","Tipo de lugar de residencia","REC0111","DatosBasicos","Explicativa",
"V024","Región","REC0111","DatosBasicos","Explicativa",
"V501","Estado civil actual","RE516171","Nupcialidad","Explicativa",
"V213","Actualmente embarazada","RE223132","Nacimientos","Selección",
"V536","Actividad sexual reciente","RE516171","Nupcialidad","Selección",
"CASEID","Identificación individual","Todos","Todos","Técnica",
"HHID","Identificación hogar","Todos","Todos","Técnica",
"V022","Estratos","REC0111","DatosBasicos","Técnica",
"V001","Conglomerado","REC0111","DatosBasicos","Técnica",
"NCONGLOME","Número del conglomerado","REC0111","DatosBasicos","Técnica",
"V004","Unidad de área final","REC0111","DatosBasicos","Técnica",
"V005","Factor de ponderación (Mujer)","REC0111","DatosBasicos","Técnica")
varsel %>% knitr::kable()| Variable | Significado | Módulo | Fichero | Tipo |
|---|---|---|---|---|
| V313 | Uso actual por tipo de método | RE223132 | Nacimientos | Explicada |
| V201 | Total de niños nacidos | RE223132 | Nacimientos | Explicativa |
| V218 | Número de niños vivos | RE223132 | Nacimientos | Explicativa |
| V106 | Nivel educativo más alto | REC0111 | DatosBasicos | Explicativa |
| V190 | Índice de riqueza | REC0111 | DatosBasicos | Explicativa |
| V012 | Edad actual - Entrevistada | REC0111 | DatosBasicos | Explicativa |
| V013 | Edad actual por grupos de 5 años | REC0111 | DatosBasicos | Explicativa |
| V025 | Tipo de lugar de residencia | REC0111 | DatosBasicos | Explicativa |
| V024 | Región | REC0111 | DatosBasicos | Explicativa |
| V501 | Estado civil actual | RE516171 | Nupcialidad | Explicativa |
| V213 | Actualmente embarazada | RE223132 | Nacimientos | Selección |
| V536 | Actividad sexual reciente | RE516171 | Nupcialidad | Selección |
| CASEID | Identificación individual | Todos | Todos | Técnica |
| HHID | Identificación hogar | Todos | Todos | Técnica |
| V022 | Estratos | REC0111 | DatosBasicos | Técnica |
| V001 | Conglomerado | REC0111 | DatosBasicos | Técnica |
| NCONGLOME | Número del conglomerado | REC0111 | DatosBasicos | Técnica |
| V004 | Unidad de área final | REC0111 | DatosBasicos | Técnica |
| V005 | Factor de ponderación (Mujer) | REC0111 | DatosBasicos | Técnica |
Se trata de ficheros comprimidos que incluyen documentación en formato pdf además de los datos en formato sav de SPSS. Para extraer los sav que nos interesan procedemos como en la primera sesión: descomprimimos los ficheros que incluyen lo que nos interesa en una carpeta temporal y copiamos de vuelta a data/ los ficheros sav que nos interesan así como el cuestionario.
temp <- tempfile()
varsel %>% filter(Fichero!="Todos") %>% pull(Fichero) %>% unique() %>% paste0("data/",.,".zip") %>% walk(unzip,junkpaths=TRUE,exdir=temp) #Descomprime todo en la misma carpeta
.m = varsel %>% filter(Módulo!="Todos") %>% pull(Módulo) %>% unique()
dir(temp) %>% as_tibble() %>% filter(value %>% str_detect(.m %>% paste(collapse="|"))) %>% mutate(value=paste(temp,value,sep="/")) %>%
pull(value) %>%
file.copy("data/")
rm(.m)En cuanto a la estrategia de lectura, vamos a utilizar haven, que ya tenemos cargado, para leer primero el fichero REC0111 que es el que tiene más variables de las que me interesan.
DHS=read_sav("data/REC0111.sav")Hemos leído el objeto que tiene 38335 observaciones (mujeres) y 105 variables. Vamos a comprobar que la lectura ha ido bien utilizando las tablas de frecuencias de sjmisc::frq
DHS %>% frq(V025)##
## Tipo de lugar de residencia (V025) <numeric>
## # total N=38335 valid N=38335 mean=1.28 sd=0.45
##
## Value | Label | N | Raw % | Valid % | Cum. %
## -------------------------------------------------
## 1 | Urbano | 27427 | 71.55 | 71.55 | 71.55
## 2 | Rural | 10908 | 28.45 | 28.45 | 100.00
## <NA> | <NA> | 0 | 0.00 | <NA> | <NA>A diferencia de la semana pasada, en la que nuestro objetivo era trabajar con funciones “generales” de R, y preferíamos convertir las variables a factores “normales”, hoy vamos a utilizar casi exclusivamente paquetes diseñados para trabajar con encuestas y que entienden los formatos de etiqueta. Por ello es más cómodo hacerlo así. Se podrían cambiar los nombres de los variables, eso sí. En cualquier caso ayuda que funciones como sjmisc::find_var buscan tanto en los nombres de las variables como en las descripciones.
DHS %>% find_var("Edad") %>% kable()| col.nr | var.name | var.label |
|---|---|---|
| 13 | V012 | Edad actual - entrevistada |
| 14 | V013 | Edad actual por grupos de 5 años |
| 75 | V152 | Edad del jefe del hogar |
Vamos ahora a pegar en este objeto base el resto de variables. Las variables compartidas son dos:
ID1: la variable de año de la encuesta, que no nos interesa porque en todos los casos en 2019, pero que es útil si nos interesa hacer un análisis pooled combinando observaciones de los distintos años.CASEID: el identificador de individuo.
DHS= DHS %>%
left_join(read_sav("data/RE223132.sav"),by=c("ID1","CASEID")) %>%
left_join(read_sav("data/RE516171.sav"),by=c("ID1","CASEID"))Con lo que el nuevo objeto tiene ahora 361 variables. Vamos a ver, por ejemplo, las variables relacionadas con uso de anticonceptivos
DHS %>% find_var("anticoncep") %>% kable| col.nr | var.name | var.label |
|---|---|---|
| 152 | V312 | Método anticonceptivo actual |
| 165 | V326 | Fuente para obtener el actual método anticonceptivo |
| 166 | V327 | Fuente para obtener el actual método anticonceptivo (Grupos) |
| 173 | V364 | Uso e intención de anticonceptivos |
| 250 | V305_03 | Alguna vez ha usado: Inyección anticonceptiva |
| 278 | QI302_05A | Uso inyección anticonceptiva de 1 mes |
| 279 | QI302_05B | Uso inyección anticonceptiva de 3 meses |
| 323 | V632 | Tomo la decisión para el uso de anticonceptivos |
| 331 | V634 | Su esposo/compañero sabe que está usando anticonceptivos |
Selección de datos y especificación del diseño de la encuesta.
La información técnica de la encuesta está contenida en la ficha técnica. Nos da información de que los estratos en la encuesta son los departamentos y el área de residencia con tres tipos: área sede (Lima y capitales departamentales), resto urbano y rural. En cuanto a los conglomerados hay 3254 seleccionados. Las variables correspondientes son V001 para el conglomerado y V022 para el estrato.
En la ficha técnica se explica de dónde salen los pesos de la encuesta. El factor básico de muestreo o peso de diseño surge como la inversa de la probabilidad de selección de la vivienda. A posteriori, este se ajusta por la no respuesta para llegar al peso definitivo. Pero este es un peso para viviendas. Se deriva un peso específico para mujeres que tiene en cuenta la distribución de mujeres 12-49 años dentro de las viviendas y que son los que nos interesan. Pero se especifica que debemos dividirlo por 1000000 para obtener el peso:
DHS = DHS %>% mutate(peso=V005/1000000)
DHS %>% descr(peso)Vemos que estos son pesos ya estandarizados con una media cercana a 1, algo que confirma la documentación de DHS (normalization of weights). Esto es importante, puesto que como expusimos en la sesión 1, nos permitirá estar más tranquilos a la hora de calcular errores estándar al evitar la confusión entre pesos de muestreo y pesos de frecuencia. Podéis ver una discusión al respecto aquí (Lumley, 2020)
Una vez tenemos definidos los pesos podemos obtener tablas de frecuencias como, por ejemplo,la de anticonceptivo de uso actual de este modo:
DHS %>% frq(V312,weights=peso,sort.frq="desc") %>% kable()
|
Vemos que la principal categoría son las mujeres que no están usando ningún método, y entre las que usan los principales métodos son la inyección, condón y abstinencia periódica. Nuestra variable explicada es el uso de métodos modernos que agrupa a las que usan métodos distintos de la abstinencia periódica y el retiro. Se trata del nivel 3 de la variable V313 que vamos a llamar moderno:
DHS %>% frq(V313,weights=peso)##
## Uso actual por tipo de método (V313) <numeric>
## # total N=33381 valid N=33381 mean=1.44 sd=1.41
##
## Value | Label | N | Raw % | Valid % | Cum. %
## -------------------------------------------------------------
## 0 | No hay método | 15916 | 47.68 | 47.68 | 47.68
## 1 | Método folclórico | 102 | 0.31 | 0.31 | 47.99
## 2 | Método tradicional | 4237 | 12.69 | 12.69 | 60.68
## 3 | Método moderno | 13126 | 39.32 | 39.32 | 100.00
## <NA> | <NA> | 0 | 0.00 | <NA> | <NA>DHS=DHS %>%
mutate(moderno=factor((V313==3), labels=c("No","Si")))
DHS %>% frq(moderno,weights=peso) %>% kable()
|
Por otro lado, nos interesa estudiar el uso de métodos modernos entre las mujeres que estén expuestas al riesgo de embarazo. Esto nos lleva a seleccionar las mujeres que no están embarazadas (V213) y que hayan tenido relaciones sexuales en el último mes (V536).
Como sugerimos en la primera sesión, si todo nuestro análisis va a ser realizado con un subconjunto de datos es más cómodo definir un objeto nuevo, que vamos a llamar mex (de mujeres expuestas) que será el que estudiemos a partir de ahora. En nuestros datos el sí se refiere a las mujeres que están usando métodos modernos, y el no engloba dos categorías: las que no usan ningún método y las que usan métodos tradicionales. Podríamos haber hecho el análisis condicional a uso, o analizado la variable con tres categorías con un modelo multinomial similar.
Una vez seleccionado el subconjunto estudiado definimos en mex las mujeres expuestas al riesgo de embarazo y, por tanto, las que pueden estar interesadas en el uso de métodos anticonceptivos modernos.
DHS %>% frq(V213,V536,weights=peso) %>% kable()
|
|
# Definición de variable de exposición
mex = DHS %>% filter(V213==0,V536==1)Ya tenemos preparado el objeto de datos que queremos analizar. Únicamente nos queda definir el diseño de la encuesta para su uso en el paquete survey. Como comentamos la filosofía de survey es tener un objeto que contiene la declaración del diseño de la encuesta de modo que posteriormente no sea necesario volverlo a hacer en cada modelo, algo parecido a lo que hemos preferido hacer con los datos. Esta fase es sencilla: hay que especificar las variables de estratos, conglomerados y pesos. Con eso, el paquete ya sabe lo que hay que hacer. En nuestro caso, y para poder estimar el efecto del diseño, vamos a crear tres objetos distintos: el que contiene el diseño correcto (mexdis), el que sólo utiliza pesos pero no diseño complejo (mexdis.w), y, por último, el que trata los datos como si fuesen un muestro aleatorio, sin utilizar pesos (mexdis.unw).
## Diseño de la encuesta para mujeres expuestas
# Recomendado por DHS, id: conglomerados, strata: estratos.
mexdis=svydesign(ids=~V021,
strata=~V022,
weights=~peso,
data=mex,
nest=TRUE)
## Para comparación y cálculo del efecto del muestreo.
# Diseño sin uso de pesos
mexdis.unw=svydesign(ids=~1,strata=NULL,weights=NULL,
data=mex)
# Diseño con pesos pero sin estructura compleja
mexdis.w=svydesign(ids=~1,strata=NULL,weights=~peso,
data=mex)La función genérica summary nos proporciona para un objeto survey.design las características del muestreo, incluyendo el número de mujeres y de conglomerados dentro de cada uno de los estratos, así como el nombre de las variables. Nótese que nos dice que se trata de un muestreo con reemplazo. No lo es: cada mujer sólo puede entrevistarse una vez. Tenerlo en cuenta complica notablemente los cálculos sin impacto en los resultados cuando la población de cada conglomerado es relativamente grade. Es posible exigir una corrección de población finita, pero para ello es necesario especificar el argumento fpc proporcionando la población en cada estrato, algo que no nos proporciona la DHS.
summary(mexdis)## Stratified 1 - level Cluster Sampling design (with replacement)
## With (3242) clusters.
## svydesign(ids = ~V021, strata = ~V022, weights = ~peso, data = mex,
## nest = TRUE)
## Probabilities:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.04204 1.08534 2.21995 3.46643 4.01508 94.97578
## Stratum Sizes:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
## obs 8 27 21 23 30 44 47 79 56 486 17 60 93 67 58 22 15 23 30 16 283 19
## design.PSU 2 4 4 4 5 8 8 13 9 62 4 11 15 12 10 4 2 5 4 3 40 4
## actual.PSU 2 4 4 4 5 8 8 13 9 62 4 11 15 12 10 4 2 5 4 3 40 4
## 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
## obs 22 34 20 17 52 42 22 34 364 56 84 139 116 76 13 34 15 41 48 73 21 45
## design.PSU 4 6 4 3 8 6 5 7 58 11 18 23 20 15 3 6 2 6 8 12 4 8
## actual.PSU 4 6 4 3 8 6 5 7 58 11 18 23 20 15 3 6 2 6 8 12 4 8
## 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
## obs 55 91 35 32 31 39 15 63 331 8 19 16 34 24 32 22 39 17 17 524 88 146
## design.PSU 9 13 6 6 6 8 3 12 50 2 3 2 5 4 6 4 7 3 3 66 19 27
## actual.PSU 9 13 6 6 6 8 3 12 50 2 3 2 5 4 6 4 7 3 3 66 19 27
## 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
## obs 256 109 144 14 38 39 53 36 15 33 28 29 20 319 9 31 23 10 15 30 32
## design.PSU 42 17 24 4 8 9 9 6 3 6 6 6 4 46 2 6 4 2 4 5 7
## actual.PSU 42 17 24 4 8 9 9 6 3 6 6 6 4 46 2 6 4 2 4 5 7
## 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106
## obs 401 38 28 30 43 19 30 73 19 7 443 26 42 54 37 21 12 74 138
## design.PSU 73 6 7 7 7 4 5 11 4 2 64 5 8 9 7 4 2 12 22
## actual.PSU 73 6 7 7 7 4 5 11 4 2 64 5 8 9 7 4 2 12 22
## 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123
## obs 136 157 75 12 49 43 61 29 75 39 16 80 64 244 34 42 79
## design.PSU 22 24 9 2 8 9 10 5 13 6 3 16 10 32 7 11 14
## actual.PSU 22 24 9 2 8 9 10 5 13 6 3 16 10 32 7 11 14
## 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
## obs 90 60 35 21 71 27 82 184 43 51 98 87 74 27 59 55 104
## design.PSU 13 9 7 3 11 6 13 23 8 10 16 13 9 4 10 8 16
## actual.PSU 13 9 7 3 11 6 13 23 8 10 16 13 9 4 10 8 16
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157
## obs 66 121 133 228 537 551 421 77 187 172 21 95 170 47 125 76 78
## design.PSU 12 16 42 49 93 92 63 15 30 28 4 14 25 8 17 11 11
## actual.PSU 12 16 42 49 93 92 63 15 30 28 4 14 25 8 17 11 11
## 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174
## obs 44 12 34 18 57 44 281 62 116 149 60 37 72 192 43 36 66
## design.PSU 8 2 4 3 10 8 32 12 24 25 9 7 14 26 8 9 13
## actual.PSU 8 2 4 3 10 8 32 12 24 25 9 7 14 26 8 9 13
## 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191
## obs 33 82 48 116 35 81 76 40 48 30 19 19 115 61 57 186 12
## design.PSU 6 15 9 23 6 14 17 8 10 6 5 3 22 11 13 34 2
## actual.PSU 6 15 9 23 6 14 17 8 10 6 5 3 22 11 13 34 2
## 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208
## obs 17 36 55 61 54 107 135 163 184 4 15 14 14 15 10 47 61
## design.PSU 4 6 8 10 8 16 18 25 22 2 3 2 2 2 2 8 10
## actual.PSU 4 6 8 10 8 16 18 25 22 2 3 2 2 2 2 8 10
## 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225
## obs 39 72 273 11 59 74 55 40 35 52 42 67 143 287 66 114 171
## design.PSU 7 13 52 2 9 10 6 6 5 8 6 11 19 32 14 27 35
## actual.PSU 7 13 52 2 9 10 6 6 5 8 6 11 19 32 14 27 35
## 226 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243
## obs 154 11 67 19 64 97 117 66 35 49 93 161 50 83 24 64 138
## design.PSU 28 5 14 4 9 16 18 10 5 7 14 23 9 10 4 9 19
## actual.PSU 28 5 14 4 9 16 18 10 5 7 14 23 9 10 4 9 19
## 244 245 246 247 248 249
## obs 155 227 15 44 41 183
## design.PSU 21 29 2 6 7 22
## actual.PSU 21 29 2 6 7 22
## Data variables:
## [1] "ID1" "HHID" "CASEID" "V001" "V002"
## [6] "V003" "V004" "V007" "V008" "V009"
## [11] "V010" "V011" "V012" "V013" "V014"
## [16] "V015" "V017" "V018" "V019" "V019A"
## [21] "V020" "V021" "V023" "V024" "V025"
## [26] "V026" "V027" "V028" "V029" "V030"
## [31] "V031" "V032" "V033" "V034" "V040"
## [36] "V042" "V043" "V044" "V000" "Q105DD"
## [41] "V101" "V102" "V103" "V104" "V105"
## [46] "V106" "V107" "V113" "V115" "V116"
## [51] "V119" "V120" "V121" "V122" "V123"
## [56] "V124" "V125" "V127" "V128" "V129"
## [61] "V130" "V131" "V133" "V134" "V135"
## [66] "V136" "V137" "V138" "V139" "V140"
## [71] "V141" "V149" "V150" "V151" "V152"
## [76] "V153" "AWFACTT" "AWFACTU" "AWFACTR" "AWFACTE"
## [81] "AWFACTW" "V155" "V156" "V157" "V158"
## [86] "V159" "V160" "V161" "V166" "V167"
## [91] "V168" "ML101" "QD333_1" "QD333_2" "QD333_3"
## [96] "QD333_4" "QD333_5" "QD333_6" "UBIGEO" "V022"
## [101] "V005" "V190" "V191" "mujeres12a49" "NCONGLOME"
## [106] "V201" "V202" "V203" "V204" "V205"
## [111] "V206" "V207" "V208" "V209" "V210"
## [116] "V211" "V212" "V213" "V214" "V215"
## [121] "V216" "V217" "V218" "V219" "V220"
## [126] "V221" "V222" "V223" "V224" "V225"
## [131] "V226" "V227" "V228" "V229" "V230"
## [136] "V231" "V232" "V233" "V234" "V235"
## [141] "V237" "V238" "V239" "V240" "V241"
## [146] "V242" "V243" "V301" "V302" "V310"
## [151] "V311" "V312" "V313" "V315" "V316"
## [156] "V317" "V318" "V319" "V320" "V321"
## [161] "V322" "V323" "V323A" "V325A" "V326"
## [166] "V327" "V337" "V359" "V360" "V361"
## [171] "V362" "V363" "V364" "V367" "V372"
## [176] "V372A" "V375A" "V376" "V376A" "V379"
## [181] "V380" "V384A" "V384B" "V384C" "V393"
## [186] "V394" "V395" "V3A00A" "V3A00B" "V3A00C"
## [191] "V3A00D" "V3A00E" "V3A00F" "V3A00G" "V3A00H"
## [196] "V3A00I" "V3A00J" "V3A00K" "V3A00L" "V3A00M"
## [201] "V3A00N" "V3A00O" "V3A00P" "V3A00Q" "V3A00R"
## [206] "V3A00S" "V3A00T" "V3A00U" "V3A00V" "V3A00W"
## [211] "V3A00X" "V3A00Y" "V3A00Z" "V3A01" "V3A02"
## [216] "V3A03" "V3A04" "V3A05" "V3A06" "V3A07"
## [221] "V3A08A" "V3A08B" "V3A08C" "V3A08D" "V3A08E"
## [226] "V3A08F" "V3A08G" "V3A08H" "V3A08I" "V3A08J"
## [231] "V3A08K" "V3A08L" "V3A08M" "V3A08N" "V3A08O"
## [236] "V3A08P" "V3A08Q" "V3A08R" "V3A08S" "V3A08T"
## [241] "V3A08U" "V3A08V" "V3A08W" "V3A08X" "V3A08Z"
## [246] "V3A09A" "V3A09B" "V305_01" "V305_02" "V305_03"
## [251] "V305_05" "V305_06" "V305_07" "V305_08" "V305_09"
## [256] "V305_10" "V305_11" "V305_13" "V305_14" "V305_15"
## [261] "V305_16" "V307_01" "V307_02" "V307_03" "V307_04"
## [266] "V307_05" "V307_06" "V307_07" "V307_08" "V307_09"
## [271] "V307_10" "V307_11" "V307_12" "V307_13" "V307_14"
## [276] "V307_15" "V307_16" "QI302_05A" "QI302_05B" "V501"
## [281] "V502" "V503" "V504" "V505" "V506"
## [286] "V507" "V508" "V509" "V510" "V511"
## [291] "V512" "V513" "V525" "V527" "V528"
## [296] "V529" "V530" "V531" "V532" "V535"
## [301] "V536" "V537" "V538" "V539" "V540"
## [306] "V541" "V602" "V603" "V604" "V605"
## [311] "V613" "V614" "V616" "V621" "V623"
## [316] "V624" "V625" "V626" "V627" "V628"
## [321] "V629" "V631" "V632" "V633A" "V633B"
## [326] "V633C" "V633D" "V633E" "V633F" "V633G"
## [331] "V634" "V701" "V702" "V704" "V705"
## [336] "V714" "V714A" "V715" "V716" "V717"
## [341] "V719" "V721" "V729" "V730" "V731"
## [346] "V732" "V739" "V740" "V741" "V743A"
## [351] "V743B" "V743C" "V743D" "V743E" "V743F"
## [356] "V744A" "V744B" "V744C" "V744D" "V744E"
## [361] "V746" "peso" "moderno"Estimación de errores estándar: efecto del diseño
Una vez definido el objeto de la encuesta, podemos analizar los resultados con las distintas funciones del paquete survey que utilizan una sintaxis parecida: requieren especificar el objeto que incluye el diseño, y son funciones generalmente precedidas por svy. En esta sección nos vamos a concentrar en dos de las funciones más sencillas, svymean para el cálculo de medias y svyciprop para el cálculo de intervalos de la proporción, analizando la influencia del diseño muestral utilizado sobre los estimadores puntuales y sobre los errores estándar estimados.
Por ejemplo, para estimar un intervalo de confianza de la proporción de mujeres usando anticonceptivos modernos en la población podemos emplear:
svyciprop(~moderno,mexdis)## 2.5% 97.5%
## moderno 0.635 0.621 0.65svyciprop(~moderno,mexdis.w)## 2.5% 97.5%
## moderno 0.635 0.622 0.65svyciprop(~moderno,mexdis.unw)## 2.5% 97.5%
## moderno 0.663 0.657 0.67Esta tabla nos muestra las diferencias básicas entre usar el diseño complejo, usar pesos y no usarlos.
La estimación puntual es la misma para los análisis que usan el diseño complejo y los que usan sólo pesos. Serán algo distintos los errores estándar estimados.
Si no se utilizan pesos, los estimadores no serán representativos de la población, y, en general, se infravalorará la dispersión de los estimadores, llevando a intervalos de confianza excesivamente estrechos y a niveles de significación reales mayores que los nominales.
A la relación entre el error estándar que tiene en cuenta el diseño complejo de la encuesta y el que asume muestreo aleatorio simple le llamamos efecto del diseño. Este será específico para cada variable, y está ligado al coeficiente de correlación intraclase.
En este caso en concreto vemos que la diferencia entre la inferencia con estructura compleja y con pesos es mínima. La estimación de la proporción sin pesos es muy diferente, ni siquiera cae en el intervalo al 95% de los otros dos, y subestima la variabilidad al ser mucho más reducido
Como mínimo siempre debemos utilizar pesos con encuestas complejas.
Para calcular el efecto de diseño necesitamos calcular el efecto de diseño. Para ello tenemos que calcular los errores estándar de cada una de las proporciones. En survey se puede aplicar la función SE a los objetos que nos dan estimaciones y nos proporciona los errores estándar.
SE_comp=svyciprop(~moderno,mexdis) %>% SE()
SE_w= svyciprop(~moderno,mexdis.w) %>% SE()
SE_unw= svyciprop(~moderno,mexdis.unw) %>% SE()
tribble(~Método,~SE,
"Encuesta Compleja",SE_comp,
"Con pesos",SE_w,
"Sin pesos",SE_unw) %>%
mutate(var=SE^2,deff=first(var)/var) %>% kable(digits=6)| Método | SE | var | deff |
|---|---|---|---|
| Encuesta Compleja | 0.007050 | 5.0e-05 | 1.000000 |
| Con pesos | 0.006738 | 4.5e-05 | 1.094999 |
| Sin pesos | 0.003356 | 1.1e-05 | 4.413595 |
Vemos como el efecto de diseño es bastante importante cuando no se utilizan pesos, en este caso la varianza se subestima 4.4 veces. Por cierto: la función svyciprop proporciona una estimación de la proporción basada en un modelo logit que funciona mejor cuando la proporción está cerca de 0 o de 1.
Si queremos obtener una tabla que recoja la estimación y el error estándar tanto para la proporción (cuando la variable es un factor), o para la media (cuando la variable es cuantitativa), tenemos la función svymean. Podemos separar las variables con +. En este caso tabulamos la proporción de usuarios y no usuarios de métodos modernos, y el número medio de hijos nacidos vivos (V201). El único problema es que no utiliza etiquetas. Una ventaja es que también le podemos solicitar que calcule el efecto de diseño. Si bien el exacto no es calculable salvo que los pesos sean la inversa de la probabilidad de selección, que no es el caso, especificando el valor (replace) lo calcula basándose en las fórmulas válidas con reemplazo. Podemos ver que el efecto de diseño es razonablemente parecido al que habíamos obtenido anteriormente para la proporción. La diferencia se debe al distinto método de cálculo de la varianza.
svymean(~moderno+V201,mexdis,deff="replace")## mean SE DEff
## modernoNo 0.3650127 0.0070504 4.2517
## modernoSi 0.6349873 0.0070504 4.2517
## V201 2.0602905 0.0212410 3.3878Podemos también comprobar que el efecto de diseño que calcula coincide aproximadamente con la comparación en base a los diseños sin pesos, y ver qué ocurre con el diseño sólo con pesos.
svymean(~moderno+V201,mexdis.unw)## mean SE
## modernoNo 0.33653 0.0034
## modernoSi 0.66347 0.0034
## V201 2.35670 0.0117svymean(~moderno+V201,mexdis.w)## mean SE
## modernoNo 0.36501 0.0067
## modernoSi 0.63499 0.0067
## V201 2.06029 0.0215De nuevo las diferencias son relativamente pequeñas para los estimadores que utilizan pesos y grandes para los que no los utilizan.
En ocasiones puede ser útil obtener algunos de los resultados asociados a svymean. Existen varias funciones para extraer resultados o reelaborarlos en forma de intervalos de confianza.
medias=svymean(~moderno+V201,mexdis,deff="replace")
coef(medias)## modernoNo modernoSi V201
## 0.3650127 0.6349873 2.0602905SE(medias)## modernoNo modernoSi V201
## 0.007050363 0.007050363 0.021240995confint(medias,level=.9)## 5 % 95 %
## modernoNo 0.3534159 0.3766095
## modernoSi 0.6233905 0.6465841
## V201 2.0253522 2.0952289Vamos ahora a comprobar las formas alternativas de calcular estas cantidades con las funciones que vimos la semana pasada. No integran la posibilidad de análisis complejo pero si la del uso de pesos. La estimación de proporciones, por ejemplo, la podemos hacer con un modelo de probabilidad lineal. Sabemos que estos presentan heterocedasticidad asi que utilizamos el estimador del SE de White consistente con heterocedasticidad.
Estas son las estimaciones de la proporción y de la media en base al modelo lineal utilizando pesos:
m.prop=lm(I(moderno=="Si")~1,weights=mex$peso,data=mex)
coeftest(m.prop,vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6349873 0.0067403 94.207 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m.mean=lm(V201~1,weights=mex$peso,data=mex)
coeftest(m.mean,vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.060291 0.021459 96.012 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Vemos que, efectivamente coinciden con los valores calculados por svymean con mexdis.w tanto en la estimación puntual como en el error estándar estimado.
Una alternativa intermedia entre el uso de pesos y los estimadores de survey es la utilización del estimador de varianzas-covarianzas robusto por conglomerados, vcovCL. Este también tiene en cuenta la variabilidad de los conglomerados, únicamente no tiene en cuenta los estratos.
coeftest(m.prop,vcov=vcovCL,cluster=mex$V021)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6349873 0.0071863 88.361 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coeftest(m.mean,vcov=vcovCL,cluster=mex$V021)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.060291 0.024613 83.709 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Vemos que en este caso los SE son incluso algo mayores a los calculados por survey por lo que no tenemos el problema de la infravaloración. Las diferencias son muy pequeñas y esto nos ofrece una alternativa a survey para calcular de manera correcta la precisión de nuestro estimadores.
Veamos que ocurre cuando no utilizamos pesos:
m.prop0=lm(I(moderno=="Si")~1,data=mex)
coeftest(m.prop0,vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.663472 0.003356 197.7 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m.mean0=lm(V201~1,data=mex)
coeftest(m.mean0,vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.356703 0.011699 201.45 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Vemos que los valores coinciden con los que había calculado survey con el diseño sin pesos. El estimador puntual bastante diferente nos recuerda que no debemos de estimar sin pesos en encuestas complejas. Este método, que es aún habitual en economía, es el único claramente incorrecto. En nuestro caso concreto hemos visto que siempre que al menos utilicemos pesos no hay grandes diferencias entre estimaciones de errores estándar.
Tabulaciones cruzadas de encuestas complejas
Como hemos visto la función svymean nos sirve para calcular proporciones para una variable o medias. Existen varios modos de obtener tabulaciones cruzadas con survey. La primera es la función svytable. El único inconveniente es que no utiliza las etiquetas lo que dificulta la interpretación.
Vamos a crear la tabla de uso de anticonceptivos modernos por estado civil. Se utiliza la función svytable, de la que si pedimos un summary nos proporciona un test de independencia de la F por defecto (también se puede hacer chi cuadrado) que tiene en cuenta la estructura compleja.
tabEC=svytable(~moderno+V501,design=mexdis)
tabEC## V501
## moderno 0 1 2 3 4 5
## No 615.141868 2128.950681 3516.030972 3.803241 17.851755 332.262810
## Si 1402.535887 3238.787541 6145.075941 3.935345 22.509699 693.143262summary(tabEC,statistic="Chisq")## V501
## moderno 0 1 2 3 4 5
## No 615 2129 3516 4 18 332
## Si 1403 3239 6145 4 23 693
##
## Pearson's X^2: Rao & Scott adjustment
##
## data: svychisq(~moderno + V501, design = mexdis, statistic = "Chisq")
## X-squared = 69.664, df = 5, p-value = 0.007007Podemos comparar estos resultados con las tablas que produce sjPlot::tab_xtab que produce un resultado análogo aunque no exactamente igual (por la estructura compleja) pero con etiquetas y formato html
with(mex,sjPlot::tab_xtab(moderno,V501,weight=peso))| moderno | Estado civil actual | Total | |||||
|---|---|---|---|---|---|---|---|
| Nunca casada | Casado | Viviendo juntos | Viuda | Divorciada | No viven juntos | ||
| No | 615 | 2129 | 3516 | 4 | 18 | 332 | 6614 |
| Si | 1403 | 3239 | 6145 | 4 | 23 | 693 | 11507 |
| Total | 2018 | 5368 | 9661 | 8 | 41 | 1025 | 18121 | χ2=63.855 · df=5 · Cramer’s V=0.059 · Fisher’s p=0.000 |
Otra función útil para producir tablas es svyby. Esto equivale a declarar grupos en el tidyverse, y se utiliza en combinación con las funciones elementales como svymean. Vamos a ver la misma tabla hecha de este modo, que nos proporciona las proporciones usando y no usando en cada categoria. Una ventaja de svyby es que el resultado es un data.frame, que se puede convertir en tibble y manipular. Además, preserva las etiquetas de las variables labelled, que se puede hacer con as_factor que muestren las etiquetas
svyby(~moderno,~V501,mexdis,svymean) %>% as_tibble() %>%
mutate(V501=as_factor(V501)) %>% kable(digits=3)| V501 | modernoNo | modernoSi | se.modernoNo | se.modernoSi |
|---|---|---|---|---|
| Nunca casada | 0.305 | 0.695 | 0.023 | 0.023 |
| Casado | 0.397 | 0.603 | 0.013 | 0.013 |
| Viviendo juntos | 0.364 | 0.636 | 0.009 | 0.009 |
| Viuda | 0.491 | 0.509 | 0.197 | 0.197 |
| Divorciada | 0.442 | 0.558 | 0.183 | 0.183 |
| No viven juntos | 0.324 | 0.676 | 0.032 | 0.032 |
Vemos como las proporciones de uso son mayores en las nunca casadas y las que no viven juntos, presumiblemente porque quieren evitar el embarazo. Los SE nos dan una idea de la incertidumbre, siendo estimaciones poco precisas las de categorías como viudas o divorciadas.
Si preferimos, se pueden obtener los intervalos que también se pueden manipular como tibbles
svyby(~moderno,~V501,mexdis,svymean) %>% confint()## 2.5 % 97.5 %
## 0:modernoNo 0.25976429 0.3499880
## 1:modernoNo 0.37034333 0.4228961
## 2:modernoNo 0.34642279 0.3814505
## 3:modernoNo 0.10450791 0.8784213
## 4:modernoNo 0.08348453 0.8011097
## 5:modernoNo 0.26158863 0.3864723
## 0:modernoSi 0.65001195 0.7402357
## 1:modernoSi 0.57710394 0.6296567
## 2:modernoSi 0.61854947 0.6535772
## 3:modernoSi 0.12157873 0.8954921
## 4:modernoSi 0.19889028 0.9165155
## 5:modernoSi 0.61352769 0.7384114La tabla se puede ampliar a más de dos variables: por ejemplo, nivel educativo. También podemos pedir intervalos en vez de SE. Esto nos muestra, en este caso, las limitaciones de esta función para calcular intervalos que se pueden ir del intervalo 0,1. Las que se obtienen con svyciprop no tienen esos problemas.
svyby(~moderno,~V106+V501,mexdis,svymean,vartype="ci") %>% as_tibble() %>%
mutate(across(starts_with("V"),as_factor)) %>% kable(digits=3)| V106 | V501 | modernoNo | modernoSi | ci_l.modernoNo | ci_l.modernoSi | ci_u.modernoNo | ci_u.modernoSi |
|---|---|---|---|---|---|---|---|
| Sin educación | Nunca casada | 0.675 | 0.325 | 0.146 | -0.204 | 1.204 | 0.854 |
| Primario | Nunca casada | 0.478 | 0.522 | 0.274 | 0.318 | 0.682 | 0.726 |
| Secundario | Nunca casada | 0.263 | 0.737 | 0.191 | 0.666 | 0.334 | 0.809 |
| Mayor | Nunca casada | 0.313 | 0.687 | 0.258 | 0.632 | 0.368 | 0.742 |
| Sin educación | Casado | 0.646 | 0.354 | 0.546 | 0.253 | 0.747 | 0.454 |
| Primario | Casado | 0.537 | 0.463 | 0.494 | 0.420 | 0.580 | 0.506 |
| Secundario | Casado | 0.365 | 0.635 | 0.325 | 0.596 | 0.404 | 0.675 |
| Mayor | Casado | 0.345 | 0.655 | 0.302 | 0.612 | 0.388 | 0.698 |
| Sin educación | Viviendo juntos | 0.532 | 0.468 | 0.431 | 0.368 | 0.632 | 0.569 |
| Primario | Viviendo juntos | 0.451 | 0.549 | 0.417 | 0.515 | 0.485 | 0.583 |
| Secundario | Viviendo juntos | 0.333 | 0.667 | 0.309 | 0.644 | 0.356 | 0.691 |
| Mayor | Viviendo juntos | 0.344 | 0.656 | 0.309 | 0.620 | 0.380 | 0.691 |
| Sin educación | Viuda | 0.000 | 1.000 | 0.000 | 1.000 | 0.000 | 1.000 |
| Primario | Viuda | 0.607 | 0.393 | 0.034 | -0.180 | 1.180 | 0.966 |
| Secundario | Viuda | 0.549 | 0.451 | -0.049 | -0.148 | 1.148 | 1.049 |
| Mayor | Viuda | 1.000 | 0.000 | 1.000 | 0.000 | 1.000 | 0.000 |
| Primario | Divorciada | 1.000 | 0.000 | 1.000 | 0.000 | 1.000 | 0.000 |
| Secundario | Divorciada | 0.068 | 0.932 | -0.078 | 0.786 | 0.214 | 1.078 |
| Mayor | Divorciada | 0.665 | 0.335 | 0.390 | 0.061 | 0.939 | 0.610 |
| Sin educación | No viven juntos | 0.323 | 0.677 | -0.089 | 0.265 | 0.735 | 1.089 |
| Primario | No viven juntos | 0.316 | 0.684 | 0.177 | 0.544 | 0.456 | 0.823 |
| Secundario | No viven juntos | 0.355 | 0.645 | 0.266 | 0.557 | 0.443 | 0.734 |
| Mayor | No viven juntos | 0.295 | 0.705 | 0.195 | 0.605 | 0.395 | 0.805 |
En este caso hemos transformado las dos variables en factores con etiquetas con el uso de across y as_factor.
Estimación de modelos con survey
Existen en el paquete survey funciones específicas para los métodos estadísticos más usuales como modelos lineales generalizados -que incluyen los lineales, logit, poisson, …- (svyglm), análisis de supervivencia (svysurvreg y svycoxph) o variables instrumentales (svyivreg). También existen dos extensiones del mismo autor, svyVGAM para el cálculo de modelos VGAM, y svy2lme para modelos lineales mixtos (multinivel). Existe una viñeta de svyVGAM que explica cómo se puede extender survey para la estimación de todo tipo de modelos especificando la función de verosimilitud, como por ejemplo la regresión de hurdle o inflada en cero. Por otro lado, como hemos visto, una alternativa es la utilización de las funciones especificas. svyVGAM, en cualquier caso, es muy interesante puesto que los modelos VGAM es una clase muy general de modelos que extiende los modelos GAM en varias direcciones (modelos vectoriales para cópulas, modelización conjunta de esperanas y varianzas para muchas familias de distribuciones, …). Lumley y Scott (2017) presenta una revisión útil incluyendo un ejemplo sobre regresión con diseño complejo. El paquete svydiags tiene funciones para diagnóstico de modelos estimados con survey
Por otro lado, como alternativa hemos visto que también se pueden utilizar los paquetes de modelización de R con datos de encuesta. Lo más importante es no olvidarnos de utilizar pesos, y que estos pesos estén estandarizados de modo que tengan media 1. Para tener en cuenta la estructura del diseño se puede emplear vcovCL como estimador de la matriz de varianzas-covarianzas.
Vamos a ver un ejemplo de aplicación: la estimación de un modelo para la proporción de usuarios de métodos modernos. Vamos a estimar un modelo logit, pero adelantándonos a uno de los problemas típicos de estos modelos, la sobredispersión, especificamos como familia la quasibinomial, que permite la estimación de un parámetro de sobredispersión que ajusta la variabilidad respecto a la fórmula de la binomial. En caso de que el parámetro sea cercano a uno, siempre podemos estimar el modelo logit estándar.
Antes de hacer la estimación, y puesto que las funciones de regresión y demás no nos dan los resultados con etiquetas, tenemos que transformar a variables con etiquetas a partir del objeto mex, y crear un nuevo objeto de diseño. De nuevo empleamos as_factor en combinación con across y una de las funciones auxiliares de selección, one_of. Otra alternativa sería transformar variables dentro del objeto de diseño, pero eso es más complejo a no ser que utilicemos el paquete srvyr, que permite operar con mutate dentro de objetos de diseño (luego lo mostraremos).
En el primer modelo introducimos la edad por grupos de edad para identificar si la relación es más o menos lineal.
mexdis2 = mex %>% mutate(across(one_of("V013","V501","V106","V190",
"V025"),as_factor)) %>%
svydesign(ids=~V021,
strata=~V022,
weights=~peso,
data=.,
nest=TRUE)
mod1 = svyglm(moderno~V013+V501+V106+V190+V025+V218,design=mexdis2,
family=quasibinomial())
summary(mod1)##
## Call:
## svyglm(formula = moderno ~ V013 + V501 + V106 + V190 + V025 +
## V218, design = mexdis2, family = quasibinomial())
##
## Survey design:
## svydesign(ids = ~V021, strata = ~V022, weights = ~peso, data = .,
## nest = TRUE)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.400117 1.362786 -0.294 0.769083
## V013De 15 a 19 años de edad 0.133970 1.362273 0.098 0.921667
## V013De 20 a 24 años de edad 0.272805 1.359246 0.201 0.840944
## V013De 25 a 29 años de edad -0.126293 1.358204 -0.093 0.925921
## V013De 30 a 34 años de edad -0.613522 1.358727 -0.452 0.651632
## V013De 35 a 39 años de edad -0.589294 1.359962 -0.433 0.664816
## V013De 40 a 44 años de edad -1.108058 1.359591 -0.815 0.415141
## V013De 45 a 49 años de edad -2.006914 1.361574 -1.474 0.140597
## V501Casado -0.002925 0.139848 -0.021 0.983312
## V501Viviendo juntos -0.015449 0.133127 -0.116 0.907622
## V501Viuda 0.308249 0.931452 0.331 0.740718
## V501Divorciada 0.169157 1.018097 0.166 0.868050
## V501No viven juntos 0.134473 0.192015 0.700 0.483778
## V106Primario 0.208748 0.179219 1.165 0.244208
## V106Secundario 0.613878 0.183341 3.348 0.000823 ***
## V106Mayor 0.569735 0.193999 2.937 0.003342 **
## V190Pobrer 0.368877 0.079797 4.623 3.95e-06 ***
## V190Medio 0.495442 0.107518 4.608 4.24e-06 ***
## V190Rico 0.697926 0.117240 5.953 2.94e-09 ***
## V190Más rico 0.996429 0.141409 7.046 2.27e-12 ***
## V025Rural -0.079816 0.079773 -1.001 0.317133
## V218 0.310764 0.027366 11.356 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 1.001175)
##
## Number of Fisher Scoring iterations: 4Anova(mod1)Vemos que el parámetro de dispersión es muy cercano a 1. Eso quiere decir que no es necesario utilizar la familia quasibinomial(). Respecto a las variables, la edad es muy significativa, pero no se aprecia que sea lineal: alcanza un pico a los 20-24 años. Podemos utilizar funciones suaves tipo spline para modelizarla. En estado civil, la categoría de referencia omitida, las nunca casadas, serían las que menos utilizan, a diferencia de las tabulaciones iniciales. Sin embargo, la variable no es significativa. Otra posibilidad es recodificar la variable en sólo 2 categorías, nunca casadas y alguna vez casadas, o en 3, nunca unidas, unidas y alguna vez casadas. La variable de educación sí que es altamente significativa siendo las mujeres que más usan las de educación secundaria. La riqueza también es importante usando más cuanto más ricas. En medio rural se usa menos pero la diferencia no es significativa. La variable de número de hijos vivos es altamente significativa, usando más las mujeres con más hijos, aunque tampoco sabemos si la relación es lineal.
Estimamos un segundo modelo permitiendo un spline cúbico para captar el efecto de la edad, y para la variable número de hijos utilizamos un factor con la categoría abierta de 4 y más. Quitamos la variable de lugar de residencia. Lo mostramos con tab_model, que saca una tabla publicable (salvo por lo que hemos hecho con el número de hijos).
library(splines)
mexdis2=mexdis2 %>% as_survey %>% mutate(HVivos=factor(pmin(V218,4)))
mod2 = svyglm(moderno~ns(V012,4)+V501+V106+V190+HVivos,design=mexdis2, family=quasibinomial())
tab_model(mod2,show.reflvl=TRUE,show.aic=TRUE,show.aicc=TRUE)| Dependent variable | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 1.12 | 0.59 – 2.11 | 0.727 |
| HVivos1 | 2.76 | 2.10 – 3.63 | <0.001 |
| HVivos2 | 5.17 | 3.87 – 6.90 | <0.001 |
| HVivos3 | 7.34 | 5.44 – 9.91 | <0.001 |
| HVivos4 | 9.14 | 6.65 – 12.57 | <0.001 |
| ns(V012, 4)1 | 0.20 | 0.11 – 0.35 | <0.001 |
| ns(V012, 4)2 | 0.30 | 0.19 – 0.48 | <0.001 |
| ns(V012, 4)3 | 0.07 | 0.02 – 0.25 | <0.001 |
| ns(V012, 4)4 | 0.05 | 0.03 – 0.07 | <0.001 |
| Nunca casada | Reference | ||
| Casado | 0.60 | 0.43 – 0.84 | 0.003 |
| Viviendo juntos | 0.61 | 0.44 – 0.84 | 0.002 |
| Viuda | 0.77 | 0.15 – 4.10 | 0.763 |
| Divorciada | 0.79 | 0.09 – 6.65 | 0.828 |
| No viven juntos | 0.75 | 0.49 – 1.14 | 0.183 |
| Sin educación | Reference | ||
| Primario | 1.07 | 0.76 – 1.50 | 0.703 |
| Secundario | 1.52 | 1.07 – 2.15 | 0.019 |
| Mayor | 1.60 | 1.11 – 2.32 | 0.012 |
| El más pobre | Reference | ||
| Pobrer | 1.46 | 1.26 – 1.69 | <0.001 |
| Medio | 1.75 | 1.47 – 2.10 | <0.001 |
| Rico | 2.10 | 1.71 – 2.58 | <0.001 |
| Más rico | 3.01 | 2.33 – 3.88 | <0.001 |
| Observations | 19826 | ||
| R2 / R2 adjusted | 0.086 / -5.096 | ||
Esta tabla muestra los intervalos de confianza así como los coeficientes medidos en riesgo relativo, que es igual a la exponencial del coeficiente. El uso de etiquetas y el mostrar las categorías de referencia hacen mucho más claro lo que tenemos. Hemos añadido también criterios de selección que nos pueden permitir elegir entre modelos, pero desafortunadamente no los muestra la tabla. Las funciones splines y polinómica no son interpretables por sus coeficientes, aunque la alta significatividad estadística muestra que son importantes. Necesitamos paquetes que muestren cómo cambia el valor predicho al cambiar el resto de variables. Del resto de variables ha cambiado el patrón del estado civil: ahora las nunca casadas aparecen de nuevo como las que más usan. Esto puede deberse al tratamiento polinomial del número de hijos. Habría que explorar también las interacciones entre el número de hijos y el estado civil.
Puesto que la tabla no muestra los criterios de selección de modelos lo haremos manualmente con AIC: escogemos el modelo con menor AIC, que en nuestro caso vemos que es el modelo 2. Esta misma estrategia la podemos emplear para explorar los modelos con interacciones. El modelo con menor AIC tiende a coincidir con el que mejor predice fuera de la muestra.
AIC(mod1, mod2)## eff.p AIC deltabar
## [1,] 79.76155 24389.92 3.798169
## [2,] 81.80227 23958.19 4.090113El efecto de la edad lo podemos representar con algún paquete auxiliar como visreg
library(visreg)
mod2 %>% visreg("V012",scale="response")
que me muestra que los jóvenes utilizan más métodos anticonceptivos modernos para unos mismos valores del resto de explicativas. Esto puede parece poco intuitivo, pero debemos tener en cuenta que hay otra variable, el número de hijos vivos, que aumenta con la edad y está también en el modelo. La caída en edades más avanzadas puede tener que ver con mujeres que se saben no fértiles. Podríamos haber filtrado las mujeres postmenopáusicas en los datos, pero no lo hemos hecho. Para poder interpretar bien el patrón por edad es necesario representarlo conjuntamente con el número de hijos vivos:
mod2 %>% visreg2d("HVivos","V012",scale="response")
Vemos como las que utilizan poco son las que han tenido pocos hijos, y sobre todo ninguno, a edades tardías. También podemos comprender ahora mejor lo que representa la mayor probabilidad de No casadas: la manera de ajustar al alza el uso para las mujeres jóvenes sin hijos que están en esa situación y que son distintas de las que se han casado, quieren tener su primer hijo, y por tanto no usan.
Ejercicios
- Explorar modelos con interacciones entre variables.
- Recodificar la variable de estado civil en 2 o 3 categorías y reestimar.
- Contrastar el modelo estimado y los SE con los obtenidos utilizando la función
glm, con pesos y sin pesos, y las matrices de varianzas-covarianzas habitual yvcovCL.
Conclusiones
Este curso ha tratado de mostrar las características especiales de trabajar con datos de encuesta en R, cubriendo todas las fases del análisis empírico y solucionando los problemas según se nos presentan en ejemplos reales.
Las principales conclusiones que obtenemos son:
- Es importante que los ficheros de datos estén bien definidos. Buena parte de los problemas de trabajar con los datos de la Encuesta de Fecundidad, por ejemplo, se debían a no tener declarado el tipo de codificación y tener errores de tabulación, con la encuesta DHS no hemos encontrado ese tipo de problemas.
- Es posible trabajar con las variables en formato
labelledcon paquetes como sjmisc, memisc o sjPlot, aunque en ocasiones puede ser mejor transformar todas las variables en factores ordinarios para poder trabajar con funciones que no entienden las etiquetas. - Las encuestas complejas difieren en gran medida de los supuestos de muestreo aleatorio simple que suele asumirse en los cursos de microeconometría. Esto tiene consecuencias importantes para nuestro análisis que de no ser consideradas invalidan nuestro trabajo.
- Es importante utilizar pesos en el análisis. De otro modo los estimadores puntuales no serán válidos, y tampoco la inferencia debido a que sobreestimamos la precisión de nuestros estimadores.
- Es importante definir los pesos como pesos estandarizados si no vamos a utilizar paquetes no especializados en encuestas como survey. Algunos paquetes interpretarían los pesos como frecuencias llevando a inferencia errónea.
- A la relación entre la varianza de los estimadores al tener en cuenta la estructura compleja y los basados en muestreo aleatorio simple se le conoce efecto del diseño. En nuestro caso estaba en torno a 4.
- Si utilizamos pesos, no es tan necesario tener en cuenta el efecto del diseño, aunque es posible con
survey. Una alternativa son las covarianzas robustas por conglomerados devcovCL. - Las funciones del tidyverse y las tibbles son muy útiles tanto en la fase inicial de transformación de los datos como en la de análisis de resultados.
- Para gráficos y tablas de publicación ggplot es una opción excelente, así como las tablas y gráficos de sjPlot. Todos ellos dentro de documentos Rmd permite realizar un análisis completo y replicable.
Referencias
Cabrera León, A. (2016) “Nuevas técnicas de investigación por muestreo aplicadas a Encuestas de Salud: Calibración de estimadores”, Tesis Doctoral, U. de Granada.
CDC (2016) “Youth Risk Behavior Surveillance System (YRBSS) Software for Analysis of YRBS Data”
Damico, A. J. (2020) “Analyze survey data for free”, asdfree.com
Gutiérrez, A. (2017) “Muestreo y análisis de estudios educacionales con R”
Jiménez, F. (2014) “Tutorial de muestreo con R”. Más un ejemplo de aplicación
Lumley, T. (2010) “Complex Surveys: A Guide to Analysis Using R”, Wiley. Ver también página web del paquete survey y materiales de un curso
Lumley, T., y Scott, A. (2017). “Fitting Regression Models to Survey Data”. Statistical Science, 32 (2), 265-278. 10.1214/16-STS605
Mayor Gallego (2015) “Análisis de muestras obtenidas mediante diseños complejos mediante R y el package ‘SURVEY’”, Univ. de Sevilla.
UCLA: Statistical Consulting Group (2020) “Survey Data Analysis with R” (accedido el 20-10-2020)
UN Statistics Division (2005) “Household Sample Surveys in Developing and. Transition Countries”, United Nations Publication, Sales No. E.05.XVII.6