Recodificación de variables categóricas ordinales

Los objetos de las clases factor y ordered se utilizan en R para manejar variables categóricas nominales y ordinales respectivamente. Se trata de objetos cuyo tratamiento es complicado por la suma de atributos que los componen (niveles, etiquetas, orden) y son una fuente de dolores de cabeza para usuarios principiantes y avanzados. Las funciones de la librería forcats::, incluída en tidyverse reducen algunos de estos problemas. Aquí presentamos una estrategia basada en la creación de una lista relacional de valores viejos y nuevos. Moviendo la definición de las relaciones de recodificación a una estructura de datos generamos un código más simple y fácil de mantener.¹

Cálculo de coeficientes de correlación para variables ordinales

Nada impide que analicemos la relación entre dos variables ordinales usando tablas de contingencia, sin embargo al hacerlo desperdiciamos un nivel adicional de información: el orden. Las buenas prácticas impiden tajantemente, en cambio, que analicemos esa misma relación calculando el coeficiente de correlación producto-momento entre los códigos numéricos de las variables. Como alternativa se utiliza el coeficiente de correlación policórica.

Visualización de datos ordinales

Una vez más, nada impide que utilicemos gráficos de barras apiladas de proporciones o gráficos de mosaico para analizar variables ordinales. Sin embargo podemos aprovechar el orden de la variable y presentar esas proporciones de manera acumulada.

Ajuste de modelos con \(y\) ordinal

Algunas preguntas pueden responderse con 2 variables y en ese caso nada mejor que gráficos o estadísticos simples. Cuando el número de variables aumenta debemos recurrir a herramientas más sofisticadas como los modelos lineales. Ubicando una variable dependiente ordinal en el modelo simplifica mucho la interpretación de los coeficientes en comparación a un modelo logit multinomial para \(y\) nominal.

library(haven)
library(tidyverse)
library(coefplot)

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")

unida <- left_join(personas, hogares, by = c("nocues", "nhog")) %>% 
  rename(fcalib_hog = f_calib3.y, 
         fcalib_per = f_calib3.x)
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)
renombrar <- function(x) {
  rename_all(x, ~diccionario$etiqueta[match(., diccionario$nombre)])
}

Niveles educativos por región (datos originales)

personas %>% 
  select(edad = v108, sexo = v109, region, t_aglo, v118a, nivel_ed, f_calib3) %>% 
  filter(edad >= 30) %>% 
  mutate_at(vars(-edad, -f_calib3), as_factor) %>% 
  mutate(region = str_remove(region, " \\(([^\\)]+)\\)"), 
         region = fct_recode(region, "GBA" = "Gran Buenos Aires")) %>% 
  count(region, nivel_ed, wt = f_calib3) %>% 
  group_by(region) %>% 
  mutate(prop = n / sum (n), 
         prop = scales::percent(prop)) %>% 
  select(-n) %>% 
  ungroup() %>% 
  spread(region, prop, fill = 0) %>% 
  knitr::kable(caption = "Porcentaje de nivel educativo por región para personas de 30 o más años")

La categoría “NS/NR” es residual, su frecuencia es prácticamente 0 en todos los casos. La categoría “Educación especial” es una categoría diferente de nivel educativo diferente a las demás y no podrías considerarla un nivel más alto que Universitario Completo, como implícitamente está en la base. Para ajustar un modelo ordinal eliminamos ambas categorías.

Más allá de estas cuestiones de tipo más técnico, el patrón no es completamente claro. La región “GBA” tiene el porcentaje más alto de educación universitaria completa, pero un porcentaje relativamente bajo de universitaria incompleta. Patagonia la más alta de Terciario completa, pero niveles de población Universitaria completa similares a los del NEA. Parte del problema tiene que ver con el gran número de categorías educativas y con la presencia de grados incompletos. Si bien son adecuados en la escala ordinal y contar con más información es mejor que contar con menos información en este caso complican la comprensión. Tomaré una decisión importante: recodificar los niveles educativos compactando las categorías incompletas a la inmediatamente anterior. Este decisión se funda en la noción de que el logro del título tiene implicaciones sociológicas que ván más allá de la educación que se obtuvo en el proceso y que opera como un capital educativo objetivado.⁵

