Los datos

Utilizaremos el bloque p7 de la Encuesta Nacional de Migración. La base está disponible en línea en varios formatos, vamos a usar la versión de SPSS porque conserva muy bien los metadatos de los reactivos (etiquetas, códigos). Naturalmente R puede leer archivos de SPSS.

library(foreign)
library(tidyverse)
library(knitr)
library(reshape)
library(poLCA)
library(parallel)
library(broom)

migracion <- read.spss("http://www.losmexicanos.unam.mx/migracion/encuesta_nacional/base_datos/Encuesta_Nacional_de_Migracion.sav", 
                       to.data.frame = TRUE)
diccionario <- tibble(nombre = names(migracion), 
                    etiquetas = str_trim(attr(migracion, "variable.labels")))
migracion <- as.tibble(migracion)

Con el script anterior se carga la base de datos en el entorno de R como un objeto de la clase tibble. Además se genera un diccionario con los nombres largos (las preguntas), que podemos unir con muestros datos y usar para presentar resultados.

Ajuste de un modelo básico con dos clases

En esta primera prueba ajustaremos un modelo con dos clases latentes. En el ajuste de LCA es necesario indicar cuantas clases queremos que se estimen. Se trata de una de las decisiones más importantes en este tipo de análisis¹ y más adelante veremos como aprovechar las cantidades de interés que arroja el LCA para informar esa decisión.

Para ajustar el modelo usamos la función poLCA del paquete del mismo nombre. Esta función recibe cuatro argumentos importantes:

• formula, que es la formula que especifica las variables que entran al modelo y su función como observadas o covariadas.
• data, un data frame con columnas de la clase factor.
• nclass, el número de clases que queremos que estime.
• nrep, el número de veces que se repetirá el ajuste. Como los modelos LCA utilizan el criterio de máxima verosimilitud es posible que, dados unos numeros iniciales “malos”, el modelo converja en un óptimo local y los resultados sean deficientes. Es buena práctica correr varias veces el modelo para controlar este posible problema.

migracion %>% 
  dplyr::select(starts_with("p7_")) %>% 
  names() 

##  [1] "p7_1"  "p7_2"  "p7_3"  "p7_4"  "p7_5"  "p7_6"  "p7_7"  "p7_8" 
##  [9] "p7_9"  "p7_10"

Modelo de dos clases

Ajustamos el modelo y asignamos un nombre al resultado. poLCA produce una salida en consola de todas maneras (mala práctica…).

library(poLCA)
f <- cbind(p7_1,p7_2,p7_3,p7_4,p7_5,p7_6,p7_7,p7_8,p7_9,p7_10)~1

migracion %>% 
  dplyr::select(starts_with("p7_")) %>%  
  poLCA(f, ., nclass = 2, nrep = 4) -> LCA_2clases

## Model 1: llik = -10640.35 ... best llik = -10640.35
## Model 2: llik = -10619.9 ... best llik = -10619.9
## Model 3: llik = -10619.9 ... best llik = -10619.9
## Model 4: llik = -10619.9 ... best llik = -10619.9
## Conditional item response (column) probabilities,
##  by outcome variable, for each class (row) 
##  
## $p7_1
##            Pr(1)  Pr(2)  Pr(3)  Pr(4)  Pr(5)
## class 1:  0.9051 0.0389 0.0498 0.0062 0.0000
## class 2:  0.2536 0.4037 0.2805 0.0553 0.0069
## 
## $p7_2
##            Pr(1)  Pr(2)  Pr(3)  Pr(4)  Pr(5)
## class 1:  0.9385 0.0375 0.0191 0.0049 0.0000
## class 2:  0.1941 0.4568 0.2615 0.0773 0.0103
## 
## $p7_3
##            Pr(1)  Pr(2)  Pr(3)  Pr(4)  Pr(5)
## class 1:  0.9475 0.0430 0.0095 0.0000 0.0000
## class 2:  0.2328 0.4801 0.2150 0.0653 0.0069
## 
## $p7_4
##            Pr(1)  Pr(2)  Pr(3)  Pr(4)  Pr(5)
## class 1:  0.7510 0.1365 0.0949 0.0176 0.0000
## class 2:  0.0665 0.3347 0.4990 0.0844 0.0155
## 
## $p7_5
##            Pr(1)  Pr(2)  Pr(3)  Pr(4)  Pr(5)
## class 1:  0.8983 0.0651 0.0319 0.0031 0.0016
## class 2:  0.1320 0.4790 0.3098 0.0672 0.0121
## 
## $p7_6
##            Pr(1)  Pr(2)  Pr(3)  Pr(4)  Pr(5)
## class 1:  0.7008 0.1461 0.1279 0.0236 0.0016
## class 2:  0.0612 0.3451 0.4778 0.1021 0.0137
## 
## $p7_7
##            Pr(1)  Pr(2)  Pr(3)  Pr(4)  Pr(5)
## class 1:  0.9356 0.0490 0.0122 0.0032 0.0000
## class 2:  0.2986 0.4463 0.1984 0.0464 0.0103
## 
## $p7_8
##            Pr(1)  Pr(2)  Pr(3) Pr(4)  Pr(5)
## class 1:  0.9328 0.0587 0.0085 0.000 0.0000
## class 2:  0.1436 0.4978 0.2813 0.067 0.0103
## 
## $p7_9
##            Pr(1)  Pr(2)  Pr(3)  Pr(4)  Pr(5)
## class 1:  0.9524 0.0418 0.0049 0.0010 0.0000
## class 2:  0.1589 0.5227 0.2405 0.0677 0.0103
## 
## $p7_10
##            Pr(1)  Pr(2)  Pr(3)  Pr(4)  Pr(5)
## class 1:  0.7814 0.0996 0.1062 0.0112 0.0016
## class 2:  0.0601 0.3188 0.5128 0.0929 0.0155
## 
## Estimated class population shares 
##  0.5145 0.4855 
##  
## Predicted class memberships (by modal posterior prob.) 
##  0.5113 0.4887 
##  
## ========================================================= 
## Fit for 2 latent classes: 
## ========================================================= 
## number of observations: 1199 
## number of estimated parameters: 81 
## residual degrees of freedom: 1118 
## maximum log-likelihood: -10619.9 
##  
## AIC(2): 21401.81
## BIC(2): 21814.04
## G^2(2): 9525.526 (Likelihood ratio/deviance statistic) 
## X^2(2): 1.301352e+18 (Chi-square goodness of fit) 
## 

