Importación

#### MUY IMPORTANTE ####
#### Modificar el path relativo a la carpeta del proyecto en su computadora ####

#base de hogares
hogares <- read_sav("../../static/docs/datos/ENES_Hogares_version_final.sav")

#base de personas
personas <- read_sav("../../static/docs/datos/ENES_Personas_version_final.sav")

Diccionario de variables

Los nombres de las variables son cortos y fáciles de manejar en R o cualquier otro software en el que usamos un script o consola como forma de acceder a los datos. Sin embargo a menos que estemos muy familiarizados con estos nombres será difícil recordar a que pregunta del cuestionario corresponden. Lo mismo si queremos comunicar los resultados de un análisis de estos datos. Es buena idea crear un diccionario de variables: una estructura de datos relacional que vincule a los nombres cortos de las variables con su definición larga. Para eso se crean dos diccionarios de variables, uno por cada base. Al ser objetos de la clase data.frame podemos verlos y explorarlos con la función View() o directamente usarlos en el código para renombrar las variables.

La función renombrar() utiliza al diccionario de variables para cambiar los nombres cortos por los nombres largos cuando es necesario. Esta función recibe como input un data.frame y regresa otro data.frame en el que solo cambian los nombres. Por lo tanto no es necesario iterarla explícitamente a todas las columnas y funciona sin problemas con el operador %>%.

diccionario_hogares <- tibble(nombre = names(hogares), 
                                etiqueta = map_chr(hogares, attr, "label"))
diccionario_personas  <- tibble(nombre = names(personas), 
                                etiqueta = map_chr(personas, attr, "label"))
diccionario <- bind_rows(diccionario_hogares, diccionario_personas)

# Función renombrar
# Requiere nombres fijos de diccionario y columnas!

renombrar <- function(x) {
  rename_all(x, ~diccionario$etiqueta[match(., diccionario$nombre)])
}

# Para ver el diccionario quite el símbolo # de la siguiente línea y ejecútela
# View(diccionario)

Unión de las bases de datos

Por diseño los datos se ubican en dos bases, una con atributos de los hogares y otra con atributos de las personas. Dado que todas las personas encuestadas pertenecen a alguno de los hogares y que los códigos de vivienda y hogar están presentes en las dos bases de datos es posible hacer un join y reunir toda la información.

Obtendremos una base de datos en la que cada fila será una persona (n=27610). Evitar confusiones con la unidad de análisis en cada caso corre por cuenta de quién analice la base de datos. Valga recordar que no es lo mismo habitar una vivienda con techo de chapas que tener techo de chapas.

unida <- left_join(personas, hogares, by = c("nocues", "nhog")) %>% 
  rename(fcalib_hog = f_calib3.y, 
         fcalib_per = f_calib3.x)

## Warning: Column `nocues` has different attributes on LHS and RHS of join

## Warning: Column `nhog` has different attributes on LHS and RHS of join

Proporción de personas entrevistadas por sexo

La primera tabla que realizaremos utiliza información sobre las personas encuestadas. Dicho de otro modo, no buscamos inferir un parámetro de la población de la que proviene la muestra ni utilizamos el ponderador muestral para corregir sesgos de selección por diseño o problemas de no respuesta. Por lo tanto las proporciones serán proporciones de personas a las que se aplicó la entrevista. Más adelante buscaremos soluciones para ambos problemas.

Usando la función count() obtenemos tablas de frecuencias en forma “largo”: una o más columnas indican el grupo sobre el que se calculó la frecuencia y la columna n indica la frecuencia calculada. En este caso (univariado) este formato es obvio al punto de parecer natural. Sin embargo veremos que para tablas de más de una dimensión diverge del formato más habitual de la tabla de contingencia. De todos modos podemos pasarlo a a tabla de contingencia en caso de ser necesario.

personas %>% 
  select(v109) %>% 
  as_factor() %>%    #La paso desde numérica etiquetada a factor
  count(v109) %>%    #Indico la variable a contar
  mutate(prop = n / sum(n))  #Agrego una columna con las proporciones

## # A tibble: 3 x 3
##   v109      n     prop
##   <fct> <int>    <dbl>
## 1 Varón 13258 0.480   
## 2 Mujer 14346 0.520   
## 3 Otro      6 0.000217