Recodificación de las variables Nivel educativo e Idioma extranjero

Para recodificar la variable nivel_ed utilizaremos la función fct_collapse, de la librería forcats::. Esta función tiene como primer argumento un objeto de la clase factor (en este caso la columna nivel_ed del data frame) y luego, separados por comas, una definición de los nuevas categorías y un vector definido con c() de las categorías originales que se colapsarán a la nueva. Es posible definir estas relaciones en la llamada a la función, pero el código se hace complicado. Es mejor práctica especificar esta definición en una lista y luego “abrirla” como código a evaluar en la llamada a la función. Para eso utilizaremos el operador !!!⁶ en la llamada a fct_collapse.

A esa lista la llamaremos reco_niveled y tiene una estructura simple: cada elemento de la lista es un vector de la clase character con un nombre dentro de la lista. El nombre del vector será la nueva categoría y los elementos del vector las viejas categorías que la compondrán. Es fundamental observar las convenciones sintácticas de R y ubicar correctamente las comas y las comillas. Es buena práctica insertar saltos de línea en cada definición e identar el código para facilitar la legibilidad.

reco_niveled <- list(
  ninguno = c("Sin instrucción (incluye nunca asistió o sólo asistió a sala de 5)", 
              "Primaria/EGB incompleto"),
  primaria = c("Primaria/EGB completo", "Secundario/Polimodal incompleto"), 
  secundaria = c("Secundario/Polimodal completo", "Terciario incompleto"), 
  terciario = c("Terciario completo", "Universitario incompleto"), 
  universitario = "Universitario completo", 
  otro = c("Educación especial", "NS/NR"))
reco_niveled

## $ninguno
## [1] "Sin instrucción (incluye nunca asistió o sólo asistió a sala de 5)"
## [2] "Primaria/EGB incompleto"                                           
## 
## $primaria
## [1] "Primaria/EGB completo"           "Secundario/Polimodal incompleto"
## 
## $secundaria
## [1] "Secundario/Polimodal completo" "Terciario incompleto"         
## 
## $terciario
## [1] "Terciario completo"       "Universitario incompleto"
## 
## $universitario
## [1] "Universitario completo"
## 
## $otro
## [1] "Educación especial" "NS/NR"

Una de las ventajas de definir a las relaciones entre categorías como una estructura de datos en la sesión de R es que la podemos manejar como cualquier otro objeto. Por ejemplo, crear una tabla bien formateada que relaciona categorías nuevas y viejas. Como la tabla se genera a partir de la misma estructura de datos no es necesario actualizar la tabla si hacemos alguna modificación a la definición de recodificación. Así es R, lo simple es tedioso y lo tedioso simple.

reco_niveled %>% 
  map_df(~`length<-`(.x, 2)) %>% #Porque los largos son desiguales y falla la conversión a df
  gather("Código nuevo", "Códigos originales") %>% 
  drop_na() %>%                  #Porque antes los alargué a 2
  knitr::kable(caption = "Recodificación de la variable Nivel Educativo")

Niveles educativos por región (recodificado)

personas %>% 
  select(v118a, region, f_calib3, t_aglo,  v108,  v109, nivel_ed, aglo) %>%
  filter(v108 >= 30) %>%
# En algún momento tendré que tomar una decisión más responsable sobre el manejo de missing. 
  drop_na %>%        
  mutate_at(vars(-v108, - f_calib3), as_factor) %>% 
# Aplicar la lista de recodificaciones: lista complicada, código simple
  mutate(nivel_ed = fct_collapse(nivel_ed, !!!reco_niveled)) %>% 
#Eliminar el texto entre paréntesis de las etiquetas de región
  mutate(region = fct_relabel(region, str_remove, " \\(([^\\)]+)\\)")) %>% 
  filter(nivel_ed != "otro", 
         v109 != "Otro") %>% 
  mutate(v118a = ordered(v118a)) %>% 
  droplevels() -> activo2