El resultado es algo parco y puede parecer decepcionante. En realidad la información más valiosa está en LCA_2clases. De todos modos con la salida en consola ya podemos ir viendo algunas características de la variable latente.

El primer bloque menciona 4 modelos (los que especificamos con nrep =) e indican el logaritmo de verosimilitud. El primero es un poco más alto, luego se estabiliza en el valor -10619.9. Esto es una buena señal –es muy probable que ese valor refleje el óptimo global- y un recordatorio de por qué es buena idea repetir el ajuste. Si nos hubiéramos quedado con el primero modelo hubiéramos analizado un modelo suboptimo.

El segundo bloque es más interes y presenta probabilidades de respuesta a los ítems –categorías de respuesta en el formulario de encuesta- para cada clase latente. Para la pregunta p7_1 definida como “¿Estaría dispuesto o no estaría dispuesto a permitir que en su casa vivieran personas de otra religión?” y la clase 1 la probabilidad de respuesta Sí es baja, mientras que es mucho más alta la probabilidad de responder `“Sí, en parte” o “No”. Mirando sólo esta pregunta encontramos que la clase 2 corresponde al grupo con actitudes de mayor tolerancia, mientras que la Clase 1 corresponde al grupo de menor tolerancia. Ese patrón se repite en cada pregunta.

El tercer dato a tener en cuenta es la proporción estimada de membresía a cada clase. En este caso la clase 1 cubre a un 48.55% de la muestra, mientras que la clase 2 abarca al 51.45% restante.

De todas maneras y a esta altura del análisis quizás el dato más relevante sea uno de los últimos en aparecer, el BIC. Se trata del criterio de inferencia bayesiano y es muy útil para informar la decisión sobre el número de clases óptimo. Dicho valor no es interpretable de manera absoluta, es necesario ajustar más de un modelo y comparar los BIC resultantes.

Comparar resultados con diferentes números de clases

La selección de número óptimo de clases se facilita cuando podemos comprar el ajuste –manteniendo constantes las variables y su codificación- que obtienen modelos con distinto número de clases. En este contexto un menor BIC se interpreta como el resultado de un modelo que obtiene un mejor ajuste sin agregar demasiados grados de libertad.²

En términos prácticos lo que buscamos es una lista con los resultados de 9 modelos, especificados con 1 a 9 clases. Con esa lista podemos explorar los estadísticos de cada modelo (entre ellos el BIC) para decidir cuál es el número óptimo.

Técnicamente podemos resolverlo de diferentes maneras. Una aproximación es especificar manualmente cada modelo, asignarle un nombre y luego reunirlo en una lista. Otra es generar un bucle for y generar iterativamente los modelos. R tiene una alternativa mucho mejor: funciones auxiliares que se encargan de aplicar una función a una lista, con la posibilidad de hacerlo de manera paralela aprovechando los múltiples núcleos de los procesadores actuales. El ajuste de modelos de clases latentes es un proceso que puede hacerse muy lento cuando aumentamos el número de clases, la cantidad de variables, las categorías de las variables y el número de filas. Es el tipo de tareas ideal para paralelizar. Por eso en lugar de las más conocida lapply utilizaremos mclapply(), del paquete parallel, para aplicar poLCA a una lista de clases, ya que paraleliza de manera automática la aplicación de la función.

mclapply aprovecha una característica de los sistemas operativos UNIX³ para abrir múltiples hilos, por lo que no funciona en entornos Windows. Se presenta una alternativa con lapply() para hacer reproducible al análisis en Windows.

migracion %>% 
  dplyr::select(starts_with("p7_")) -> bloque_p7

# Para UNIX
mclapply(1:9, function(x) {poLCA(f, bloque_p7, nclass = x, nrep = 5)}) -> LCA_1a9

# Para Windows: 
# lapply(1:9, function(x) {poLCA(f, bloque_p7, nclass = x, nrep = 5)}) -> LCA_1a9

En términos de R LCA_1a9 es una lista dentro de la que hay… más listas. Exactamente 9, una por cada modelo. Tenemos toda la información, solo falta consultarla.

map_df(LCA_1a9, glance) %>% 
  kable(caption = "Comparación de modelos con 1 a nueve clases")

map_df(LCA_1a9, glance) %>% 
  dplyr::select(-df, -df.residual, -chi.squared) %>% 
  mutate(modelo = as.character(1:nrow(.))) %>%
  gather(cantidad, valor, -modelo) %>% 
  ggplot(aes(x= modelo, y = valor, group = 1)) + 
    geom_line() + 
    geom_point() + 
    facet_wrap(~cantidad, scales = "free") +
    labs(title = "Comparación de modelos para el bloque p7", 
         subtitle = "La pendiente se aplana a partir de modelos con cinco clases") + 
    theme_minimal()

Análisis de un modelo con 5 clases

Una vez seleccionado el modelo podemos analizarlo con mayor profundidad. En primer lugar nos interesa conocer la conformación de cada clase ¿cómo responden a las preguntas los miembros de cada clase?

Por supuesto, podemos revisar la salida en consola, pero un gráfico facilita las cosas. Para consultar las cantidades en el modelo usamos la función tidy del paquete broom, que nos regresa un data frame en el que cada fila es una combinación de pregunta, respuesta, clase y probabilidad. Con esos elementos podemos hacer nuestro gráfico.

Sugerencia los joins son una excelente alternativa cuando estamos trabajando con etiquetas cortas y etiquetas largas. Si tenemos un diccionario de variables como pares de claves y valores (variables y etiquetas) podemos usar los nombres cortos durante el manejo de datos y, justo antes de graficar o generar una tabla, unir los nombres largos y usarlos para el output.

etiquetas_item <- tibble(outcome = 1:length(levels(migracion$p7_1)), 
                         etiqueta_item = levels(migracion$p7_1)
)

LCA_1a9[[5]] %>% 
  tidy %>% 
  left_join(dplyr::rename(diccionario, variable = nombre)) %>% 
  left_join(etiquetas_item) %>% 
  mutate(etiquetas = str_remove(etiquetas, 
                                "7 ¿Estaría dispuesto o no estaría dispuesto a permitir que en su casa vivieran personas "), 
         etiquetas = str_wrap(etiquetas, 30),
         etiqueta_item = factor(etiqueta_item, 
                                levels = c("No", "Sí, en parte", "Sí", "NS", "NC")), 
         class = paste("Clase", class)) %>%
  ggplot(aes(x = etiquetas, fill = etiqueta_item, y = estimate)) + 
  geom_col() + 
  facet_wrap(~class, nrow = 1) + 
  scale_fill_viridis_d() + 
  theme_minimal() +
  theme(legend.position = "bottom") + 
  labs(fill = "Respuesta", 
       x = "¿Estaría dispuesto o no estaría dispuesto a
       permitir que en su casa vivieran personas ", 
       y = "Probabilidad de respuesta por clase") + 
  coord_flip()

Figure 1: Actitudes de tolerancia. Probabilidades de item/respuesta para 5 clases latentes.

Proporciones de membresía

tibble(clase = paste("Clase", 1:5), 
       prop = LCA_1a9[[5]]$P) %>% 
  ggplot(aes(x = clase, y = prop)) + 
    geom_col(fill = "white", color = "black") + 
    scale_y_continuous(labels = scales::percent) + 
    theme_minimal() + 
    labs(x = NULL, 
         y = NULL, 
         caption = "Elaboración propia
         Datos: Encuesta Nacional de Migración", 
         title = "Mayoría de actitudes tolerantes", 
         subtitle = "Excepto en diversidad de género")

Figure 2: Proporciones de pertenencia a clases latentes

La clase 3, aquella presenta mayor apertura, es la clase modal, como poco del 33% de la muestra, seguida de las clases 2, 1 y 5. La suma de proporciones de la clases 2 y 5, en ambos casos poco tolerantes a la diversidad de género y estigmas asociados alcanza un 40%.

El siguiente paso será incorporar covariadas a este análisis para explorar si estas actitudes varían en dependencia del género, edad, educación y ubicación geográfica de los y las respondentes. Eso será en otra entrega.