> library(patchwork)
> library(tidyverse)
> library(tidymodels)
> library(GDAtools)
La clase pasada habíamos trabajado sobre PCA que operaba sobre variables cuantitativas. Ahora, ¿qué pasa si tenemos variables cualitativas? Si bien no podemos usar PCA, sí podemos echar mano a otras técnicas más adecuadas y con una lógica bastante similar.
Seguramente, la aplicación más famosa de la técnica que vamos a trabajar sea el siguiente gráfico:
El famoso mapa de los espacios y capitales sociales que Bourdieu despliega en “La Distinción”. Luego podemos discutir sobre si el gráfico es claro o no, pero lo que nos interesa destacar aquí es que el mismo fue construido con el llamado Análisis de Correspondencias Múltiples (o MCA por sus sigrlas en inglés)
Tenemos, ahora, una matriz de datos con variables cualitivas. Si recordamos lo que trabajamos el cuatrimestre anterior, podemos retomar la idea de que es posible representar una variable cualitativa en una matriz, como un conjunto de columnas “dummy”: tomábamos una categoría base y podíamos representar la variable como k-1 columnas (k es la cantidad de categorías) dicotómicas.
Por ejemplo, esta variable
> ex <- tibble(id = 1:7, confianza = as.factor(c("Media", "Alta", "Baja", "Media",
+ "Baja", "Baja", "Alta")), voto = as.factor(c("Sí", "No", "Sí", "Sí", "Sí",
+ "No", "Sí")), niv_educ = as.factor(c("Prim.", "Sec.", "Sec.", "Terc. o univ.",
+ "Prim.", "Sec.", "Terc. o univ.")), )
podemos representarla como
> ex_dum <- ex %>%
+ recipe(~.) %>%
+ step_dummy(confianza, voto, niv_educ) %>%
+ prep(training = ex) %>%
+ bake(new_data = ex)
>
> ex_dum
## # A tibble: 7 × 6
## id confianza_Baja confianza_Media voto_Sí niv_educ_Sec.
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 1 0 1 1 0
## 2 2 0 0 0 1
## 3 3 1 0 1 1
## 4 4 0 1 1 0
## 5 5 1 0 1 0
## 6 6 1 0 0 1
## 7 7 0 0 1 0
## # ℹ 1 more variable: niv_educ_Terc..o.univ. <dbl>
Actividad ¿Qué pasó en el código anterior?
Se observa, por ejemplo, que la categoría base que se tomó fue
“Alta”. Así, el perfil
confianza_Baja == 0 & confianza_Media == 0
que el
registro 5 tiene confianza “Alta”.
Pero supongamos que queremos retener todos los niveles, es decir, no tomar ninguna categoría como base:
> ex_dummy2 <- ex %>%
+ recipe(~.) %>%
+ step_dummy(confianza, voto, niv_educ, one_hot = TRUE) %>%
+ prep(training = ex) %>%
+ bake(new_data = ex)
>
>
> ex_dummy2
## # A tibble: 7 × 9
## id confianza_Alta confianza_Baja confianza_Media voto_No voto_Sí
## <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 0 0 1 0 1
## 2 2 1 0 0 1 0
## 3 3 0 1 0 0 1
## 4 4 0 0 1 0 1
## 5 5 0 1 0 0 1
## 6 6 0 1 0 1 0
## 7 7 1 0 0 0 1
## # ℹ 3 more variables: niv_educ_Prim. <dbl>, niv_educ_Sec. <dbl>,
## # niv_educ_Terc..o.univ. <dbl>
Vamos a llamar a esta tabla “tabla disyuntinva completa” y va a ser uno de los insumos fundamentales del MCA. En general, vamos a ver que a esta tabla se la suele llamar \(\textbf{Z}\). Siendo la matriz \(\textbf{Z}\) una matriz de ceros y unos, se debe tomar en cuenta este hecho para deducir la forma particular que tienen las nubes de puntos sobre las que se hará el análisis.
Tenemos, entonces, \(n\) individuos y \(q\) variables. Los marginales de \(\textbf{Z}\) son
Así, la suma total de la tabla es
\[\sum_{i}\sum_{j} z_{ij}= nq\]
A partir \(\textbf{Z}\) vamos a construir la llamada “Tabla de Burt” \(\textbf{B}\) de la siguiente forma:
\[\textbf{B} = \textbf{Z}^T \textbf{Z}\]
Más simple: \(\textbf{B}\), que contiene por bloques las tablas de contingencia de cada pareja de variables cualitativas. Por ejemplo, para el caso de las tablas anteriores, entonces la tabla de Burt puede tener el aspecto
> as_tibble(burt(data.frame(ex %>%
+ select(-id))), rownames = "var")
## # A tibble: 8 × 9
## var confianza.Alta confianza.Baja confianza.Media voto.No voto.Sí
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 confianza.Alta 2 0 0 1 1
## 2 confianza.Baja 0 3 0 1 2
## 3 confianza.Media 0 0 2 0 2
## 4 voto.No 1 1 0 2 0
## 5 voto.Sí 1 2 2 0 5
## 6 niv_educ.Prim. 0 1 1 0 2
## 7 niv_educ.Sec. 1 2 0 2 1
## 8 niv_educ.Terc. … 1 0 1 0 2
## # ℹ 3 more variables: niv_educ.Prim. <dbl>, niv_educ.Sec. <dbl>,
## # `niv_educ.Terc. o univ.` <dbl>
Si llamamos \(D = diag(z_{.j})\) a la diagonal de los márgenes columna de \(\textbf{Z}\), entonces, hacer MCA es diagonalizar la matriz \(\textbf{Z}\) o la matriz \(\textbf{B}\)
\[ \textbf{S} = \frac{1}{q} \textbf{Z}^T \textbf{Z} \textbf{D}^{-1} = \frac{1}{q} \textbf{B} \textbf{D}{^-1}\]
1. Matriz de indicadores (matriz \(Z\))
Se construye una matriz binaria:
2. Matriz de Burt (\(B = Z^t Z\)) La matriz de Burt es una matriz cuadrada y simétrica:
3. Matriz de perfiles y centro Se transforma la matriz para analizar las desviaciones respecto a la independencia:
4. Diagonalización / Descomposición en Valores Singulares (SVD) Se aplica un procedimiento de análisis de componentes principales con métrica chi-cuadrado:
Este paso es lo que permite representar los datos en un espacio reducido.
5. Cálculo de coordenadas y visualización - Se calculan las coordenadas de cada categoría en el nuevo espacio. - Se puede representar también a los individuos.
El resultado es un mapa factorial donde:
Basta de abstracciones. Vamos a trabajar con el siguiente problema: queremos a partir de la Encuesta Permanente de Hogares del III trimestre de 2022 poder construir índices que nos hablen de las condiciones contractuales y laborales en que se desarrolla el trabajo de los asalariados ocupados en la semana de referencia.
Si revisamos la documentación y el diseño de registro vemos que hay muchas variables que podrían servir para tal objetivo. En particular, vamos a ver las siguientes:
Así, podemos ver que hay cierta estructura en las variables:
PP07F1
); uso gratuito de
vivienda (PP07F2
); pago en productos (PP07F3
)
y otros (PP07F4
)PP07G1
), aguinaldo (PP07G2
), días
por enfermedad (PP07G3
), obra social (PP07G4
),
descuentos jubilatorios y/o aportes realizados de forma autónoma por le
trabajador (un monotributista del Estado sería una situación típica)
(PP07H
, PP07I
)PP07J
) y la percepción de un recibo de sueldo
(PP07K
)De esta forma, la idea va a ser combinar estas variables en por lo menos un índice que nos permita resumir dicha información.
Ahora sí, vamos a los datos. En primer lugar, cargamos los paquetes.
Además de tidyverse
vamos a usar los siguientes:
eph
: para la carga y etiquetado de los datos de la
EPHFactoMineR
: para la estimación del modelo de MCAfactoextra
: para algunas visualizacionesggrepel
: para mejorar las visualizaciones> library(tidyverse)
> library(eph)
> library(FactoMineR)
> library(factoextra)
> library(ggrepel)
En primer lugar, cargamos los datos y los etiquetamos:
> df <- get_microdata(year = 2021, trimester = 3) %>%
+ organize_labels() %>%
+ organize_caes() %>%
+ organize_cno()
Luego, vamos a realizar algunas operaciones de preprocesamiento
importantes. Transformamos dos variables a formato
character
.
> df <- df %>%
+ mutate(ESTADO = as.character(ESTADO), CAT_OCUP = as.character(CAT_OCUP), PONDERA = as.numeric(PONDERA))
Luego, recodificamos las variables que vamos a usar y las dicotomizamos. Además, las renombramos para poder interpretarlas mejor.
> df <- df %>%
+ mutate(across(PP07F1:PP07I, ~as.factor(if_else(. == 2, "No percibe", "Percibe")))) %>%
+ mutate(across(PP07J:PP07K, ~as.factor(if_else(. == 1, "Si", "No")))) %>%
+ rename(PP07F1_comida = PP07F1, PP07F2_vivienda = PP07F2, PP07F3_producto = PP07F3,
+ PP07F4_otro_especie = PP07F4, PP07G1_vacaciones = PP07G1, PP07G2_aguinaldo = PP07G2,
+ PP07G3_dias_enf = PP07G3, PP07G4_obra_social = PP07G4, PP07H_descuento_jub = PP07H,
+ PP07I_aporta_si = PP07I, PP07J_turno_diurno = PP07J, PP07K_recibo = PP07K)
Filtramos a las personas ocupadas y a las personas asalariadas:
> df <- df %>%
+ filter(ESTADO == "Ocupado" & CAT_OCUP == "Obrero o empleado")
Luego, pasamos a esitmar el modelo. En las primeras líneas
seleccionamos las 13 variables que vamos a trabajar
(PP07F1_comida:PP07K_recibo
) y eliminamos las que solo
codifican las no respuestas en ninguna (PP07F5
PP07G_59
). Ese dataset de 13 columnas será el input para la
función MCA()
. La misma tiene cuatro argumentos importantes
(aunque tiene varios otros que pueden chequear en la documentación):
X
: los datos de inputncp
: la cantidad de componentes a preservar en el
resultado (en este caso, 2)row.w
: un vector de pesos de las observacionesgraph
: un booleano que es FALSE
si no
queremos que se despliegue otro gráficoVamos a estimar, entonces, un MCA con las 13 columnas mencionadas,
sobre los 14831 registros de asalariados ocupados, nos vamos a quedar
con 2 componentes y vamos a pasarle la variable PONDERA
como vector de pesos.
> mca_ict <- df %>%
+ select(PP07F1_comida:PP07K_recibo, -PP07F5, -PP07G_59) %>%
+ MCA(X = ., ncp = 2, row.w = df$PONDERA, graph = FALSE)
Hemos guardado en el objeto mca_ict
los resultados.
Vamos a explorarlos un poco. Se trata de una lista con varios elementos.
Lo primero que podemos ver son los autovalores y la varianza
explicada.
> mca_ict$eig
## eigenvalue percentage of variance cumulative percentage of variance
## dim 1 0.536465883 53.6465883 53.64659
## dim 2 0.149298171 14.9298171 68.57641
## dim 3 0.076400089 7.6400089 76.21641
## dim 4 0.069641983 6.9641983 83.18061
## dim 5 0.065538250 6.5538250 89.73444
## dim 6 0.055833717 5.5833717 95.31781
## dim 7 0.019356036 1.9356036 97.25341
## dim 8 0.012187432 1.2187432 98.47216
## dim 9 0.005609748 0.5609748 99.03313
## dim 10 0.004359425 0.4359425 99.46907
## dim 11 0.003157818 0.3157818 99.78486
## dim 12 0.002151449 0.2151449 100.00000
Podemos ver un restulado gráfico de la siguiente forma:
> fviz_eig(mca_ict)
Se observa que tenemos una dimensión (la primera) que explica más de un 50% de la varianza y la segunda un 15%. Por lo cual, manteniendo solamente los dos primeros componentes estaríamos reteniendo un 68% de la variabilidad contenida.
Pero, además, es importante poder interpretar los resultados. Para ello vamos a genear un gráfico de las categorías de cada variable.
> mca_ict$var$coord %>%
+ as_tibble(rownames = "var") %>%
+ ggplot(aes(x = `Dim 1`, y = `Dim 2`)) + geom_point() + geom_text_repel(aes(label = var),
+ max.overlaps = 20, size = 3) + ylim(-3, 6) + theme_minimal()
Obtendríamos un resultados similar utilizando la función
fviz_mca_ind()
. Si intentamos analizar qué expresa cada uno
de los ejes, podemos ver que a la izquierda del eje X () tenemos
respuestas positivas en derechos laborales (vacaciones, aguinaldo, días
por enfermedad, tenencia de aportes jubilatorios, ya sean realzados por
el empleador o realizados por el asalariado mismo). A la derecha, en
cambio, se observa la ausencia de estos conceptos. De esta forma,
podemos interpretar que el eje X es algo así como un índice de
percepción de derechos laborales: valores negativos implican mayor
percepción de derechos. Lógicamente, este esquema es algo confuso,
quisiéramos un índice que tenga valores altos cuando la percepción de
derechos es alta. Para resolver esto, simplemente, multiplicamos los
valores de la Dim 1
por -1.
> mca_ict$var$coord %>%
+ as_tibble(rownames = "var") %>%
+ ggplot(aes(x = `Dim 1` * -1, y = `Dim 2`)) + geom_point() + geom_text_repel(aes(label = var),
+ max.overlaps = 20, size = 3) + ylim(-3, 6) + labs(x = "Dim 1", y = "Dim 2") +
+ theme_minimal()
A su vez, el eje Y muestra la percepción de algunos beneficios que podríamos llamar informales/en especie: valores altos en el eje Y muestran respuestas positivas en beneficios en especie (casa, comida, etc.). Podríamos interpretar esto diciendo que la percepción de varios de esos beneficios (casa, comida, etc.) es independiente de la percepción de otros derechos.
Un punto importate es entender cuál es la “importancia” de las diferentes variables en cada dimensión. Una primera intuició podemos obtenerla a partir de ver la cercanía de cada etiqueta sobre cada eje. Cuanto más cerca esté cada etiqueta de cada eje, más importnate será la variable, más correlacionada estará con esa dimensión.
Actividad Nombrar las dos variables más importantes de la dimensión 1 a partir del gráfico anterior.
Pero podemos utilizar otra función para generar un gráfico
fviz_contrib()
. Le pedimos que use el objeto
mca_ict
que contiene nuestras estimaciones, que calcule la
contribución de cada variable de la dimensión 1.
> fviz_contrib(mca_ict, "var", axes = 1)
Podemos ver, entonces que las variables más importantes son: - obra social - aguinaldo - días por enfermedad - descuento jubilatorio
es decir, aquellas variables vinculadas a derechos formales.
|Actividad |¿Cuáles son las variables más relevantes de la Dimensión 2?
Como dijimos más arriba, podemos pensar a la dimensión 1 como un indicador que expresa una variable “latente”: el nivel de condiciones de trabajo, particularmente de aquellos beneficios de carácter más formal.
Tenemos, entonces, una variable cuantitativa que expresa la
co-ocurrencia de respuestas positivas (o negativas) vinculadas a
beneficios laborales formales. Podemos, entonces, calcular el valor de
cada individuo en ese indicador: extraemos del objeto
mca_ict
las coordenadas individuales
mca_ict$ind$coord
, los transformamos en una tibble y los
agregamos como columnas al dataset original.
> df <- df %>%
+ bind_cols(mca_ict$ind$coord %>%
+ as_tibble() %>%
+ rename(dim_1 = `Dim 1`, dim_2 = `Dim 2`))
Luego, lo multiplicamos por -1 para que quede en la escala correcta:
> df <- df %>%
+ mutate(dim_1 = dim_1 * -1)
Entonces, podemos trabajar este índice como una variable cuantitativa cualquiera. Podemos calcular diferentes métricas, generar gráficos, etc. A título de ejemplo, podemos pensar en analizar las diferencias entre los diferentes grados de formalidad del trabajo entre diferentes calificaciones. Vamos a calcular la media (ponderada por la varibale de pesos que la EPH usa):
> df %>%
+ mutate(CALIFICACION = factor(case_when(CALIFICACION %in% c("falta informacion",
+ "Ns.Nc", "otro") ~ "Ns/Nc", TRUE ~ CALIFICACION), levels = c("Profesionales",
+ "Técnicos", "Operativos", "No calificados", "Ns/Nc"))) %>%
+ group_by((CALIFICACION)) %>%
+ summarise(mean = weighted.mean(dim_1, PONDERA), w_median = Hmisc::wtd.quantile(dim_1,
+ PONDERA, probs = 0.5), median = median(dim_1), n = sum(PONDERA), nn = n())
## # A tibble: 5 × 6
## `(CALIFICACION)` mean w_median median n nn
## <fct> <dbl> <dbl> <dbl> <dbl> <int>
## 1 Profesionales 0.271 0.513 0.513 875144 1144
## 2 Técnicos 0.320 0.513 0.513 1805848 2775
## 3 Operativos 0.0432 0.513 0.513 4078965 6989
## 4 No calificados -0.482 -1.13 -1.13 2089994 3834
## 5 Ns/Nc 0.189 0.513 0.277 86291 89
Se observa cómo a mayor calificación las ocupaciones de mayor calificaión aparecen más asociadas a mejores condiciones de contratación.