Confusoras

En el análisis bivariado podemos rechazar -informalmente- la hipótesis de independencia estadística entre los tres pares de variables. Hay una relación entre el idioma extranjero y la región en la que se vive, pero a su vez esta variable está relacionada con el nivel educativo alcanzado. Como el nivel educativo no se distribuye de manera homogénea en las regiones cabe la posibilidad de que sea esa distribución desigual del nivel educativo la que nos haga ver, quizás erróneamente, que hay una relación entre idioma extranjero y región, cuando en realidad la relación se da entre nivel educativo e idioma extranjero.

Sería fantástico poder cuantificar la relación entre región y conocimiento de idioma extranjero en unas condiciones ideales en las que el nivel educativo está fijo. Sería fantástico, pero no lo es. Es una posibilidad muy real, es exactamente lo que podemos hacer con un modelo lineal multivariado.

Modelos logit ordinales

El Modelo 1 (en el script mordi1) es un modelo de razones de probabilidad proporcionales con función de enlace canónica, mejor conocido como logit ordinal. La variable dependiente será el nivel educativo medido ordinalmente con las categorías y el orden definidos en la Tabla 2. Aunque la variable define categorías el modelo se ajusta con respecto a una variable latente continua⁷.

Las variables independientes o predictores serán la región (categoría de referencia: Gran Buenos) y nivel educativo (categoría de referencia: Ninguna o Primaria incompleta).

Tanto los coeficientes de cada categoría de los predictores como la propia variable dependiente se calculan como logaritmos de razón de probabilidad, un tipo de magnitud bastante complicado de entender.⁸ Para mantener las cosas simples analizaremos tres características de los coeficientes del modelo:

1. El signo. Coeficientes con signo positivo indican que el predictor⁹ está asociado a un incremento del conocimiento de un idioma extranjero y coeficientes de signo negativo a menor nivel educativo. Esto nos ahorro el problema de interpretar al modelo como lo haríamos con un logit multinomial propiamente dicho, como un sistema de diferencias con respecto a una categoría de referencia.

2. La magnitud. Coeficientes más grandes indican mayores diferencias.

3. Relacionado con lo anterior, la cercanía de un coeficiente a 0 proporcional a su error estándar. Si el coeficiente es pequeño o el error es muy grande no podríamos rechazar la hipótesis de que el coeficiente poblacional o paramétrico es cero. Consideraremos entonces que ese coeficiente es efectivamente 0 y que no hay asociación entre el predictor y el resultado, es decir, que hay independencia.

Esto parece ser innecesariamente complicado y tal vez lo sea.¹⁰ La ventaja es que esos coeficientes se calculan en un nivel fijo de las demás variables. Si toda la asociación entre conocimiento de idioma extranjero y región se explicara por las diferencias educativas entre regiones podríamos esperar que los coeficientes de región se redujeran o su error se incrementara, de tal modo que los interpretaríamos como si fueran 0. En este caso aceptaríamos la hipótesis del nivel educativo como confusora.

Exploración del modelo 3

MASS::polr(v118a~region+t_aglo+nivel_ed+v108+v109, 
           data = activo2,  
           Hess = TRUE) -> mordi3
coefplot(mordi3, 
         innerCI = 2,
         coefficients = names(mordi3$coefficients)) +
  theme_minimal() + 
  labs(x = "Coeficiente estimado + error", 
       y = "Predictor", 
       title = "Modelo 3",
       caption = "Intervalo de confianza a 2 desviaciones estándar") + 
  scale_colour_manual(values = "black")

Con el ajuste por el tamaño del aglomerado la región pampena deja de tener una diferencia significativa con Gran Buenos Aires. En el caso de NEA el resultado es similar, aunque tiene una pendiente negativa importante el error de estimación es suficientemente grande como para que consisderemos que bien podría ser 0. El sexo y la edad no hay nada destacable para reportar: el conocimiento de idioma extranjero es linealmente independiente de estas variables. El tamaño del aglomerado sí es importante y a mayor tamaño, mayor nivel de idiomas extranjeros.