Proporción de entrevistas por sexo y estado civil

En este caso obtenemos las frecuencias condicionales de dos variables: sexo y una pregunta sobre si la persona entrevistada planchó o no planchó ropa la semana pasada.

Modificamos levementa la sintaxis anterior:

1. Renombramos las variables para aprovechar las etiquetas largas
2. En lugar de definir al grupo para las frecuencias en las llamada misma a count() lo definimos previamente con group_by_all(). Este atajo permite definir grupos para todas las variables disponibles. NUNCA usarlo sin hacer previamente un subconjunto de variables con select().
3. Pasamos la tabla a formato “ancho”, ubicando a las categorías de una variable como filas y las categorías de la otra como columnas.
4. Usamos la función kable() para mejorar el formato de salida de la tabla.

personas %>% 
  select(v109, v154) %>% 
  renombrar() %>% 
  as_factor() %>% 
  group_by_all() %>% 
  count() %>% 
  spread(Sexo, n, fill = 0) %>% 
  knitr::kable(caption = "Tareas del hogar y género")

## Warning: Factor `Planchar durante la semana pasada` contains implicit NA,
## consider using `forcats::fct_explicit_na`

## Warning: Factor `Planchar durante la semana pasada` contains implicit NA,
## consider using `forcats::fct_explicit_na`

## Warning: Factor `Planchar durante la semana pasada` contains implicit NA,
## consider using `forcats::fct_explicit_na`

Los resultados, lamentablemente, sorprenden poco. Más allá de la interpretación de la tabla ha emergido un dato técnico de las variables. La función count() emite unas advertencias algo extrañas:

1: Factor `Planchar durante la semana pasada` contains implicit NA, consider using `forcats::fct_explicit_na` 

Simplemente nos advierte que una de las variables (Planchar durante la semana pasada) contiene casos perdidos (missing), que en R se codifican como NA. No sabemos mucho sobre el orígen de esos missing, pero deberemos tomar una decisión: los consideramos una categoría de respuesta en sí misma (hacerlos explícitos) o los eliminamos de la base de datos.²

Seguiremos la práctica utilizada por autores y autoras del libro “Cómo somos, vivimos y convivimos en una sociedad desigual” y eliminaremos los missing con drop_na().

personas %>% 
  select(v109, v154) %>% 
  renombrar() %>% 
  as_factor() %>% 
  drop_na() %>% 
  group_by_all() %>% 
  count() %>% 
  
  #Pasar a formato ancho para facilitar la lectura
  
  spread(key = Sexo, value = n, fill = 0) %>% 
  knitr::kable(caption = "Tareas del hogar y género sin missing")

Proporciones de conocimiento de idioma extranjero - Nivel Nacional

Este es el caso más simple, ya que calcularemos las frecuencias de una sola variable. Sin embargo es necesario incluir otras dos en los datos que utilizaremos: la edad para poder aplicar el filtro de mayores de 29 años y f_calib3, la columna del data.frame personas en la que se ubica el ponderador muestral.

Procedimiento:

• seleccionar las variables de interés
• filtrar las filas con edades mayores o iguales a 30 años
• eliminar los missing de la base de datos
• recodificar la variable v118a compactando a las categorías “Sí, nivel avanzado” y “Sí, bilingüe” en la nueva categoría “Avanzado o bilingüe”. Para mantener los nombres idénticos a la tabla 7.9 recodificamos “Sin conocimiento” como “No”.
• contar utilizando el ponderador muestral como argumento wt= en la función count()
• Calcular proporciones por cada nivel conocimiento de idioma extranjero

personas %>% 
  select(v118a, f_calib3, v108) %>% 
  filter(v108 >= 30) %>% 
  drop_na %>% 
  as_factor() %>% 
  mutate(v118a = fct_recode(v118a, "Avanzado o bilingüe" = "Sí, nivel avanzado", 
                            "Avanzado o bilingüe" = "Sí, bilingüe",
                            "Sin conocimiento" = "No")) %>% 
  count(v118a, wt = f_calib3) %>% 
  mutate (Nacional = n/sum(n)) %>% 
  select(-n) -> Nacional

Nacional

## # A tibble: 4 x 2
##   v118a               Nacional
##   <fct>                  <dbl>
## 1 Sin conocimiento      0.727 
## 2 Sí, nivel básico      0.178 
## 3 Sí, nivel medio       0.0505
## 4 Avanzado o bilingüe   0.0450

