Sesión 1 - Lectura de encuestas INE y análisis microeconométrico
EncuestaR - Análisis de Encuestas con R. Curso del Doctorado en Economía 2020-21.
José Antonio Ortega
08 de noviembre de 2020
El minicurso
Esta es la 4ª edición de cursos de R dentro del programa de Doctorado en Economía de la Universidad de Salamanca. En esta ocasión lo dedicaremos a temas de gran utilidad potencial para los estudiantes de doctorado, como es el análisis de encuestas. El curso se basa en el desarrollo de ejemplos completos replicables que muestran las posibilidades de R para el análisis de encuestas y, al cubrir todas las fases del análisis desde la consecución de los datos hasta la redacción de informes, enseña a los estudiantes a abordar los problemas que surgen en todo análisis empírico.
En particular cubriremos dos tipos de temas: La primera sesión ilustrará técnicas básicas de manipulación de encuestas desde el formato típico en el que el INE proporciona las encuestas recientes, presentando técnicas microeconométricas útiles para el análisis y para la generación de informes. La segunda sesión mostrará el impacto que tiene el diseño de la encuesta sobre la variabilidad de los estimadores con la ayuda del paquete survey. Las encuestas internacionales generalmente se basan en un diseño complejo polietápico por conglomerados, con selección de unidades primarias de muestreo. Esta estructura debe tenerse en cuenta en el análisis, tanto por la probabilidad de selección diferente para cada individuo como por las posibilidades de variabilidad diferente en cada unidad. Se proporcionarán ejemplos completos que trabajaremos juntos, bibliografía y referencias.
Aunque el curso está pensado para estudiantes del doctorado en economía puede ser de utilidad para investigadores que trabajen con encuestas, profesores y estudiantes de doctorado, tanto en el mundo de la economía como en el de la empresa, ciencias sociales y del comportamiento, o ciencias de la salud.
Si no tienes una buena base de R con el tidyverse deberías repasar referencias como r for data science, traducido al español, o Manipulación de datos e investigación reproducible en R de Corcoran
Datos de encuestas del INE: La Encuesta de Fecundidad 2018
Para la investigación empírica en nuestro país probablemente el principal proveedor de datos sea el INE. Aunque proporciona y lleva proporcionando datos mucho tiempo, los formatos en que aparecen muchas veces no son todo lo amigables que pudiéramos desear.
Muchas de las encuestas más antiguas están codificadas en el formato de ancho fijo (FWF, fixed-width format). Este es particularmente complicado de leer porque requiere definir previamente la correspondencia entre variables y columnas a partir del libro de códigos. Una vez hecho esto, eso sí, esto no es mayor problema utilizando la función read_fwf del paquete readr. Hay un paquete de R que puede ayudar en algunos casos al tener predefinidos los libros de códigos (ejemplo: para las encuestas de presupuestos familiares) MicroDatosES. Yo no lo voy a cubrir porque no está en desarrollo activo y porque no nos ayuda con las encuestas más recientes.
En las encuestas más recientes, el INE proporciona un único fichero comprimido en formato zip que, en principio, contiene los ficheros en formatos múltiples para ser leídos desde R, SAS, SPSS o STATA, y en formato CSV. En la práctica, las cosas son un poco más complicadas: si bien para SAS, SPSS y STATA los datos están en formato interno de dichos programas y se pueden cargar directamente, en el de R no es el caso: se proporciona un código R para ejecutarlo que permite leer desde R los ficheros a partir del formato csv. Sin embargo, no recomiendo proceder de este modo puesto que las variables categóricas aparecen codificadas con códigos numéricos cuyo significado sólo aparece en el libro de códigos. Es mejor trabajar a partir del formato sav de SPSS. Pero esto sólo resuelve parcialmente el problema, como veremos.
Como ejemplo de análisis vamos a adoptar la Encuesta de Fecundidad 2018. Podemos acceder a las características técnicas a partir de aquí. Aquí están enlazados el informe metodológico, el cuestionario de hombres y el de mujeres. Esta es la página con los microdatos.
Se trata realmente de dos encuestas: una a hombres (N=2619) y otra a mujeres (N=14556). El diseño es bietápico: en la primera etapa se seleccionan secciones censales (510 para hombres, 1825 para mujeres), y dentro de cada sección censal, se selecciona aleatoriamente una muestra de la población de 18 a 55 años de edad (10 hombres o 14 mujeres por sección).
Puesto que nuestro interés es mostrar un ejemplo, vamos a trabajar con los datos de hombres que, al ser una muestra más pequeña, nos permite hacer los análisis más rápido y ocupando menos memoria, aunque esto no debería ser un problema en la mayoría de los casos.
Asimismo, para cada sexo tenemos tres objetos distintos: el fichero individual, el fichero de hogar y el fichero de hijos. Nosotros nos vamos a concentrar en trabajar con el fichero individual, pero se pueden combinar los distintos ficheros con left_join a partir del identificador de cuestionario IDEN. En todos los casos la clave es trabajar con el diseño de registro que nos especifica qué mide cada variable, qué valores puede tomar, y qué significan.
Aunque podríamos bajar “a mano” los ficheros, vamos a hacerlo desde R para mantener el análisis reproducible. Utilizamos la función walk del paquete purrr. Este paquete dispone de herramientas de análisis funcional que reemplazan lo que en otros lenguajes se haría con bucles. Las funciones que producen resultados en nuestro espacio de trabajo de R son variantes de map. Las que no producen como resultado un objeto de R son variantes de walk. Es la que utilizamos: descargamos varios ficheros con download.file cuyas url hemos guardado en el objeto fich. Hemos descomprimido los dos ficheros zip y guardado en la carpeta data, que debemos tener creada de antemano, únicamente los ficheros .sav
fich=c("ftp://www.ine.es/temas/fecun/Hombres_Tabla_Hijos.zip",
"ftp://www.ine.es/temas/fecun/Hombres_Individual.zip",
"ftp://www.ine.es/temas/fecun/dr_EFhijos_2018.xlsx",
"ftp://www.ine.es/temas/fecun/dr_EFindividual_2018.xlsx"
)
library(purrr)
walk(fich,~download.file(.x,destfile=paste0("data/",basename(.x)), mode="wb"))
temp <- tempfile()
dir("data",pattern="zip") %>% paste0("data/",.) %>% walk(unzip,exdir=temp)
file.copy( paste0(temp,"/SPSS/",dir(paste0(temp,"/SPSS"))),
"data/" )
## Para guardar también el fichero .R y el .txt,
# aunque no los utilizaremos.
f_keep <- c("/md_Hombres_Individual.txt",
"/R/MD_Hombres_Individual.R")
file.copy(paste0(temp,f_keep),"data/")Lectura de datos en formato SPSS
Ya tenemos en nuestra carpeta los ficheros .sav que nos interesan. Existen, básicamente, cuatro alternativas para leer los ficheros desde R. Todos ellos también tienen funciones equivalente para fichero en formato STATA:
| Paquete | Función | Resultado |
|---|---|---|
| foreign | read.spss |
list o data.frame con atributo labels |
| haven | read_sav |
tibble clase haven_labelled |
| sjlabelled | read_spss |
tibble clase labelled |
| memisc | spss.system.file + as.data.set |
data.set convertible a data.frame |
En principio los cuatro preservan la información de las etiquetas, aunque son bastante distintos en la práctica. Estos son algunos de los factores para elegir:
- Tanto memisc como foreign son capaces de leer bien las etiquetas en español, algo que no he conseguido con haven ni con sjlabelled. Ambos tienen un argumento
encoding, pero no he podido conseguir que funcione bien con ninguna especificación. - Los formatos
labelledyhaven_labelledfuncionan como un factor que preserva los niveles originales pero que mantiene las etiquetas que entienden algunos paquetes, como el propio haven y sjmisc. Para pasarlo a un factor que use las etiquetas necesitamos utilizarto_factorde sjmisc. - El formato
data.setdememiscse puede convertir tanto alhaven_labelledcomo a undata.frameque usa las etiquetas como niveles. Posiblemente esto es lo más cómodo. - foreign permite elegir si se importa a una lista o a un data.frame. En este último caso se puede elegir si se utilizan las etiquetas o los niveles originales. Además en el caso de variables que contengan códigos repetidos, algo que ocurre en este fichero, foreign lo reasigna automáticamente dando un mensaje de advertencia. memisc necesita que se reasignen manualmente, aunque proporciona las herramientas para hacerlo.
- Otra posibilidad es la de leer las etiquetas y las preguntas desde el diseño de registro excel y añadírselos a los datos leídos en uno de los formatos anteriores, incluido el
csv.
Vamos a leerlo desde los tres primeros paquetes. En un Rmd externo (codebook.Rmd) se presenta cómo hacerlo con memisc. También lo leo desde foreign de dos maneras: con valores y con etiquetas. De este modo, si se desea cambiar a la otra opción es posible variable a variable reasignando la correspondiente del otro objeto.
library(tidyverse)
# EF_txt=readr::read_table("data/Hombres_Individual.csv")
# Con foreign, como etiquetas y como valores
EF=foreign::read.spss("data/Hombres_Individual.sav",
to.data.frame=TRUE,
encoding="iso8859-1")
EF_val=foreign::read.spss("data/Hombres_Individual.sav",
to.data.frame=TRUE,
use.value.labels = FALSE,
encoding="iso8859-1")
# Con haven
EF_haven=haven::read_sav("data/Hombres_Individual.sav",
encoding="iso8859-1")
# Con sjlabelled
EF_sj=sjlabelled::read_spss("data/Hombres_Individual.sav",enc="iso8859-1")Veamos ahora las diferencias de formato haciendo una tabulación elemental con la función frq del paquete sjmisc, que permite hacer tabulaciones cruzadas para cualquiera de los formatos.
library(sjmisc)
## Con frq. La opción out="v", permite que salga en html bonito
sjmisc::frq(EF$VIVIR15P1,out="v")| val | label | frq | raw.prc | valid.prc | cum.prc | |
|---|---|---|---|---|---|---|
| Sí |
|
1582 | 60.4 | 96.64 | 96.64 | |
| No, había fallecido |
|
17 | 0.65 | 1.04 | 97.68 | |
| No, por otros motivos |
|
38 | 1.45 | 2.32 | 100 | |
| NA | NA | 982 | 37.5 | NA | NA | |
| total N=2619 · valid N=1637 · x̄=1.06 · σ=0.32 | ||||||
sjmisc::frq(EF_val$VIVIR15P1)##
## x <character>
## # total N=2619 valid N=2619 mean=1.66 sd=0.57
##
## Value | N | Raw % | Valid % | Cum. %
## ---------------------------------------
## | 982 | 37.50 | 37.50 | 37.50
## 1 | 1582 | 60.40 | 60.40 | 97.90
## 2 | 17 | 0.65 | 0.65 | 98.55
## 3 | 38 | 1.45 | 1.45 | 100.00
## <NA> | 0 | 0.00 | <NA> | <NA># Aunque en pantalla sale mal, con out="v" EF_haven funciona bien.
sjmisc::frq(EF_haven$VIVIR15P1,out="v")| val | label | frq | raw.prc | valid.prc | cum.prc | |
|---|---|---|---|---|---|---|
| 1 | 982 | 37.5 | 37.5 | 37.5 | ||
| 2 | Sí | 1582 | 60.4 | 60.4 | 97.9 | |
| 3 | No, había fallecido | 17 | 0.65 | 0.65 | 98.55 | |
| 4 | No, por otros motivos | 38 | 1.45 | 1.45 | 100 | |
| NA | NA | 0 | 0 | NA | NA | |
| total N=2619 · valid N=2619 · x̄=1.66 · σ=0.57 | ||||||
# Lo mismo ocurre con EF_sj. Vemos que como texto mal.
sjmisc::frq(EF_sj$VIVIR15P1)##
## x <character>
## # total N=2619 valid N=2619 mean=1.66 sd=0.57
##
## Value | Label | N | Raw % | Valid % | Cum. %
## ----------------------------------------------------------------
## 1 | | 982 | 37.50 | 37.50 | 37.50
## 2 | S<ed> | 1582 | 60.40 | 60.40 | 97.90
## 3 | No, hab<ed>a fallecido | 17 | 0.65 | 0.65 | 98.55
## 4 | No, por otros motivos | 38 | 1.45 | 1.45 | 100.00
## <NA> | <NA> | 0 | 0.00 | <NA> | <NA>Vamos ahora a leer los metadatos desde el fichero excel que contiene el diseño de registro. No sólo leemos los nombres de las variables en la primera hoja, sino que a cada una le asignamos los valores que toma la variable según se especifica en las hojas sucesivas. El objeto general formado, meta, lo podemos utilizar como referencia o para reasignar etiquetas de variable o nivel a las variables. En nuestro caso hemos asignado etiquetas de variable al objeto EF que ya teníamos disponible.
Comprobamos de nuevo la utilidad de funciones del tidyverse como pueden ser mutate para transformar variables, rename, para cambiar el nombre, fill, para rellenar columnas copiando de celdas superiores, y las funciones map de purrr, como map_df, para organizar los resultados de aplicar una función a un conjunto de datos.
## Leyendo los metadatos desde excel.
# - Usamos una manipulación "sencilla" del excel en tidyverse para
# cargar los metadatos basada también en `purrr`.
# - Nos sirve de diccionario
## PRECAUCIÓN: Da error de lectura si tenemos abierto el fichero en excel al mismo tiempo.
library(readxl)
metaind="data/dr_EFindividual_2018.xlsx"
metasheets=readxl::excel_sheets(metaind)
meta0=read_excel(metaind,metasheets[1],skip=1) %>%
dplyr::rename(Diccionario=`Diccionario de la variable`,
Bloque=`Bloque del cuestionario`) %>%
tidyr::fill(Bloque)
metaall=map_df(metasheets[-1],read_excel,path=metaind,
.id="tabla",col_names=c("cod","descr"),
col_type="text",range=cell_cols(1:2)) %>%
mutate(Diccionario=if_else((lead(cod)=="Código"),cod,
NA_character_)) %>%
tidyr::fill(Diccionario) %>% filter(!is.na(descr),
cod!="Código")
meta= left_join(meta0,metaall,by="Diccionario")
# Poniendo etiquetas de variable a partir de meta0
# Hay que quitar la última "variable", total.
# Esto asigna las preguntas a cada una de las variables
labelled::var_label(EF) <- meta0$Descripción %>% setNames(meta0$Variable) %>% as.list %>% .[-486]
# Esto pondría los niveles correctos en cada variable,
# para los objetos que no los tienen, pero
# también falla por el problema de las tres variables con
# etiquetas no únicas.
# labelled::val_labels(EF_sj) <- meta %>% filter(!is.na(Diccionario)) %>% mutate(descr=setNames(descr,cod)) %>% pull(descr) %>% split(meta %>% filter(!is.na(Diccionario)) %>% pull(Variable))Vamos a ver un ejemplo de table con frq después de incorporadas las etiquetas. Así salen en pantalla:
EF %>% frq(VIVIR15P1)##
## Vivía su progenitor 1 con usted cuando tenía 15 años (VIVIR15P1) <categorical>
## # total N=2619 valid N=1637 mean=1.06 sd=0.32
##
## Value | N | Raw % | Valid % | Cum. %
## -------------------------------------------------------
## Sí | 1582 | 60.40 | 96.64 | 96.64
## No, había fallecido | 17 | 0.65 | 1.04 | 97.68
## No, por otros motivos | 38 | 1.45 | 2.32 | 100.00
## <NA> | 982 | 37.50 | <NA> | <NA>Y otro del uso de pesos en la tabulación, que permite replicar los resultados que da el INE, con la opción out="v".
EF %>% frq(USOANTICONCEP,weights=FACTOR, out="v")| val | label | frq | raw.prc | valid.prc | cum.prc | |
|---|---|---|---|---|---|---|
| Si |
|
6275880 | 52.17 | 52.17 | 52.17 | |
| No |
|
5754191 | 47.83 | 47.83 | 100 | |
| NA | NA | 0 | 0 | NA | NA | |
| total N=12030071 · valid N=12030071 · x̄=1.48 · σ=0.50 | ||||||
Tabulaciones cruzadas y otras utilidades
flat_table
Una de las principales técnicas de análisis con encuestas son las tabulaciones cruzadas. En la segunda sesión haremos un tratamiento más detallado sobre su estimación en encuestas complejas. Una función sencilla pero potente es sjmisc::flat_table. Si queramos únicamente el recuento sencillamente tenemos que especificar las variables que queremos tabular separadas por comas:
EF %>% flat_table(NHIJOBIO,USOANTICONCEP)## USOANTICONCEP Si No
## NHIJOBIO
## 0 710 692
## 1 218 240
## 2 362 270
## 3 49 48
## 4 12 15
## 5 2 1El resultado es de clase ftable, pero se puede convertir en tibble para su manipulación posterior. Uno de los motivos que la hacen útil con encuestas es que permite especificar un vector de pesos. Como vemos aquí
EF %>% flat_table(NHIJOBIO,USOANTICONCEP) %>%
as_tibbleSi especificamos como vector de pesos el factor de elevación, FACTOR, obtenemos la estimación del efectivo total a partir de la encuesta que coincide con las cifras que proporciona el INE.
EF %>% flat_table(NHIJOBIO,USOANTICONCEP,weights=FACTOR)## USOANTICONCEP Si No
## NHIJOBIO
## 0 3346631 3160395
## 1 975630 1073269
## 2 1650733 1215763
## 3 233779 215344
## 4 57416 84373
## 5 11692 5046Podemos ver la diferencia de utilizar pesos o no en los porcentajes calculados:
EF %>% flat_table(NHIJOBIO,USOANTICONCEP,margin="col")## USOANTICONCEP Si No
## NHIJOBIO
## 0 52.48 54.66
## 1 16.11 18.96
## 2 26.76 21.33
## 3 3.62 3.79
## 4 0.89 1.18
## 5 0.15 0.08EF %>% flat_table(NHIJOBIO,USOANTICONCEP,weights=FACTOR,margin="col")## USOANTICONCEP Si No
## NHIJOBIO
## 0 53.33 54.92
## 1 15.55 18.65
## 2 26.30 21.13
## 3 3.73 3.74
## 4 0.91 1.47
## 5 0.19 0.09No son excesivamente diferentes, pero si queremos que las cifras se correspondan con las publicadas por el INE y que sean representativas de la población española, debemos usar el factor de elevación como peso. Podemos ver, con un gráfico ggplot, que los factores de elevación de los distintos individuos pueden ser muy distintos:
EF %>% ggplot(aes(y=FACTOR),alpha=0.6) +geom_boxplot() +
theme_minimal()
Podemos ver como hay mucha diferencia entre factores de elevación, que pueden llegar a un rango de 10 a 1. Con la función sjmisc::descr podemos obtener estadísticos descriptivos de la variable que nos lo confirman.
EF %>% descr(FACTOR)En la segunda sesión hablaremos más de cómo se calculan estos pesos y qué significa lo de que estén calibrados por sexo, edad y provincia. Eso sí: nosotros no tenemos información en los datos sobre la provincia, ni siquiera sobre la CCAA.
Tabulaciones cruzadas con frq + group_by
Cuando utilizamos group_by en nuestro objeto de datos lo dividimos en subtablas asociadas a cada nivel de la combinación de variables generadas. Esto lo podemos hacer, por ejemplo, para ver la distribución de la población por número de hijos para cada ámbito residencial con frq
EF %>% filter(EDAD>45) %>% group_by(DEGURBA) %>%
frq(NHIJOBIO,weights = FACTOR,out="v",encoding="latin1")| val | label | frq | raw.prc | valid.prc | cum.prc | |
|---|---|---|---|---|---|---|
| 0 |
|
621081 | 33.02 | 33.02 | 33.02 | |
| 1 |
|
404943 | 21.53 | 21.53 | 54.55 | |
| 2 |
|
706070 | 37.54 | 37.54 | 92.08 | |
| 3 |
|
85203 | 4.53 | 4.53 | 96.61 | |
| 4 |
|
51773 | 2.75 | 2.75 | 99.37 | |
| 5 |
|
11925 | 0.63 | 0.63 | 100 | |
| NA | NA | 0 | 0 | NA | NA | |
| total N=1880995 · valid N=1880995 · x̄=1.24 · σ=1.08 | ||||||
| val | label | frq | raw.prc | valid.prc | cum.prc | |
|---|---|---|---|---|---|---|
| 0 |
|
383526 | 29.86 | 29.86 | 29.86 | |
| 1 |
|
302761 | 23.57 | 23.57 | 53.43 | |
| 2 |
|
509009 | 39.63 | 39.63 | 93.06 | |
| 3 |
|
60110 | 4.68 | 4.68 | 97.74 | |
| 4 |
|
24185 | 1.88 | 1.88 | 99.63 | |
| 5 |
|
4813 | 0.37 | 0.37 | 100 | |
| NA | NA | 0 | 0 | NA | NA | |
| total N=1284404 · valid N=1284404 · x̄=1.26 · σ=1.02 | ||||||
| val | label | frq | raw.prc | valid.prc | cum.prc | |
|---|---|---|---|---|---|---|
| 0 |
|
204134 | 39.99 | 39.99 | 39.99 | |
| 1 |
|
96528 | 18.91 | 18.91 | 58.91 | |
| 2 |
|
176933 | 34.67 | 34.67 | 93.57 | |
| 3 |
|
23541 | 4.61 | 4.61 | 98.18 | |
| 4 |
|
9270 | 1.82 | 1.82 | 100 | |
| NA | NA | 0 | 0 | NA | NA | |
| total N=510406 · valid N=510406 · x̄=1.09 · σ=1.04 | ||||||
donde se aprecia que la probabilidad de no tener hijos tras los 45 años es mínima en la zona intermedia. Nos puede interesar más ver el número medio de hijos. Este podríamos calcularlo con summarize utilizando la función `memisc::Weighted.Mean:
EF %>% filter(EDAD>45) %>% group_by(DEGURBA) %>%
summarize(Descendencia=memisc::Weighted.Mean(NHIJOBIO,w=FACTOR))O, si queremos estadísticos descriptivos de la variable, podemos emplear sjmisc::descr (esto no permite pesos)
EF %>% filter(EDAD>45) %>% group_by(DEGURBA) %>%
descr(NHIJOBIO)##
## ## Basic descriptive statistics
##
##
## Grouped by: Urbano
##
## var type label n NA.prc mean sd se md
## NHIJOBIO numeric Número de hijos biológicos 420 0 1.25 1.08 0.05 1
## trimmed range iqr skew
## 1.16 5 (0-5) 2 0.51
##
##
## Grouped by: Intermedio
##
## var type label n NA.prc mean sd se md
## NHIJOBIO numeric Número de hijos biológicos 272 0 1.24 1.01 0.06 1
## trimmed range iqr skew
## 1.19 5 (0-5) 2 0.34
##
##
## Grouped by: Rural
##
## var type label n NA.prc mean sd se md
## NHIJOBIO numeric Número de hijos biológicos 121 0 1.09 1.05 0.1 1
## trimmed range iqr skew
## 1.01 4 (0-4) 2 0.43El paquete weights proporciona funciones para el cálculo de estadísticos descriptivos, tablas y contrastes de hipótesis utilizando pesos, pero no lo vamos a utilizar aquí.
Transformaciones de variables
Antes de poder pasar a estimar modelos, es muy posible que tengamos que transformar las variables originales. Los casos más típicos de transformaciones requiridas son:
- Calcular nuevas variables numéricas: Esto en general lo haremos con
mutatedel tidyverse. - Recodificar variables categóricas, combinando variables, reduciendo niveles, … Para esto son muy útiles los paquetes
forcats,sjmisc(funciónrec), y el paquetestringrpara manipulaciones de cadenas de texto. - Para factores con etiquetas muy largas puede ser conveniente asignar nuevos niveles con nombres más cortos. En este caso en concreto los niveles numéricos originales no son útiles, por eso nos hemos quedado con las etiquetas como niveles, pero si las etiquetas son demasiado largas podemos reasignarlas como etiquetas y crear nuevos niveles más cortos. Para ello utilizaríamos el paquete
sjlabelled. - Manejar variables de tipo fecha. Para ello el paquete más útil es
lubridate.
Vamos a ver algunos ejemplos típicos de recodificación, algunos los haremos, otros los dejamos indicados como necesarios si se quieren analizar esas variables
- Hay variables continuas, por ejemplos las de edades a las que se hacen cosas, que tienen una categoría especial que se suele codificar como
99. En general habrá que recodificarlo como NA. Ejemplo:
EF %>% frq(EDADPRS)##
## Edad a la que tuvo su primera relación sexual (EDADPRS) <numeric>
## # total N=2619 valid N=2475 mean=18.97 sd=10.40
##
## Value | N | Raw % | Valid % | Cum. %
## --------------------------------------
## 11 | 7 | 0.27 | 0.28 | 0.28
## 12 | 14 | 0.53 | 0.57 | 0.85
## 13 | 41 | 1.57 | 1.66 | 2.51
## 14 | 108 | 4.12 | 4.36 | 6.87
## 15 | 228 | 8.71 | 9.21 | 16.08
## 16 | 412 | 15.73 | 16.65 | 32.73
## 17 | 483 | 18.44 | 19.52 | 52.24
## 18 | 515 | 19.66 | 20.81 | 73.05
## 19 | 191 | 7.29 | 7.72 | 80.77
## 20 | 171 | 6.53 | 6.91 | 87.68
## 21 | 82 | 3.13 | 3.31 | 90.99
## 22 | 50 | 1.91 | 2.02 | 93.01
## 23 | 26 | 0.99 | 1.05 | 94.06
## 24 | 24 | 0.92 | 0.97 | 95.03
## 25 | 27 | 1.03 | 1.09 | 96.12
## 26 | 14 | 0.53 | 0.57 | 96.69
## 27 | 8 | 0.31 | 0.32 | 97.01
## 28 | 10 | 0.38 | 0.40 | 97.41
## 29 | 6 | 0.23 | 0.24 | 97.66
## 30 | 11 | 0.42 | 0.44 | 98.10
## 32 | 1 | 0.04 | 0.04 | 98.14
## 33 | 2 | 0.08 | 0.08 | 98.22
## 34 | 1 | 0.04 | 0.04 | 98.26
## 35 | 1 | 0.04 | 0.04 | 98.30
## 37 | 1 | 0.04 | 0.04 | 98.34
## 39 | 2 | 0.08 | 0.08 | 98.42
## 40 | 1 | 0.04 | 0.04 | 98.46
## 99 | 38 | 1.45 | 1.54 | 100.00
## <NA> | 144 | 5.50 | <NA> | <NA>En este caso tendríamos que recodificar el 99 como NA. Esto marca, en primer lugar, la importancia de mirar las variables con frq antes de hacer el análisis. Una vez identificado el problema es fácil con la función na_if. Como el problema no es sólo de una variable, para todas las variables que contengan la cadena EDAD, vamos a recodificar los 99 como NA. Para ver qué variables son estas podemos verlo con find_var
EF %>% find_var("EDAD")Sirve para identificar las variables. En nuestro caso constatamos que no todas las variables que contienen EDAD son edades: solo las que empiezan por edad y son numéricas. Hacemos por lo tanto el siguiente cambio
EF = EF %>% mutate(across(starts_with("EDAD")&(where(is.numeric)),na_if,y=99,.names="{col}_R"))La función across sirve para controlar dónde se hace la mutación. Debemos de ser conscientes que al hacer estos cambios hemos perdido las variables originales sino ponemos otros nombres. El 99 tiene un significado distinto del NA: en general se refiere a que aún no se ha producido la transición de la que se habla. Si queremos mantener ambas hay que poner un sufijo, por ejemplo _R. No era necesario hacerlo para todas las variables, pero nos es más complicado encontrar cuáles. Es buena práctica mantener las dos variables puesto que el 99 puede tener sentido en algunas de ellas, por ejemplo EDADP1 y EDADP2, las edades de los padres. De este modo tenemos las dos variables preparadas y en el momento del análisis podemos ver cuál nos interesa, la original o la recodificada con el sufijo _R.
- La variable
EDAD, en las tabulaciones del INE, aparece codificada por tramos. Existen varias funciones que nos pueden ayudar a recodificarla por tramos:sjmisc::recnos permite mantenerla como labelled,cutes una función de RBase que también funciona, y también es útildplyr::case_when. Vamos a hacerlo con cut. La función permite especificar las etiquetas con el argumentolabels, pero a mí personalmente me parecen claras las etiquetas por defecto en forma de intervalo.
EF = EF %>%
mutate(
EDADG=cut(EDAD,breaks=c(18,30,35,40,45,50,55),right=FALSE,
include.lowest=TRUE))
EF %>% frq(EDADG)##
## EDADG <categorical>
## # total N=2619 valid N=2619 mean=3.46 sd=1.81
##
## Value | N | Raw % | Valid % | Cum. %
## ----------------------------------------
## [18,30) | 612 | 23.37 | 23.37 | 23.37
## [30,35) | 271 | 10.35 | 10.35 | 33.72
## [35,40) | 403 | 15.39 | 15.39 | 49.10
## [40,45) | 451 | 17.22 | 17.22 | 66.32
## [45,50) | 389 | 14.85 | 14.85 | 81.18
## [50,55] | 493 | 18.82 | 18.82 | 100.00
## <NA> | 0 | 0.00 | <NA> | <NA>Podemos comprobar ahora si replicamos exactamente la tabla de uso de anticonceptivos:
EF %>% flat_table(EDADG,USOANTICONCEP,weights=FACTOR)## USOANTICONCEP Si No
## EDADG
## [18,30) 1742088 1128589
## [30,35) 750164 637473
## [35,40) 950100 798384
## [40,45) 1099602 906901
## [45,50) 985523 913497
## [50,55] 748405 1369345Podemos comprobar que es así.
- Pasar a formato fecha las que vienen determinadas por año y mes (con prefijos ANO y MES. Podéis buscarlas con
find_vars). Aquí también tenemos varias opciones: podemos pasarlas a fecha con lubridate utilizando, por ejemplo, el primer día del mes. También podemos codificarla en un formato de año-mes: por ejemployearmondel paquete zoo oyearmonthdel paquete tsibble. Esto es útil si queremos hacer después, por ejemplo, análisis de supervivencia. Voy a hacer tres ejemplos: una fecha de nacimiento (sufijo N) que imputa el 14 del mes a cada caso, una fecha de emancipación (sufijo EMANCIPA) , y una diferencia entre las edades del ego y de la pareja actual en meses (sufijos N y NPAR). Este último lo utilizaremos en un análisis de regresión posterior.
library(zoo) # Para que entienda los yearmon
EF = EF %>%
mutate(
FECHAN=lubridate::ymd(paste(ANON,as.numeric(MESN),"1",sep="-")),
YMEMANCIPA=as.yearmon(paste(ANOEMANCIPA,as.numeric(MESEMANCIPA),sep="-")),
DIFEDAD=(ANONPAR-ANON)+(as.numeric(MESNPAR)-as.numeric(MESN))/12)
## Vemos que se pueden calcular también "medias" de fechas
EF$FECHAN %>% summary()## Min. 1st Qu. Median Mean 3rd Qu. Max.
## "1962-03-01" "1970-06-01" "1978-02-01" "1979-03-14" "1987-03-16" "2000-01-01"## Las fechas de emancipación quedan en formato yearmon
EF$YMEMANCIPA %>% head()## [1] NA NA NA "may. 1987" NA NA## Las diferencias de edades nos muestran el hombre mayor
## que su pareja. Puede haber parejas del mismo sexo,
## están señaladas con SEXOPAR.
EF$DIFEDAD %>% summary()## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -60.0000 -0.3125 1.4167 1.7587 3.6667 27.9167 1117EF %>% ggplot(aes(x=FECHAN,y=NHIJOBIO),alpha=0.5)+geom_hex() + geom_smooth() + scale_fill_steps2() + labs(title="Número de hijos en función de fecha de nacimiento")
EF %>% ggplot(aes(x=FECHAN,y=DIFEDAD)) + geom_jitter(aes(size=FACTOR,color=SEXOPAR),alpha=0.2)+ geom_smooth() + scale_fill_steps2() + labs(title="Diferencia de edad con pareja en función de fecha de nacimiento")
Una de las utilidades de hacer gráficos y resúmenes de los datos es la capacidad de poder detectar posibles atípicos que puedan estar ligados a errores en los datos. Esto es el caso de la diferencia de edad entre cónyuges. En particular hay un par de observaciones sospechosas de contener errores que corresponden a hombres cuya pareja es otro hombre hasta 60 años mayor. Vamos a mirar los datos de más de 30 años de diferencia
EF %>% filter(abs(DIFEDAD)>30) %>%
select(FECHAN,ANONPAR,SEXOPAR,EDAD,DIFEDAD,YMEMANCIPA)Los dos casos parecen sospechosos de ser errores en los datos, en particular la fecha de emancipación parece sospechosa: en un caso no la proporciona y en otro coincide con la de nacimiento. Se podrían estudiar más en detalle las observaciones, pero podría optarse por quitarlas, al menos para el análisis de diferencia de edad con la pareja.
Otra variable que podemos redefinir es la de educación, tanto para el ego, ESTUDIOSA, como para la pareja, ESTUDIOSPAR. Los estudios terminados, por ejemplo, tienen 9 niveles
EF %>% frq(ESTUDIOSA)##
## Nivel de estudios alcanzados (ESTUDIOSA) <categorical>
## # total N=2619 valid N=2619 mean=4.86 sd=2.03
##
## Value | N | Raw % | Valid % | Cum. %
## ---------------------------------------------------------------------------------------------
## Menos que primaria | 64 | 2.44 | 2.44 | 2.44
## Educación primaria | 346 | 13.21 | 13.21 | 15.65
## Primera etapa de educación secundaria y similar (con orienta | 344 | 13.13 | 13.13 | 28.79
## Segunda etapa de educación secundaria y similar (con orienta | 464 | 17.72 | 17.72 | 46.51
## Educación postsecundaria no superior | 295 | 11.26 | 11.26 | 57.77
## Enseñanzas de formación profesional, artes plásticas y diseñ | 494 | 18.86 | 18.86 | 76.63
## Grados universitarios de hasta 240 créditos ECTS, diplomados | 266 | 10.16 | 10.16 | 86.79
## Grados universitarios de hasta 240 créditos ECTS, licenciado | 319 | 12.18 | 12.18 | 98.97
## Enseñanzas de doctorado | 27 | 1.03 | 1.03 | 100.00
## <NA> | 0 | 0.00 | <NA> | <NA>Podemos redefinirlo en tres niveles, por ejemplo con sjmisc::rec que me permite trabajar con etiquetas
EF = EF %>%
mutate(
ESTUD=rec(as.numeric(ESTUDIOSA),
rec="1:2=1;3:6=2;7:9=3",
as.num=FALSE,
val.labels = c("Primaria","Secundaria","Universitaria"),
var.label = "Nivel de estudios"),
ESTUDPAR=rec(as.numeric(ESTUDIOSPAR),
rec="1:2=1;3:6=2;7:9=3",
as.num=FALSE,
val.labels = c("Primaria","Secundaria","Universitaria"),
var.label = "Nivel de estudios de la pareja")
)
EF %>% flat_table(ESTUDIOSA,ESTUD)## ESTUD Primaria Secundaria Universitaria
## ESTUDIOSA
## Menos que primaria 64 0 0
## Educación primaria 346 0 0
## Primera etapa de educación secundaria y similar (con orienta 0 344 0
## Segunda etapa de educación secundaria y similar (con orienta 0 464 0
## Educación postsecundaria no superior 0 295 0
## Enseñanzas de formación profesional, artes plásticas y diseñ 0 494 0
## Grados universitarios de hasta 240 créditos ECTS, diplomados 0 0 266
## Grados universitarios de hasta 240 créditos ECTS, licenciado 0 0 319
## Enseñanzas de doctorado 0 0 27EF %>% flat_table(ESTUDPAR,ESTUD)## ESTUD Primaria Secundaria Universitaria
## ESTUDPAR
## Primaria 135 87 16
## Secundaria 138 711 134
## Universitaria 31 279 312Vemos que ha funcionado bien en la primera tabulación cruzada, y apreciamos en la de estudios del ego y de la pareja la existencia de homogamia.
Vector de pesos estandarizados
Una transformación importante antes de meternos con el análisis econométrico es la estandarización del vector de pesos. El que proporciona la EF es realmente un factor de elevación, a cuánta gente representa un individuo. Algunas de las funciones que utilizan pesos, como por ejemplo las glm, asumen cuando se da un vector de pesos que este corresponde a un vector de frecuencias en las que se observan esos datos. Esto lleva a grandes errores en la estimación de los errores estándar puesto que se piensa que hay muchos más datos de los que realmente hay. Para evitarlo es mejor tener calculados y emplear en todo el análisis econométrico unos pesos estandarizados. Es tan sencillo como dividir FACTOR por su valor medio. De este modo los pesos calculados tienen media de 1 y lo único que representan es si algún tipo de individuos están infra o sobrerepresentados en la muestra
Por prevención debemos asegurarnos de que el vector de pesos en nuestro análisis econométricos esté estandarizado para que los errores estándar sean aproximadamente correctos.
EF = EF %>% mutate(PESO=FACTOR/mean(FACTOR))
EF %>% descr(PESO,FACTOR)Análisis microeconométrico de encuestas
Modelos lineales
En principio la estimación de modelos lineales con datos de encuesta la hacemos igual que con cualquier otro tipo de datos, con lm. Es importante no olvidarnos de incluir la variable de pesos, eso sí. Si vamos a centrarnos en una única variable explicada y esta afecta solo a parte de la muestra, puede ser buena idea crearnos otro objeto de encuesta que contenga a la submuestra que nos interesa.
Vamos a hacerlo en nuestro caso para el estudio de la diferencia de edades en parejas de distinto sexo. Para ello tenemos que seleccionar aquellos que tienen pareja y que esta es mujer.
EFedad = EF %>% filter(SEXOPAR=="Mujer")Vamos a seleccionar como variables explicativas la edad del individuo (EDAD), el nivel de estudios del ego y la pareja (las creadas anteriormente ESTUD y ESTUDPAR), el tipo de union con la pareja (TIPOUNION) y los lugares de nacimiento del ego y de la pareja (NACIM, NACIMPAR). Antes de estimar modelos está bien realizar análisis exploratorio detectando posibles problemas en los datos. Una herramienta muy útil es la función GGally::ggpairs. Vamos a verla en la práctica
library(GGally)
EFedad %>% select(DIFEDAD,EDAD,ESTUD,ESTUDPAR,TIPOUNION,NACIM,NACIMPAR) %>%
ggpairs()
Vamos ahora a estimar el primer modelo y ver sus propiedades. summary nos proporciona las estimaciones y la inferencia basada en la fórmulas con homocedasticidad. El paquete sjPlot contiene una función útil para representar los intervalos de confianza.
medad1=lm(
DIFEDAD~EDAD+ESTUD+ESTUDPAR+TIPOUNION+NACIM+NACIMPAR,
data=EFedad,weights=EFedad$PESO)
summary(medad1)##
## Call:
## lm(formula = DIFEDAD ~ EDAD + ESTUD + ESTUDPAR + TIPOUNION +
## NACIM + NACIMPAR, data = EFedad, weights = EFedad$PESO)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -22.6698 -2.0813 -0.1038 1.8930 20.5352
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.49406 0.75053 -5.988 2.67e-09 ***
## EDAD 0.13649 0.01352 10.096 < 2e-16 ***
## ESTUD2 -0.74255 0.32325 -2.297 0.0218 *
## ESTUD3 -0.43119 0.37659 -1.145 0.2524
## ESTUDPAR2 0.63062 0.36619 1.722 0.0853 .
## ESTUDPAR3 0.41148 0.39524 1.041 0.2980
## TIPOUNIONPareja de hecho registrada 0.77018 0.69466 1.109 0.2677
## TIPOUNIONPareja de hecho sin registrar 1.65994 0.26552 6.252 5.31e-10 ***
## NACIMOtro país -0.23127 0.39298 -0.589 0.5563
## NACIMPAROtro país 2.37750 0.35593 6.680 3.39e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.036 on 1469 degrees of freedom
## (338 observations deleted due to missingness)
## Multiple R-squared: 0.1007, Adjusted R-squared: 0.09516
## F-statistic: 18.27 on 9 and 1469 DF, p-value: < 2.2e-16sjPlot::plot_model(medad1)
La mayor parte de los coeficientes es interpretable. No es un modelo causal sino que trata de analizar la asociación entre distintas dimensiones de homogamia.
Un inconveniente de lo presentado anteriormente es que los errores estándar y los intervalos se basan en las fórmulas con homocedasticidad cuando es muy posible que exista heterocedasticidad. Una primera cosa que podemos hacer es estudiar el cumplimiento de las hipótesis del modelo lineal con los gráficos de diagnóstico básicos:
library(ggfortify)
autoplot(medad1)
El gráfico muestra que hay indicios de heterocedasticidad, clara falta de normalidad con mucha probabilidad en ambas colas, y la existencia de algunas observaciones influyentes. También habría que explorar la linealidad de la relación con EDAD.
Inferencia robusta
Debido a que los errores estándar no están bien calculados a no ser que haya homocedasticidad, algo que parece no darse, siempre es buena práctica utilizar estimadores de la matriz de varianzas-covarianzas robustos frente a heterocedasticidad. El paquete sandwich contiene todos los que podemos necesitar. En particular, vcovHC proporciona el estimador de White, robusto con heterocedasticidad y vcovBS una estimación bootstrap de los errores estándar. En el caso de que sospechemos que haya grupos heterogeneos se puede también emplear vcovCL (clustered covariance matrix estimation) espeficando la variable que contiene los clusters a emplear. En nuestro caso un candidato lógico habría sido la provincia o la CCAA, pero no está incluido en los datos.
La tabla de coeficientes con la inferencia robusta es la siguiente:
library(sandwich); library(lmtest)
coeftest(medad1,vcov=vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.494062 0.892935 -5.0329 5.427e-07
## EDAD 0.136486 0.015883 8.5934 < 2.2e-16
## ESTUD2 -0.742547 0.384062 -1.9334 0.05338
## ESTUD3 -0.431190 0.414601 -1.0400 0.29851
## ESTUDPAR2 0.630615 0.499135 1.2634 0.20664
## ESTUDPAR3 0.411476 0.509203 0.8081 0.41918
## TIPOUNIONPareja de hecho registrada 0.770176 1.304673 0.5903 0.55507
## TIPOUNIONPareja de hecho sin registrar 1.659943 0.315829 5.2558 1.691e-07
## NACIMOtro país -0.231274 0.500437 -0.4621 0.64405
## NACIMPAROtro país 2.377505 0.502785 4.7287 2.477e-06
##
## (Intercept) ***
## EDAD ***
## ESTUD2 .
## ESTUD3
## ESTUDPAR2
## ESTUDPAR3
## TIPOUNIONPareja de hecho registrada
## TIPOUNIONPareja de hecho sin registrar ***
## NACIMOtro país
## NACIMPAROtro país ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Podemos hacer también contrastes de significación de grupos con Anova del paquete car:
library(car)
Anova(medad1, vcov=vcovHC)Esta me marca que únicamente el tipo de unión, la edad y el nacimiento en el extranjero de la pareja son significativas a los niveles habituales.
El paquete sjPlot tiene la función tab_model que genera una tabla en formato html y que permite inferencia robusta y mostrar intervalos de confianza. La opción vcov.fun="HC" sirve para errores robustos frente a heterodasticidad y vcov.fun="CR" para robusta por conglomerados (requiere especificar cluster)
library(sjPlot)
tab_model(medad1,vcov.fun="HC",show.reflvl = TRUE, prefix.labels = "varname")|
Año de nacimineto de la pareja |
|||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | -4.49 | -6.25 – -2.74 | <0.001 |
| Edad en años cumplidos | 0.14 | 0.11 – 0.17 | <0.001 |
|
Factor de elevación del dato muestral a la población calibrado por sexo, grupo quinquenal de edad y provincia de residencia |
Reference | ||
| ESTUD: Primaria | Reference | ||
| ESTUD: Secundaria | -0.74 | -1.50 – 0.01 | 0.053 |
| ESTUD: Universitaria | -0.43 | -1.24 – 0.38 | 0.299 |
| ESTUDPAR: ESTUD: Primaria | Reference | ||
|
ESTUDPAR: ESTUD: Secundaria |
0.63 | -0.35 – 1.61 | 0.207 |
|
ESTUDPAR: ESTUD: Universitaria |
0.41 | -0.59 – 1.41 | 0.419 |
| TIPOUNION: Matrimonio | Reference | ||
|
TIPOUNION: Pareja de hecho registrada |
0.77 | -1.79 – 3.33 | 0.555 |
|
TIPOUNION: Pareja de hecho sin registrar |
1.66 | 1.04 – 2.28 | <0.001 |
| NACIM: España | Reference | ||
| NACIM: Otro paÃÂs | -0.23 | -1.21 – 0.75 | 0.644 |
| NACIMPAR: NACIM: España | Reference | ||
|
NACIMPAR: NACIM: Otro paÃÂs |
2.38 | 1.39 – 3.36 | <0.001 |
| Observations | 1479 | ||
| R2 / R2 adjusted | 0.101 / 0.095 | ||
Funciona razonablemente bien utilizando etiquetas como vemos en las variables de nivel de estudios. Hay problemas con las fuentes en el html (no en el panel Viewer, al menos en mi caso) y alguna cosa rara como que aparezca una mención a FACTOR como variable. Una viñeta muestra las opciones disponibles.
El modelo 1 puede ser demasiado simple (por ejemplo, puede haber no linealidad en la relación con EDAD o interacciones entre variables que no hemos considerado) o demasiado complejo (estamos introduciendo variables que no podemos descartar que no tengan efecto. Esto también está ligado a la muestra pequeña). Podéis explorar esas dimensiones como ejercicio. Referencias que os pueden ser útiles para explorar más sobre modelos lineales son Aguilar (2020), Hanck y otros (2020) y Ortega(2018).
Extensiones al modelo lineal: modelo GAM
Una de las limitaciones del modelo anterior es especificar una función lineal para la relación entre la edad del ego y la diferencia de edades. Esta relación podría ser no lineal. Existen diversas maneras de acomodar esto. Podríamos estimar modelos lineales que utilicen una función no lineal para la edad, como un polinomio o un spline. Otra alternativa que hace básicamente lo mismo (utiliza distintos tipos de splines) pero que no requiere especificar a priori el tipo de función que se utiliza son los modelos GAM (Clark, 2019)
library(mgcv)
medad2=gam(
DIFEDAD~s(EDAD)+ESTUD+ESTUDPAR+TIPOUNION+NACIM+NACIMPAR,
data=EFedad,weights=EFedad$PESO)
summary(medad2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## DIFEDAD ~ s(EDAD) + ESTUD + ESTUDPAR + TIPOUNION + NACIM + NACIMPAR
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.0271 0.3666 2.802 0.00515 **
## ESTUD2 -0.7425 0.3233 -2.297 0.02175 *
## ESTUD3 -0.4312 0.3766 -1.145 0.25240
## ESTUDPAR2 0.6306 0.3662 1.722 0.08526 .
## ESTUDPAR3 0.4115 0.3952 1.041 0.29802
## TIPOUNIONPareja de hecho registrada 0.7702 0.6947 1.109 0.26774
## TIPOUNIONPareja de hecho sin registrar 1.6599 0.2655 6.252 5.31e-10 ***
## NACIMOtro país -0.2313 0.3930 -0.589 0.55627
## NACIMPAROtro país 2.3775 0.3559 6.680 3.39e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(EDAD) 1 1 101.9 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0952 Deviance explained = 10.1%
## GCV = 16.397 Scale est. = 16.286 n = 1479plot(medad2)
El gráfico de la función suave estimada parece lineal con lo que es posible que no mejore el modelo con una función suave. De hecho, si nos fijamos en los grados de libertad equivalentes asociados al componente s(EDAD) vemos que es solo 1. Esto quiere decir que de hecho es un modelo lineal. Una forma de comprobarlo es utilizando criterios de selección de modelos, sabiendo que el modelo 2 es idéntico al modelo 1
library(broom)
list(medad1=medad1,medad2=medad2) %>% map_df(glance,.id="modelo")Vemos que son indistinguibles: no hay evidencia de no linealidad en la relación entre edad y diferencia de edad.
Sin embargo, el modelo sigue sin captar bien la varianza o la distribución. Para ello podemos utilizar modelos más complejos que permiten modelizar la varianza o momentos sucesivos. Dos paquetes en esta dirección son VGAM y GAMLSS.
Modelo GAMLSS
Los modelos GAMLSS son una generalización de los modelos GAM en dos direcciones:
- En vez de limitarnos a una familia de distribuciones, la exponencial, hay una variedad casi infinita de distribuciones disponibles, tanto para variables continuas definidas en distintos rangos, como discretas, mixturas,… que vienen definidas por una serie de momentos. Podemos obtener una demostración con este código
# install.packages("gamlss.demo")
library(gamlss.demo)
gamlss.demo()- Para cada momento podemos utilizar la fórmula que deseemos de modo que captemos la dependencia de ese momento respecto a las variables que nos interesen. Las distintas distribuciones han sido recodificadas de modo que, en general, los parámetros se correspondan con los momentos de orden 1 (esperanza), 2 (varianza), 3 (asimetría) y 4 (kurtosis)
Existe también la posibilidad de buscar qué distribución se ajusta mejor dentro de las disponibles para un problema con una variable. Vamos a hacerlo para nuestra diferencia de edades utilizando las distribuciones continuas con valores en todo el rango real.
library(gamlss)
distedad=fitDist(EFedad$DIFEDAD,type="realline")
# Para ver cuáles han sido los mejores ajustes
distedad$fits## SEP2 SHASHo SHASHo2 EGB2 JSU JSUo SEP1 SEP3
## 8105.666 8105.877 8105.877 8106.234 8106.671 8106.671 8109.049 8109.956
## SHASH ST3 SST ST2 ST1 ST5 SEP4 ST4
## 8111.730 8113.021 8113.021 8115.072 8116.092 8116.950 8117.949 8124.549
## GT TF2 TF PE PE2 NET LO exGAUS
## 8142.499 8143.161 8143.161 8155.027 8155.027 8169.108 8212.991 8270.467
## SN1 SN2 NO RG GU
## 8339.870 8362.638 8408.370 8649.384 9219.532En este caso la mejor distribución es la SEP2: Skew Exponential Power type 2.
Podemos ver cómo es el ajuste
distedad##
## Family: c("SEP2", "Skew Exponential Power type 2")
## Fitting method: "nlminb"
##
## Call: gamlssML(formula = y, family = DIST[i])
##
## Mu Coefficients:
## [1] 0.2342
## Sigma Coefficients:
## [1] 1.111
## Nu Coefficients:
## [1] 0.3734
## Tau Coefficients:
## [1] -0.03382
##
## Degrees of Freedom for the fit: 4 Residual Deg. of Freedom 1475
## Global Deviance: 8097.67
## AIC: 8105.67
## SBC: 8126.86wp(distedad)
Vemos que es una distribución que tiene los cuatro momentos definidos. El diagrama de gusano (worm plot) que representa un qqplot al que le quitamos la tendencia, muestra que el ajuste es correcto.
Vamos a pasar ahora a estimar un modelo donde los cuatro momentos pueden depender de las variables más importantes que vimos anteriormente: EDAD (con posible perfil no lineal), TIPOUNION y NACIMPAREJA. Posiblemente es un modelo demasiado complejo, pero eso se puede corregir buscando los modelos con el AIC o el BIC más bajos. Un problema es que no permite NAs en los datos, así que tenemos que crear antes otro objeto de datos sin NA
EFedad2=EFedad %>% dplyr::select(DIFEDAD,EDAD,TIPOUNION,NACIMPAR,PESO) %>%
na.omit()
medad3=gamlss(DIFEDAD~pb(EDAD)+TIPOUNION+NACIMPAR,
sigma.fo =~EDAD+TIPOUNION+NACIMPAR,
nu.fo = ~EDAD+TIPOUNION+NACIMPAR,
tau.fo =~EDAD+TIPOUNION+NACIMPAR,
data=EFedad2, weight=PESO,family=SEP2)## GAMLSS-RS iteration 1: Global Deviance = 8083.419
## GAMLSS-RS iteration 2: Global Deviance = 8062.818
## GAMLSS-RS iteration 3: Global Deviance = 8046.238
## GAMLSS-RS iteration 4: Global Deviance = 8033.658
## GAMLSS-RS iteration 5: Global Deviance = 8024.598
## GAMLSS-RS iteration 6: Global Deviance = 8019.385
## GAMLSS-RS iteration 7: Global Deviance = 8016.557
## GAMLSS-RS iteration 8: Global Deviance = 8013.947
## GAMLSS-RS iteration 9: Global Deviance = 8011.385
## GAMLSS-RS iteration 10: Global Deviance = 8009.96
## GAMLSS-RS iteration 11: Global Deviance = 8008.181
## GAMLSS-RS iteration 12: Global Deviance = 8006.669
## GAMLSS-RS iteration 13: Global Deviance = 8005.723
## GAMLSS-RS iteration 14: Global Deviance = 8005.239
## GAMLSS-RS iteration 15: Global Deviance = 8004.287
## GAMLSS-RS iteration 16: Global Deviance = 8003.451
## GAMLSS-RS iteration 17: Global Deviance = 8002.075
## GAMLSS-RS iteration 18: Global Deviance = 8001.449
## GAMLSS-RS iteration 19: Global Deviance = 8001.118
## GAMLSS-RS iteration 20: Global Deviance = 8001.038# Como no ha convergido, le decimos que siga
medad3=gamlss(DIFEDAD~pb(EDAD)+TIPOUNION+NACIMPAR,
sigma.fo =~EDAD+TIPOUNION+NACIMPAR,
nu.fo = ~EDAD+TIPOUNION+NACIMPAR,
tau.fo =~EDAD+TIPOUNION+NACIMPAR,
data=EFedad2, weight=PESO,family=SEP2,
start.from=medad3)## GAMLSS-RS iteration 1: Global Deviance = 7999.819
## GAMLSS-RS iteration 2: Global Deviance = 7999.93
## GAMLSS-RS iteration 3: Global Deviance = 8000.592
## GAMLSS-RS iteration 4: Global Deviance = 8000.542
## GAMLSS-RS iteration 5: Global Deviance = 8000.515
## GAMLSS-RS iteration 6: Global Deviance = 8000.464
## GAMLSS-RS iteration 7: Global Deviance = 8000.45
## GAMLSS-RS iteration 8: Global Deviance = 8000.443
## GAMLSS-RS iteration 9: Global Deviance = 8000.418
## GAMLSS-RS iteration 10: Global Deviance = 8000.416
## GAMLSS-RS iteration 11: Global Deviance = 8000.414
## GAMLSS-RS iteration 12: Global Deviance = 8000.414summary(medad3)## ******************************************************************
## Family: c("SEP2", "Skew Exponential Power type 2")
##
## Call:
## gamlss(formula = DIFEDAD ~ pb(EDAD) + TIPOUNION + NACIMPAR, sigma.formula = ~EDAD +
## TIPOUNION + NACIMPAR, nu.formula = ~EDAD + TIPOUNION + NACIMPAR,
## tau.formula = ~EDAD + TIPOUNION + NACIMPAR, family = SEP2,
## data = EFedad2, weights = PESO, start.from = medad3)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: identity
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.072211 0.388588 5.333 1.12e-07
## pb(EDAD) -0.045538 0.008934 -5.097 3.89e-07
## TIPOUNIONPareja de hecho registrada -0.311228 0.634081 -0.491 0.623617
## TIPOUNIONPareja de hecho sin registrar 0.540958 0.193541 2.795 0.005256
## NACIMPAROtro país 1.129779 0.314965 3.587 0.000345
##
## (Intercept) ***
## pb(EDAD) ***
## TIPOUNIONPareja de hecho registrada
## TIPOUNIONPareja de hecho sin registrar **
## NACIMPAROtro país ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.366719 0.130570 -2.809 0.00504
## EDAD 0.032666 0.002822 11.577 < 2e-16
## TIPOUNIONPareja de hecho registrada 0.454616 0.177174 2.566 0.01039
## TIPOUNIONPareja de hecho sin registrar 0.340773 0.055406 6.150 9.93e-10
## NACIMPAROtro país 0.520479 0.070637 7.368 2.86e-13
##
## (Intercept) **
## EDAD ***
## TIPOUNIONPareja de hecho registrada *
## TIPOUNIONPareja de hecho sin registrar ***
## NACIMPAROtro país ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.663839 0.167912 -9.909 <2e-16
## EDAD 0.048288 0.003857 12.520 <2e-16
## TIPOUNIONPareja de hecho registrada 0.014732 0.196394 0.075 0.940
## TIPOUNIONPareja de hecho sin registrar 0.115329 0.070502 1.636 0.102
## NACIMPAROtro país 0.080202 0.091119 0.880 0.379
##
## (Intercept) ***
## EDAD ***
## TIPOUNIONPareja de hecho registrada
## TIPOUNIONPareja de hecho sin registrar
## NACIMPAROtro país
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Tau link function: log
## Tau Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.268942 0.190574 -1.411 0.1584
## EDAD 0.009440 0.004215 2.240 0.0253 *
## TIPOUNIONPareja de hecho registrada -0.118468 0.290332 -0.408 0.6833
## TIPOUNIONPareja de hecho sin registrar 0.102288 0.081017 1.263 0.2069
## NACIMPAROtro país 0.045316 0.116230 0.390 0.6967
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms may not be reliable.
## ------------------------------------------------------------------
## No. of observations in the fit: 1479
## Degrees of Freedom for the fit: 20.00038
## Residual Deg. of Freedom: 1459
## at cycle: 12
##
## Global Deviance: 8000.414
## AIC: 8040.414
## SBC: 8146.399
## ******************************************************************plot(medad3)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = -0.0088528
## variance = 0.9879008
## coef. of skewness = -0.07594307
## coef. of kurtosis = 3.061939
## Filliben correlation coefficient = 0.9995814
## ******************************************************************wp(medad3)
El modelo consigue un buen ajuste a la distribución pero tiene posiblemente más parámetros de los necesarios. Quitamos las variables distintas de la edad de los momentos de orden 3 y 4
medad4=gamlss(DIFEDAD~pb(EDAD)+TIPOUNION+NACIMPAR,
sigma.fo =~EDAD+TIPOUNION+NACIMPAR,
nu.fo = ~EDAD,
tau.fo =~EDAD,
data=EFedad2, weight=PESO,family=SEP2,
start.from=medad3)## GAMLSS-RS iteration 1: Global Deviance = 8005.854
## GAMLSS-RS iteration 2: Global Deviance = 8002.187
## GAMLSS-RS iteration 3: Global Deviance = 8002.099
## GAMLSS-RS iteration 4: Global Deviance = 8002.088
## GAMLSS-RS iteration 5: Global Deviance = 8002.053
## GAMLSS-RS iteration 6: Global Deviance = 8001.939
## GAMLSS-RS iteration 7: Global Deviance = 8001.93
## GAMLSS-RS iteration 8: Global Deviance = 8001.93summary(medad4)## ******************************************************************
## Family: c("SEP2", "Skew Exponential Power type 2")
##
## Call: gamlss(formula = DIFEDAD ~ pb(EDAD) + TIPOUNION + NACIMPAR,
## sigma.formula = ~EDAD + TIPOUNION + NACIMPAR, nu.formula = ~EDAD,
## tau.formula = ~EDAD, family = SEP2, data = EFedad2,
## weights = PESO, start.from = medad3)
##
## Fitting method: RS()
##
## ------------------------------------------------------------------
## Mu link function: identity
## Mu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.65601 0.08879 18.650 < 2e-16 ***
## pb(EDAD) -0.03920 0.00378 -10.371 < 2e-16 ***
## TIPOUNIONPareja de hecho registrada -0.08728 1.83492 -0.048 0.962
## TIPOUNIONPareja de hecho sin registrar 0.80576 0.09045 8.909 < 2e-16 ***
## NACIMPAROtro país 1.34504 0.17864 7.529 8.88e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Sigma link function: log
## Sigma Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.27052 0.15254 -1.773 0.0764 .
## EDAD 0.03112 0.00350 8.892 < 2e-16 ***
## TIPOUNIONPareja de hecho registrada 0.48073 0.15998 3.005 0.0027 **
## TIPOUNIONPareja de hecho sin registrar 0.27753 0.05824 4.765 2.07e-06 ***
## NACIMPAROtro país 0.49595 0.06747 7.351 3.26e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Nu link function: identity
## Nu Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.509887 0.150645 -10.02 <2e-16 ***
## EDAD 0.046183 0.004292 10.76 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------------------------------
## Tau link function: log
## Tau Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.119425 0.178606 -0.669 0.504
## EDAD 0.006874 0.004346 1.582 0.114
##
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas:
## i) Std. Error for smoothers are for the linear effect only.
## ii) Std. Error for the linear terms maybe are not accurate.
## ------------------------------------------------------------------
## No. of observations in the fit: 1479
## Degrees of Freedom for the fit: 14.00038
## Residual Deg. of Freedom: 1465
## at cycle: 8
##
## Global Deviance: 8001.93
## AIC: 8029.931
## SBC: 8104.12
## ******************************************************************plot(medad4)
## ******************************************************************
## Summary of the Quantile Residuals
## mean = -0.008679355
## variance = 0.988063
## coef. of skewness = -0.09147764
## coef. of kurtosis = 3.080067
## Filliben correlation coefficient = 0.9994709
## ******************************************************************wp(medad4)
Conclusiones
En una sesión es imposible cubrir en detalle todos los aspectos de la modelización de encuestas. Lo que hemos visto ha sido una selección de los aspectos que emergen en cualquier análisis empírico, desde los más prosaicos sobre cómo importar los datos, comprobar las variables, etc, hasta los más avanzados como pueden ser los modelos GAMLSS que abren nuevas fronteras a la modelización al caracterizar la dependencia de todos los momentos de la distribución respecto a las variables de interés. Podemos apreciar la distinta forma de confrontar el análisis de los estimadores robustos de varianzas-covarianzas o los modelos gamlss. Mientras que los primeros solo pretender no fallar por la posible existencia de heterocedasticidad, que no se analiza, los segundos se basan en captar toda la estructura de heterocedasticidad y de dependencia de los momentos sucesivos en los datos. El ánimo no ha sido el de ser exhaustivo sino mostrar cosas que se pueden hacer y que pueden ser útiles para nuestra investigación.
Respecto a los aspectos que creo interesantes y que no hemos tocado están los siguientes:
La modelización de la variabilidad a través de modelos mixtos de efectos fijos y aleatorios. Esta es una dimensión más que se puede integrar dentro de los modelos. Una referencia es Lüdecke (2019). No lo hemos cubierto porque la unidad natural para la modelización en nuestra encuesta serían las provincias, pero, dado el pequeño tamaño muestral, no se incluye esa información en la muestra. Sí que se incluye la CCAA para la muestra de mujeres.
Modelos para variables de recuento y discretas. Aunque no las hemos visto,
glm, análogo alm, sirve para los modelos básicos. Las dos extensiones que hemos tratado,gamygamlss, pueden acomodar también modelos para variables de estos tipos.El INE no proporciona información sobre las unidades de muestreo. Esto es importante para estimar mejor los errores estándar de los parámetros, pero no es posible. El próximo día estudiaremos una encuesta parecida a la de Fecundidad que sí que incluye esta información, la encuesta Demográfica y de Salud (DHS) de Perú. Esto nos permitirá ver en qué manera influye la estructura compleja de muestreo sobre los errores estándar. Lo más (y lo menos) que podíamos hacer en el caso de hoy es utilizar los pesos para tener en cuenta las distintas probabilidades de selección. Es también importante la advertencia de redefinir unos pesos estandarizados puesto que no todos los paquetes interpretan los pesos del mismo modo, y queremos evitar que se interpreten como frecuencias. Una forma de salir de dudas es estimar el mismo modelo utilizando los dos vectores de pesos: el factor de elevación y el peso estandarizado. Si son distintos, en general mucho menores los basados en el factor de elevación, es que se da este problema.
Referencias
Aguilar Esteva, A. (2020) Notas de Microeconometría Aplicada
Amat Rodrigo, J. (2020) GAMLSS: modelos aditivos generalizados para posición, escala y forma
Clark, M. (2019) Generalized Additive Models
Corcoran, D. (2020) Manipulación de datos e investigación reproducible en R
Fox, J. y S. Weisberg (2019) An R Companion to Applied Regression, 3ª Ed, Sage
González-Vidal, A., Maurandi-López, A., y Del Río, L. (2019) tabularCola - Tabulación de datos con R: Curso on line autónomo (Version 1.5), Zenodo.
Grolemand, G. y Wickham, H. (2020) R for data science, O’Reilly. Traducción en español.
Hanck, C. et al (2020) Introduction to econometrics with R.
Lüdecke, D. (2019) Mixed models snippets
Ortega, J. A. (2018) ModelaR: Modelización avanzada con R, Universidad de Salamanca.
Ortega, J. A. (2020) Instalación y primeros pasos con R y R Studio. EncuestaR