Proporciones de conocimiento de idioma extranjero por regiones

En este caso calculamos la proporción de personas por categoría de conocimiento de idioma extranjero dentro de cada región.

El procedimiento es el mismo que en el anterior, excepto que que utilizamos a las dos variables en group_by() previo al conteo y luego definimos un grupo por región para el cálculo de las proporciones.

spread() se encarga de pasar del formato “largo” de tríos de claves y valores al formato propiamente tabular.

personas %>% 
  select(v118a, region, f_calib3, v108) %>% 
  filter(v108 >= 30) %>% 
  drop_na %>% 
  as_factor() %>% 
  mutate(v118a = fct_recode(v118a, "Avanzado o bilingüe" = "Sí, nivel avanzado", 
                            "Avanzado o bilingüe" = "Sí, bilingüe",
                            "Sin conocimiento" = "No")) %>% 
  group_by(v118a, region) %>%
  count(wt = f_calib3) %>% 
  group_by(region) %>% 
  mutate(prop = n/sum(n)) %>% 
  select(-n) %>% 
  spread(region, prop) %>% 
  #Una expresión regular para eliminar el texto entre paréntesis de los nombres de variable
  rename_all(str_remove, " \\(([^\\)]+)\\)") -> por_region

por_region

## # A tibble: 4 x 8
##   v118a     `Gran Buenos Ai…   Cuyo Pampeana Centro    NEA    NOA Patagonia
##   <fct>                <dbl>  <dbl>    <dbl>  <dbl>  <dbl>  <dbl>     <dbl>
## 1 Sin cono…           0.691  0.792    0.704  0.712  0.797  0.805     0.760 
## 2 Sí, nive…           0.173  0.148    0.200  0.216  0.125  0.151     0.161 
## 3 Sí, nive…           0.0697 0.0294   0.0549 0.0398 0.0410 0.0280    0.0363
## 4 Avanzado…           0.0667 0.0304   0.0415 0.0327 0.0373 0.0165    0.0430

Proporciones de conocimiento de idioma extranjeo (CABA)

El caso de CABA requiere que filtremos la base de datos y dejemos solamente aquellas filas en las que el aglomerado es CABA.

personas %>% 
  select(v118a, aglo, f_calib3, v108) %>% 
  filter(v108>=30) %>% 
  as_factor() %>% 
  filter(aglo == "CABA") %>% 
  mutate(v118a = fct_recode(v118a, "Avanzado o bilingüe" = "Sí, nivel avanzado", 
                            "Avanzado o bilingüe" = "Sí, bilingüe",
                            "Sin conocimiento" = "No")) %>% 
  count(v118a, wt = f_calib3) %>% 
  mutate (CABA = n/sum(n)) %>% 
  select(-n) -> CABA

Unir las tablas

Como las tres tablas tienen en común la columna v118a es posible hacer un join para unir las tres estructuras en una sola.

list(CABA, por_region,  Nacional) %>% 
  reduce(left_join) -> tabla7.9

## Joining, by = "v118a"
## Joining, by = "v118a"

## Warning: Column `v118a` has different attributes on LHS and RHS of join

tabla7.9 %>% 
  mutate_if(is.numeric, scales::percent) %>% 
  knitr::kable(caption = "Conocimiento de idiomas extranjeros en la población de 30 años o más (en porcentajes). Total nacional, CABA y regiones (Replicación)")

El resultado replica exactamente el reportado en la tabla 7.9 del libro.

Gráfico de proporciones

Los datos de la Tabla 7.9 se pueden representar facilmente en un gráfico de barras apiladas. Se utiliza una paleta de colores de alto contraste y apta para personas con la mayoría de los tipos de daltonismo.

tabla7.9 %>% 
  gather(variable, valor, -v118a) %>% 
  mutate(variable = fct_relevel(variable, "Nacional"), 
         variable = fct_recode(variable, "GBA" = "Gran Buenos Aires")) %>% 
  ggplot(aes(x=variable, y = valor, fill = v118a)) + 
  geom_col() + 
  theme_minimal() + 
  labs(fill = "Lengua extranjera", 
       y = NULL, 
       x = "Región", 
       title = "Concentración de ventajas") + 
  theme(#axis.text.y = element_blank(), 
        legend.position = "bottom") + 
  scale_fill_viridis_d() + 
  scale_y_continuous(labels = scales::percent)

Bootstrap de proporciones de conocimiento idioma extranjero por sexo

El primer caso que tomaremos es el de la distribución de conocimiento de idioma extranjero por sexo. A modo de referencia comenzaremos calculando las proporciones por grupos con los datos tal como están en la base. Es decir, calculando es estimador puntual de proporciones.

# Preparación de los datos

personas %>% 
  select(v118a, region, f_calib3, v108, v109) %>% 
  filter(v108 >= 30) %>% 
  drop_na %>% 
  as_factor() %>% 
  mutate(v118a = fct_recode(v118a, "Avanzado o bilingüe" = "Sí, nivel avanzado", 
                            "Avanzado o bilingüe" = "Sí, bilingüe",
                            "Sin conocimiento" = "No")) %>% 
  select(-v108) -> para_bootstrap

para_bootstrap %>% 
  filter ( v109 != "Otro") %>% 
  group_by(v109) %>%
  count(v118a, v109, wt = f_calib3) %>% 
  group_by(v109) %>% 
  mutate(prop = n/sum(n), 
         prop = scales::percent(prop)) %>% 
  select(-n) %>% 
  spread(v109, prop) %>% 
  knitr::kable(caption = "Estimador puntual de proporciones de conocimiento de idioma extranjero por sexo")

Una lectura precipitada de esta tabla nos llevaria a pensar que, aunque pequeñas, hay diferencias entre las proporciones de hombres y mujeres en todos los niveles de conocimiento de idioma extranjero. Esto es verdadero para la muestra, pero no sabemos si lo es para la población de la que seleccionamos la muestra. Si el azar hubiera seleccionado otros hogares bien podríamos tener otras proporciones. El método bootstrap nos ayuda a cuantificar la probabilidad de que en una muestra al azar obtengamos una diferencia como la que estamos observando en la tabla. Este procedimiento es intrínsecamente incierto, por eso prefiero hablar de cuantificar la incertidumbre más que poner a prueba una hipótesis. La afirmación que hagamos sobre la diferencia de proporciones depende de los números que da el bootstrap, pero sobre todo de nuestra tolerancia a la incertidumbre.

# Replicaciones bootstrap

rerun(10000, 
para_bootstrap %>% 
          sample_frac(size = 1,     
                      replace = TRUE) %>% 
  group_by(v109) %>%
  count(v118a, v109, wt = f_calib3) %>% 
  mutate(prop = n /sum(n)) ) -> boot_idiomas_sexo

Con las replicaciones realizadas podemos explorar la distribución de las proporciones. Visualmente un histograma presenta nos permite dar cuenta de la variabilidad de las proporciones estimadas en las replicaciones, así como valorar visualmente si hay una diferencia significativa de en las proporciones poblacionales en alguna de las categorías.

Nótese que cada uno de los subgráficos tiene su propia escala en el eje \(x\). Se omite la etiqueta del eje \(y\) por tener poca relevancia. Su escala depende de una decisión arbitraria: el número de replicaciones.

boot_idiomas_sexo %>% 
  bind_rows() %>% 
        filter(v109 != "Otro") %>% 
    ggplot(aes(x=prop, fill = v109)) +
      geom_histogram(position = "dodge", bins = 150) + 
      facet_wrap(~v118a, scales = "free") +  
      theme_minimal() + 
  theme(legend.position = "bottom", 
        axis.text.y = element_blank()) + 
  scale_x_continuous(labels = scales::percent) + 
  scale_fill_viridis_d() +
  labs(title = "Distribución bootstrap de proporciones", 
       fill = "Sexo", 
       y = NULL, 
       x = NULL)

En todos los casos encontramos una distribución simétrica y muy concentrada en torno a la media. Podemos esperar intervalos de confianza estrechos.

Estos gráficos tienen una interpretación muy directa. Si las “campanas” de barras amarillas y violetas estuvieran completamente superpuestas interpretaríamos que no hay diferencias en la proporción de conocimiento de idioma extranjero entre varones y mujeres.

Si estuvieran completamente separadas interpretaríamos que hay una diferencia en las proporciones a nivel poblacional. Es decir, que la diferencia de proporciones en el estimador puntual no es un caso aislado, producto del azar que nos dió una mano mala. Sacando casos muy poco frecuentes o con datos simulados lo usual es que haya un poco de superposición entre las campanas: esa región entre los “picos” de las campanas en la que hay barras amarillas y violetas. Tal es el caso de las categorías “Sin conocimiento”, “Sí, nivel básico” y “Sí, nivel intermedio”. En muchas de las muestras replicadas por bootstrap la proporción de hombres es mayor que la de mujeres y en muchas otras obtenemos opuesto resultado. No podemos afirmar que no haya una diferencia, pero tampoco que la haya. Al menos no con el intrumento con el que hemos medido. Cualquiera de las afirmaciones tendría un nivel inaceptablemente alto de incertidumbre y en esos casos la costrumbre es aceptar la hipótesis de nulidad: no hay diferencia.

El caso de “Avanzado o bilingüe” es diferente. Hay una zona en la que las barras se superponen, pero es pequeña. Sin omitir que existe una probabilidad mayor que 0 de que en la población las proporciones sean iguales es muy razonable afirmar que hay más mujeres que hombres con conocimientos avanzados de un idioma extranjero o que son directamente bilingües. Lo sabemos porque en muy pocas replicaciones obtuvimos el resultado inverso.

boot_idiomas_sexo %>% 
  bind_rows() %>% 
      filter(v109 != "Otro") %>% 
  group_by(v118a, v109) %>% 
  summarise(perc1 = quantile(prop, probs = 0.025), 
            perc99  = quantile(prop, probs = 0.975))  %>% 
  mutate_at(vars(perc1, perc99), scales::percent) %>% 
  unite(intervalo, perc1, perc99, sep = " - ") %>% 
  spread(v109, intervalo) %>% 
  knitr::kable()

Conocimiento de idioma extranjero. Intervalos de confianza.

En este último caso indagamos sobre las proporciones de conocimiento de idioma extranjero por región. El procedimiento de replicación es el mismo, la diferencia en es que en lugar de presentar la información sobre la incertidumbre de la estimación en forma gráfica lo haremos a través de intervalos de confianza. Aprovechando la distribución aproximadamente normal del bootstrap, la forma más simple obtener estos intervalos es reportar los cuantiles 0.025 y 99.75 de las proporciones replicadas.

rerun(10000, 
para_bootstrap %>% 
          sample_frac(size = 1,     
                      replace = TRUE) %>% 
  count(v118a, region, wt = f_calib3) %>% 
  group_by(region) %>% 
  mutate(prop = n /sum(n)) ) -> boot_idiomas_region

boot_idiomas_region%>% 
  bind_rows() %>% 
  group_by(region, v118a) %>% 
    summarise(perc1 = quantile(prop, probs = 0.025), 
            perc99  = quantile(prop, probs = 0.975))  %>% 
  mutate_if(is.numeric, scales::percent) %>% 
  unite(intervalo, c(perc1, perc99), sep = " - ") %>% 
  spread(v118a, intervalo) %>% 
  ungroup() %>% 
  mutate(region = str_remove(region, " \\(([^\\)]+)\\)")) %>% 
  knitr::kable(caption = "Conocimiento de idioma extranjero por región. Intervalo de confianza al 99%")

## `mutate_if()` ignored the following grouping variables:
## Column `region`

Hay diferencias importantes en todas las categorías de conocimiento de idioma extrajero, sobre todo si tomamos como referencia al Gran Buenos Aires.⁵ Tiene la proporción más alta de personas con nivel avanzado o bilingüe y una de las mas bajas de personas en la categoría “Sin conocimientos”. Al mismo tiempo sabemos que en esa región los niveles educativos agregados son más altos que en las demás y que el nivel educativo a su vez está relacionado con el nivel de conocimiento de idioma extranjero. Quizas la diferencia que observamos en los niveles de conocimiento de idioma extranjero pueda ser atribuida a la distribución desigual de las personas por nivel educativo en las regiones. Para responder esta pregunta es necesario utilizar alguna herramienta más sofisticada que las tablas de contingencia. En la próxima entrada buscaremos una respuesta a esta pregunta ajustando un modelo de regresión lineal a los datos de las ENES.