Reducción de dimensionalidad y segmentación de clusters con variables socio-económicas

Resumen ejecutivo y planteamiento del problema

El siguiente trabajo pretende replicar el abordaje metodológico de una ponencia presentada inicialmente en las Terceras Jornadas de Sociología de la Universidad de Buenos Aires del año 2000.

En la misma, se plantea producir información de utilidad para la toma de decisiones a partir de la reducción de variables complejas, y la consiguiente clasificación de unidades geográficas pertenecientes a los departamentos de San Juan, Mendoza, Catamarca y La Rioja, para lograr una aproximación al impacto del desarrollo en las condiciones de vida de la población de estos departamentos que permita conformar una regionalización de áreas subprovinciales.

Objetivo general

Clasificar en clusters los factores generados y analizar la variabilidad de los mismos.

Objetivos específicos

• Reducir 25 indicadores en 6 factores a partir de un PCA

• Segmentar en 4 clusters los 6 factores según impacto en el desarrollo de condiciones de vida

• Realizar análisis de la varianza para dar cuenta de la variabilidad explicada por los cluesters generados

Anexo metodológico.

En lo que refiere a la información a utilizar, nos valdremos de 26 indicadores de diferentes dominios temáticos, sociales, demográficos, y económicos pertenecientes al Censo de Población Hogares y Vivienda del año 1991 y a Estadísticas Vitales del Ministerio de Salud.

En cuanto a las técnicas, utilizaremos el análisis de componentes princiaples para reducir el espacio multidimensional propio de los datos originales. Con este insumo, se clasificarán los departamentos en 4 estratos, vía análisis de clusters, y diferenciados según impacto del desarrollo en las condiciones de vida.

Finalmente, se evalúa el resultado de esta estratificación midiendo la asociación con el porcentaje de hogares con NBI. Para esto, se realizaron análisis de varianza para corroborar que, si bien la provincia de pertenencia de los departamentos está asociada con la incidencia del NBI, las variaciones de este indicador son más explicados por el estrato definido en este ejercicio.

Pre procesamiento de datos

En esta sección, realizamos un repaso sobre las variables disponibles, el nivel de medición de las mismas, la presencia de datos faltantes, y los valores outliers.

library(ggplot2)
library(tidyverse)
library(dplyr)
library(haven)
library(knitr)
library(kableExtra)
library(VIM)
library(naniar)
library(corrplot)
library(psych)
library(semPlot)
library(cluster)
library(factoextra)
library(lsr)

datos = read_sav("basejor3.sav")

variables1 = tribble(~Variable, ~Etiqueta,
  "N_ORDEN",     "Clave identificatoria de las observaciones",
  "DEPART",      "Departamento",
  "DEPTO",       "Código de departamento",
  "PURBANA",     "Porcentaje de Población Urbana",
  "MIGINTPR",    "Porcentaje de migrantes intraprovinciales",
  "MIGLIMIT",    "Porcentaje de migrantes países limítrofes",
  "HOGNBI",      "Hogares con necesidades básicas insatisfechas",
  "VSINAGUA",    "Viviendas sin agua corriente",
  "VSINELEC",    "Viviendas sin electricidad",
  "HSINGAS",     "Porcentaje de hogares sin gas de red o envasado",
  "CONDESAG",    "Porcentaje de viviendas con desague cloacal a red pública",
  "MORTINF",     "Tasa de mortalidad infantil")

variables2 = tribble(~Variable, ~Etiqueta,
  "MORTNEON",    "Tasa de mortalidad neonatal",
  "NACVIV20",    "Porcentaje de nacidos con madres de 20 años",
  "SINCOBSA",     "Porcentaje de población sin cobertura de salud",
  "ESCMEDIO",    "Tasa neta de escolarización del nivel medio",
  "ANALFAB",     "Tasa de analfabetismo",
  "ASISHTPI",    "Porcentaje que asistió hasta primario incompleto",
  "PRPERHOG",    "Promedio de personas por hogar",
  "JEFMUJER",    "Porcentaje de hogares con jefa mujer",
  "VIVPREC",     "Porcentaje de viviendas precarias",
  "IRRTENVI",    "Porcentaje de hogares en situaciones irregulares de tenencia",
  "TASACTTO",    "Tasa de actividad de ambos sexos",
  "TASACTMU",    "Tasa de actividad de mujeres")

variables3 = tribble(~Variable, ~Etiqueta,
  "TASDESTO",    "Tasa de desocupación de ambos sexos",
  "TASDESMU",    "Tasa de desocupación de mujeres",
  "ASALPUB",    "Porcentaje de asalariados del sector público",
  "CTAPROP",     "Porcentaje de cuentapropistas",
  "PRECASAL",    "Porcentaje de asalariados precarios (sin descuento)",
  "FAC1_1",      "Puntaje factorial 1",
  "FAC2_1",      "Puntaje factorial 2",
  "FAC3_1",      "Puntaje factorial 3",
  "FAC4_1",      "Puntaje factorial 4",
  "FAC5_1",      "Puntaje factorial 5",
  "FAC6_1",      "Puntaje factorial 6")

variables4 = tribble(~Variable, ~Etiqueta,
  "PCIA",        "Provincia",
  "D1",          "Acceso a vivienda y educación",
  "D2",          "Dimensión 2",
  "D3",          "Dimensión 3",
  "D4",          "Dimensión 4",
  "D5",          "Dimensión 5",
  "QCL_1",       "Clasificación jerárquica",
  "ESTRAT",      "Estratificación según IDECV",
  "CLU11_1",     "Enlace promedio entre grupos",
  "FILT",        "Filtro aplicado"
)

kbl(variables1, align = "c") |>  
  kable_styling(full_width = FALSE, position = "c", font_size = 9) |>  
  column_spec(1, width = "2cm") |>   
  column_spec(2, width = "7cm") |>
  row_spec(0, bold = T, background = "#87CEEB") |>
  add_header_above(c("Tabla 1 - Variables del dataset"=2), font_size = 12, align = "l")

Tabla 1 - Variables del dataset
Variable	Etiqueta
N_ORDEN	Clave identificatoria de las observaciones
DEPART	Departamento
DEPTO	Código de departamento
PURBANA	Porcentaje de Población Urbana
MIGINTPR	Porcentaje de migrantes intraprovinciales
MIGLIMIT	Porcentaje de migrantes países limítrofes
HOGNBI	Hogares con necesidades básicas insatisfechas
VSINAGUA	Viviendas sin agua corriente
VSINELEC	Viviendas sin electricidad
HSINGAS	Porcentaje de hogares sin gas de red o envasado
CONDESAG	Porcentaje de viviendas con desague cloacal a red pública
MORTINF	Tasa de mortalidad infantil

kbl(variables2, align = "c") |>  
  kable_styling(full_width = FALSE, position = "c", font_size = 9) |>  
  column_spec(1, width = "2cm") |>   
  column_spec(2, width = "7cm") |>
  row_spec(0, bold = T, background = "#87CEEB") |>
  add_header_above(c("Tabla 2 - Variables del dataset"=2), font_size = 12, align = "l")

Tabla 2 - Variables del dataset
Variable	Etiqueta
MORTNEON	Tasa de mortalidad neonatal
NACVIV20	Porcentaje de nacidos con madres de 20 años
SINCOBSA	Porcentaje de población sin cobertura de salud
ESCMEDIO	Tasa neta de escolarización del nivel medio
ANALFAB	Tasa de analfabetismo
ASISHTPI	Porcentaje que asistió hasta primario incompleto
PRPERHOG	Promedio de personas por hogar
JEFMUJER	Porcentaje de hogares con jefa mujer
VIVPREC	Porcentaje de viviendas precarias
IRRTENVI	Porcentaje de hogares en situaciones irregulares de tenencia
TASACTTO	Tasa de actividad de ambos sexos
TASACTMU	Tasa de actividad de mujeres

kbl(variables3, align = "c") |>  
  kable_styling(full_width = FALSE, position = "c", font_size = 9) |>  
  column_spec(1, width = "2cm") |>   
  column_spec(2, width = "7cm") |>
  row_spec(0, bold = T, background = "#87CEEB") |>
  add_header_above(c("Tabla 3 - Variables del dataset"=2), font_size = 12, align = "l")

Tabla 3 - Variables del dataset
Variable	Etiqueta
TASDESTO	Tasa de desocupación de ambos sexos
TASDESMU	Tasa de desocupación de mujeres
ASALPUB	Porcentaje de asalariados del sector público
CTAPROP	Porcentaje de cuentapropistas
PRECASAL	Porcentaje de asalariados precarios (sin descuento)
FAC1_1	Puntaje factorial 1
FAC2_1	Puntaje factorial 2
FAC3_1	Puntaje factorial 3
FAC4_1	Puntaje factorial 4
FAC5_1	Puntaje factorial 5
FAC6_1	Puntaje factorial 6

kbl(variables4, align = "c") |>  
  kable_styling(full_width = FALSE, position = "c", font_size = 9) |>  
  column_spec(1, width = "2cm") |>   
  column_spec(2, width = "7cm") |>
  row_spec(0, bold = T, background = "#87CEEB") |>
  add_header_above(c("Tabla 4 - Variables del dataset"=2), font_size = 12, align = "l")

Tabla 4 - Variables del dataset
Variable	Etiqueta
PCIA	Provincia
D1	Acceso a vivienda y educación
D2	Dimensión 2
D3	Dimensión 3
D4	Dimensión 4
D5	Dimensión 5
QCL_1	Clasificación jerárquica
ESTRAT	Estratificación según IDECV
CLU11_1	Enlace promedio entre grupos
FILT	Filtro aplicado

Nivel de medición de variables

De las variables recién mencionadas, solo 4 de ellas son categóricas, a saber, N_ORDEN, DEPART, DEPTO y PCIA, mientras que las restantes son de tipo de numérica.

Valores faltantes

A propósito de los valores missing, encontramos que solo la variable FILT contiene 9 registros en esta condición.

gg_miss_var(datos, show_pct = F) +
  labs(title = "Gráfico 1 - Cantidad de valores faltantes por variable") +
  scale_y_continuous(breaks = seq(0, max(colSums(is.na(datos))), by = 4))

Distribución de frecuencias de variables utilizadas

Tenemos una primera aproximación a la forma de distribución de las variables que utilizaremos para el análisis a partir de histogramas de cada una. Esto, además, nos da un indicio de posibles sesgos, valores atípicos, diferencias de escala, entre otras que pueden condicionar el armado de los clusters, generalmente a partir de la formación de grupos alejados del centro de los datos.

datos_long = pivot_longer(datos, cols = c(4:29), names_to = "Variables", values_to = "Porcentajes") |> mutate(Porcentajes=as.numeric(Porcentajes))

# Genero mis histogramas con escala ajustada para cada variable con el parámetro #"scales"


ggplot(datos_long, aes(x = Porcentajes)) +
  geom_histogram(fill = "#87CEEB", bins = 10) +
  facet_wrap(~Variables, scales = "free", ncol = 4) +
  theme_minimal() +
  labs(title = "Histogramas por variable", x = "", y = "Frecuencia absoluta")

Si bien el valor por defecto de la cantidad de barras en ggplot2 es de 30, en este caso recortamos a 10, perdiendo parte del detalle a riesgo de evitar pequeños sobresaltos entre rangos, que visualmente podrían dar apariencia de outliers cuando no es así. Sin embargo observaremos que este fenómeno se observa en variables como el Promedio de Personas por Hogar PRPERHOG que sucede por valores coincidentes para diversos departamentos.

Por otro lado, observamos variables con fuerte sesgo a la derecha, como es el caso del Porcentaje de viviendas con desagüe cloacal a red pública CONDESAG, el Porcentaje de migrantes de países limítrofes MIGLIMIT, Porcentaje de viviendas sin electricidad, entre otras.

Así también tenemos presencia de valores atípicos en algunas de ellas como es el caso de las variables de tasa de actividad.

Finalmente tenemos algunas distribuciones de más de una moda, mientras que otras de ellas parecen agruparse en torno a un centro definido, como es el caso del Porcentaje de cuentapropistas CTAPROP, la tasa de desocupación para ambos sexos TASDESTO, entre otras de ellas.

Estandarización de variables

Este proceso lo llevamos adelante para poder trabajar con variables con escalas de medición muy disímiles por distinto nivel de medición de la variable, sea ordinal, continua o discreta. Lo que aquí se hace es asignar media cero y desviación estándar de uno al conjunto de las variables de mi dataset.

De más está aclarar que esto no resuelve el problema de la presencia de outliers, o del sesgo, en este caso puede resultar pertinente realizar una transformación logarítmica que suavice estos fenómenos, y no se comprima a los valores mínimos y máximos alrededor de la media.

datos[,c(4:29)] = scale(datos[,c(4:29)])

Análisis de componentes principales

Matriz de correlación

Previo a la ejecución de tal análisis conviene generar una matriz que permita dar cuenta de la existencia de correlación para la posible y posterior conformación de los factores. En la misma, los recuadros con los colores azules más intensos representan correlación positiva y fuerte, mientras que los rojos, la negativa.

Podemos ver algunas variables como PRECASAL y ASALPUBL con correlación fuertemente negativa, y HOGNBI con VSINELEC tienen correlación fuerte y positiva. Probablemente tiendas a agruparse en factores en común.

var = datos[,c(4:29)]

matriz = cor(var)

corrplot(matriz, method = "square" , order = "hclust", tl.cex = 0.7, tl.col="black")

Utilidad del análisis factorial de componentes principales

Mediante tal técnica, no se pretende predecir valores para una variable objetivo sino más bien estimar una serie de valores -cargas factoriales- que permitan identificar que variables resultan pertienentes para identificar los factores que expliquen buena parte de la varianza de las observables. Para esto, se espera que las cargas generadas superen un 0.4

Vale decir, se pretende identificar estructuras en los datos a partir de las correlaciones entre sus variables con los factores generados, y por supuesto a partir de ahí lograr reducir la dimensionalidad original.

Análisis de componentes principales o Análisis factorial?

Si bien, siguiendo los pasos de los autores, utilizaremos el PCA (análisis de componentes principales), vale la pena preguntarse en qué casos resulta más adecuado este enfoque frente al análisis factorial.

El análisis factorial parte del supuesto de que existe varianza común compartida entre las variables, la cual se intenta modelar mediante factores latentes no observables que explican las correlaciones entre ellas -comunalidad-. Para verificar si esta estructura subyacente es plausible, se utilizan indicadores como el de Kaiser Meyer Olkin (KMO).

Por su parte, el análisis de componentes principales no parte de un modelo latente, sino que busca explicar la varianza total de los datos, sin distinguir entre varianza común y específica -no explicitada por los factores-. Su objetivo principal es reducir la dimensionalidad del conjunto de datos manteniendo la mayor cantidad posible de información.

En resumen, ambas técnicas permiten reducir la dimensionalidad, pero se diferencian en el tipo de varianza que intentan explicar: el análisis factorial modela la varianza común, mientras que el PCA trabaja con la varianza total -común y específica-.

Criterio para la elección de la cantidad de factores - Scree plot

Para determinar la cantidad de factores o componentes a extraer, nos apoyamos en el scree plot, un gráfico que muestra la cantidad de varianza explicada por cada uno. En el eje Y se representan los autovalores (eigenvalues), que indican cuánta varianza explica cada factor o componente.

Según la regla de Kaiser, se considera apropiado retener aquellos factores cuyos autovalores sean mayores a 1. A partir de esta regla, y observando el gráfico, el análisis factorial parecería justificar la extracción de solo 4 factores.

En cambio, para el caso del PCA, este mismo criterio sugiere retener 6 componentes principales, en concordancia con la elección realizada por los autores.

scree(matriz)

Resultados del PCA

Según los resultados, los seis primeros componentes explican en conjunto el 83% -casi como el 82,5% señalado por los autores- de la varianza total de las variables observadas. Se observa, no obstante, que la ganancia de varianza explicada decrece a partir del cuarto componente: mientras los primeros tres explican en conjunto el 65%, el cuarto aporta un 7% adicional, y los siguientes menos del 6% cada uno.

Interesa resaltar, que independientemente del porcentaje de la varianza explicada por los factores, estos deben tener sentido teórico en función del objetivo propuesto, con lo cual podría resultar pertinente elegir quedarse con más o menos factores de los sugeridos por el criterio de raíz latente.

pca = principal(matriz, nfactors = 6, rotate = "varimax", n.obs = 71)
pca

Principal Components Analysis
Call: principal(r = matriz, nfactors = 6, rotate = "varimax", n.obs = 71)
Standardized loadings (pattern matrix) based upon correlation matrix
           RC1   RC2   RC5   RC4   RC6   RC3   h2    u2 com
PURBANA  -0.70  0.06 -0.41  0.05 -0.16  0.36 0.81 0.185 2.4
MIGINTPR -0.12 -0.12  0.44 -0.28 -0.41  0.38 0.61 0.392 4.0
MIGLIMIT -0.23  0.63  0.14 -0.07 -0.49 -0.11 0.72 0.276 2.4
HOGNBI    0.85 -0.31  0.21 -0.11  0.06  0.09 0.90 0.103 1.5
VSINAGUA  0.90  0.27 -0.06 -0.06 -0.06 -0.06 0.90 0.099 1.2
VSINELEC  0.85 -0.33  0.23  0.04  0.05  0.24 0.95 0.050 1.6
HSINGAS   0.74 -0.37  0.39 -0.01  0.28  0.19 0.95 0.055 2.6
CONDESAG -0.46  0.06 -0.05  0.07 -0.70  0.12 0.72 0.279 1.9
MORTINF  -0.03 -0.11  0.06  0.85 -0.17  0.20 0.81 0.193 1.2
MORTNEON -0.13 -0.07  0.00  0.90 -0.18  0.12 0.88 0.118 1.2
NACVIV20  0.08 -0.19  0.06 -0.29  0.77  0.13 0.75 0.255 1.5
SINCOBSA  0.47  0.61  0.04  0.10  0.18  0.53 0.92 0.077 3.2
ESCMEDIO -0.89 -0.15 -0.16  0.05 -0.17  0.05 0.88 0.120 1.2
ANALFAB   0.62  0.48  0.16  0.04  0.43  0.00 0.83 0.174 2.9
ASISHTPI  0.67  0.49  0.11 -0.10  0.30  0.14 0.83 0.173 2.5
PRPERHOG  0.05  0.06 -0.20 -0.24  0.61 -0.40 0.63 0.368 2.4
JEFMUJER  0.03 -0.64  0.37  0.28  0.03  0.50 0.87 0.133 3.0
VIVPREC   0.84 -0.08 -0.12  0.07  0.12  0.03 0.75 0.254 1.1
IRRTENVI  0.30  0.41 -0.52 -0.19  0.22 -0.36 0.74 0.258 4.3
TASACTTO  0.17  0.16  0.74  0.53  0.03  0.00 0.88 0.118 2.1
TASACTMU  0.02 -0.27  0.70  0.55 -0.02  0.21 0.91 0.092 2.5
TASDESTO -0.35  0.18 -0.84  0.12 -0.03  0.00 0.87 0.126 1.5
TASDESMU -0.24  0.24 -0.87  0.02  0.06 -0.11 0.89 0.106 1.4
ASALPUBL  0.08 -0.83  0.39 -0.09  0.21 -0.11 0.90 0.100 1.7
CTAPROP   0.11  0.13  0.10  0.34 -0.09  0.75 0.73 0.271 1.6
PRECASAL -0.07  0.86 -0.29 -0.17 -0.01  0.09 0.87 0.131 1.4

                       RC1  RC2  RC5  RC4  RC6  RC3
SS loadings           6.50 3.97 3.89 2.69 2.49 1.96
Proportion Var        0.25 0.15 0.15 0.10 0.10 0.08
Cumulative Var        0.25 0.40 0.55 0.66 0.75 0.83
Proportion Explained  0.30 0.18 0.18 0.13 0.12 0.09
Cumulative Proportion 0.30 0.49 0.67 0.79 0.91 1.00

Mean item complexity =  2.1
Test of the hypothesis that 6 components are sufficient.

The root mean square of the residuals (RMSR) is  0.04 
 with the empirical chi square  79.48  with prob <  1 

Fit based upon off diagonal values = 0.99

Con la finalidad de facilitar la interpretación de los factores, tomaremos aquellas cargas mayores a 0.6¹, que se reflejarán en los colores más oscuros de la paleta, aplicando además, tal como realizan los autores, el método de rotación ortogonal varimax, que maximiza la suma de las varianzas de las cargas requeridas de la matriz de factores. Este resulta útil en casos donde existe una elevada intercorrelación entre las variables originales, y que por ende varios indicadores remiten a la misma dimensión, sin permitir distinguir las particularidades de cada factor. Entonces, este método, al intensificar la variabilidad al interior de cada componente sin alterar la configuración espacial original, permite aproximarse a un apoyo único de cada indicador con un factor en particular.

Analizando e interpretando los 6 factores resultantes, en el primero parece prevalecer la dimensión habitacional

Por ejemplo, en el factor 1, las variables con cargas de fuerte correlación positiva parecen estar más vinculadas a la dimensión habitacional -VSINELEC, HSINGAS, HOGNBI, VIVPREC- y educacional -ANALFAB,ASISHTPI,ESCMEDIO-

El segundo parece asimilarse a lo que es el tercer componente del trabajo, más vinculado a la dimensión de inestabilidad sociolaboral -SINCOBSA, MIGLIMIT,PRECASAL,JEFMUJER,ASALPUB-

El tercer componente está pura y exclusivamente asociado a la informalidad de forma positiva siendo el Porcentaje de cuentapropistas el indicador prevaleciente.

El cuarto componente se asocia a indicadores de mortalidad infantil, específicamente a la Tasa de mortalidad neonatal y a la Tasa de mortalidad infantil.

El quinto tiene correlación fuerte con indicadores ligados a la dimensión laboral. Estos son los indicadores de desocupación y actividad.

Mientras que el sexto puede ligarse a condiciones de riesgo para las infancias, por combinar indicadores sanitarios, de tamaño del hogar, y de nacidos vivos por madres menores a 20 años

cargas = as.data.frame(unclass(pca$loadings))
cargas$Variable = rownames(cargas)
cargas_long = pivot_longer(cargas, cols = c(1:6), names_to = "Factores",
                           values_to = "Cargas")


ggplot(cargas_long, aes(x = Factores, y = Variable, fill = Cargas)) +
  geom_tile(color = "white") +
  scale_fill_gradient2(low = "grey", high = "#87CEEB", mid = "white", midpoint = 0) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Cargas factoriales", fill = "Carga")

Estratificación de Departamentos o Análisis de Cluster.

En este eje se pretende identificar espacios subregionales que presenten condiciones equivalentes, es decir, definir subgrupos de unidades espaciales intraprovinciales caracterizados por una alta homogeneidad interna y heterogeneidad externa respecto a los demás estratos.

Para ello, se inicia con una técnica de agrupación jerárquica exploratoria que permite identificar la posible estructura de los datos y orientar la selección del número óptimo de conglomerados. Este análisis considera que los departamentos pertenecen a cuatro jurisdicciones distintas y que están atravesados de forma diferencial por el Índice de Desarrollo Económico en las Condiciones de Vida (IDECV).

A partir de la cantidad de clusters sugerida por los métodos jerárquicos y complementada con métricas como el método del codo o el coeficiente de silueta, se procede a aplicar el algoritmo no jerárquico de k-means. En este método, cada unidad es asignada al cluster cuyo centroide (el promedio de las variables en ese grupo) se encuentra a menor distancia euclídea. Este procedimiento es iterativo: tras cada reasignación de las unidades, se recalculan los centroides hasta que las asignaciones se estabilizan, logrando así una solución que minimiza la variabilidad intra-cluster.

Criterio para la agrupación de departamentos.

Para realizar la agrupación de los departamentos se recurre a la medición de similitud entre unidades a partir de sus valores en un conjunto de dimensiones comunes. En la mayoría de los métodos de clustering, esta similitud se calcula mediante alguna medida de distancia. En este caso, se utiliza la distancia euclídea, definida como la raíz cuadrada de la suma de las diferencias cuadráticas entre las observaciones en cada una de las dimensiones consideradas. Esta medida resulta especialmente apropiada cuando los datos se encuentran en una escala común.

La base para calcular esta distancia se compone de las puntuaciones factoriales obtenidas a partir de cinco componentes principales² seleccionadas por su capacidad explicativa y por facilitar la interpretación posterior de los grupos.

El uso de las puntuaciones factoriales garantiza que las variables se encuentren en una escala comparable, lo que ya había sido asegurado previamente mediante estandarización, y evita que las diferencias en magnitud entre variables influyan indebidamente en la formación de los grupos.

options(scipen = 999)

## Genero el análisis con 5 componentes, especificando rotación y puntuaciones, que ya vienen por default, dicho sea de paso. Y además, utilizamos el dataframe, y no la matriz de correlaciones. Y visualizamos los valores

d = datos[,c(4:29)]

pca = principal(d, nfactors = 5,
                rotate = "varimax", n.obs = 71,
                scores = TRUE)

## Aprovechamos para visualizar las cargas con 5 componentes

cargas = as.data.frame(unclass(pca$loadings))
cargas$Variable = rownames(cargas)
cargas_long = pivot_longer(cargas, cols = c(1:5), names_to = "Factores",
                           values_to = "Cargas")


ggplot(cargas_long, aes(x = Factores, y = Variable, fill = Cargas)) +
  geom_tile(color = "white") +
  scale_fill_gradient2(low = "grey", high = "#87CEEB", mid = "white", midpoint = 0) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Cargas factoriales", fill = "Carga")

## Obtenemos las puntuaciones

puntajes = data.frame(pca$scores)

datos1 = data.frame(d, puntajes)

Luego de visualizar las cargas y obtener las puntuaciones factoriales³, avanzamos con el método de aglomeración jerárquica, observando cómo los conglomerados se van conformando progresivamente en forma ramificada a través del dendrograma.

En particular, utilizamos el método de Ward, que minimiza la suma de cuadrados dentro de los grupos en cada etapa de fusión, lo que favorece la formación de conglomerados compactos y homogéneos.⁴

Este método resulta sensible a los valores atípicos, pero al mismo tiempo permite reducir su impacto al priorizar la homogeneidad interna de los grupos. En este caso, se visualiza el dendrograma solicitando una partición en 4 conglomerados.

agnes_methods = c("ward", "complete", "average", "single")
sapply(agnes_methods, function(met) {
  ag = agnes(scale(datos1[,c(27:31)]), method = met)
  ag$ac 
})

     ward  complete   average    single 
0.8890428 0.8554156 0.8199828 0.7770709 

aglo = agnes(datos1[,c(27:31)], metric = "euclidean",stand=TRUE, method =  "ward")

paleta = c("#999999", "#66B2FF", "#3399FF", "#0077CC")

fviz_dend(aglo, cex = 0.6, k=4, k_colors = paleta, main = "Cluster aglomerativo - Método de Ward")

Método del codo - Elección óptima del número de cluster

Esta técnica sirve para determinar el número óptimo de cluster (k) en algoritmos como k-means o en clustering jerárquico cuando se desea cortar el dendrograma. La misma consiste en un cálculo para una cantidad k de cluster que se visualizará en el eje X, y que generalmente suele ser la suma total de errores al cuadrado, o el R cuadrado, que mide la proporción de la varianza explicada; valores ambos que se verán en el eje Y.

fviz_nbclust(datos1[,c(27:31)], kmeans, method = "wss") +
  geom_vline(xintercept = 6, linetype = 2)+
  labs(subtitle = "Método del codo")

En el gráfico vemos que la cantidad óptima de cluster puede ir entre 2 y 6 de ellos. El punto de inflexión parece verse en el conglomerado 6. Pasando de 6 a 7 cluster, la ganacia en homogeneidad interna o reducción de la varianza interna, no parece ser sustancial. Aunque también de 3 a 4 la reducción no parece significativa.

Método de Silhoutte

Otra forma de validar la elección del número de clusters es a través del coeficiente de silueta. Este coeficiente evalúa la calidad del agrupamiento considerando dos aspectos:

Cohesión interna: medida como el promedio de la distancia entre cada observación y los demás puntos de su mismo cluster (qué tan bien está asignada al grupo).

Separación externa: medida como el promedio de la distancia entre esa observación y los puntos del cluster más cercano distinto (qué tan mal estaría si se asignara al segundo mejor grupo).

fviz_nbclust(datos1[,c(27:31)], kmeans, method = "silhouette") +
  labs(subtitle = "Método Silhouette")

En este caso, vemos que el número óptimo se encuentra entre 2 y 5 cluster.

El método de aglomeración no jerárquico - K Means.

Una vez identificada la cantidad óptima de conglomerados mediante métodos jerárquicos y criterios como el codo o el coeficiente de silueta, se procede a aplicar un método de agrupamiento no jerárquico, conocido como K-means. Este algoritmo busca minimizar la suma de las distancias cuadradas entre cada observación y el centroide de su cluster, logrando así una buena cohesión al interior de cada subgrupo.

En su funcionamiento, cada observación es asignada al centroide más cercano según una métrica de distancia, que en este caso es la euclídea. Luego, los centroides se recalculan como el promedio de las observaciones asignadas a cada cluster. Este proceso se repite iterativamente hasta que los centroides se estabilizan y la asignación de observaciones no cambia o hasta que se alcanza un número máximo de iteraciones.

set.seed(123)

## Generamos los cluster

d = dist(datos1[,c(27:31)], method = "euclidean")

kmedias = kmeans(d, centers = 4)

## Obtenemos el cluster de pertenencia de cada departamento y agrupamos

clusters = kmedias$cluster

datos2 = data.frame(datos[,c(2,3,36)],datos1, clusters)

agrupacion = datos2 |> group_by(clusters) |>
  summarise(Frecuencia=n()) |> ungroup() |>
  mutate(Porcentaje=round(Frecuencia/sum(Frecuencia)*100,1)) |>
  rename(Estrato=clusters)

kbl(agrupacion, align = "c") |>  
  kable_styling(full_width = FALSE, position = "c", font_size = 9) |>  
  column_spec(1:3, width = "4cm") |>
  row_spec(0, bold = T, background = "#87CEEB") |>
  add_header_above(c("Tabla 5 - Distribución de frecuencias para subgrupos"=3), font_size = 12, align = "l")

Tabla 5 - Distribución de frecuencias para subgrupos
Estrato	Frecuencia	Porcentaje
1	14	19.7
2	33	46.5
3	18	25.4
4	6	8.5

alargado = datos2 |> pivot_longer(cols = c(30:34), names_to = "Componentes", values_to = "Puntuaciones")

agrupacion2 = alargado |> group_by(clusters, Componentes) |>
  summarise(Media=mean(Puntuaciones, na.rm = T)) |> ungroup() |>
  pivot_wider(id_cols = clusters, names_from = Componentes, values_from = Media) |> rename(Estrato=clusters)

otro = datos2 |> group_by(PCIA,clusters) |> summarise(desmuj=mean(TASDESMU), desto=mean(TASDESTO))

colnames(agrupacion2) = c("Estrato","Déficit habitacional y educacional","Inestabilidad sociolaboral","Cuentapropistas","Mortalidad infantil","Participación en mercado laboral")

kbl(agrupacion2, align = "c") |>  
  kable_styling(full_width = FALSE, position = "c", font_size = 9) |>  
  column_spec(1:6, width = "3cm") |>
  row_spec(0, bold = T, background = "#87CEEB") |>
  add_header_above(c("Tabla 6 - Promedio de puntuaciones factoriales por Estrato"=6), font_size = 12, align = "l")

Tabla 6 - Promedio de puntuaciones factoriales por Estrato
Estrato	Déficit habitacional y educacional	Inestabilidad sociolaboral	Cuentapropistas	Mortalidad infantil	Participación en mercado laboral
1	0.3613513	0.7194380	-0.2695134	-0.8279861	-1.0746515
2	-0.0636063	-0.7411117	-0.2427259	-0.1495436	0.5158453
3	-0.4741980	0.7940957	0.2105587	0.4801395	-0.2039036
4	0.9292755	0.0151385	1.3321809	1.3140388	0.2820820

La caracterización de los conglomerados generados según la distribución de la puntuación media de cada componente es la siguiente:

Estrato 1: Son departamentos que tienen una elevada inestabilidad sociolaboral, aunque de baja informalidad. Tienen baja incidencia de condiciones riesgosas para la infancia y una baja accesibilidad al mercado laboral.

Estrato 2: Los departamentos de este estrato no presentan valores deficitarios -casi que neutros- en cuanto acceso a la educación y la vivienda, ni en cuanto a la mortalidad infantil. Además, ostentan un mercado laboral de accesibilidad entre moderada/alta y baja inestabilidad sociolaboral.

Estrato 3: Un cluster complejo. Pese a mejores condiciones habitacionales y educacionales, presenta alta inestabilidad sociolaboral y baja accesibilidad al mercado, con mortalidad infantil entre moderada y alta.

Estrato 4: Estrato extremadamente vulnerable estructuralmente. Condiciones habitacionales críticas y alta mortalidad infantil. Sin embargo, con buena estabilidad sociolaboral formal y moderada participación en el mercado laboral, aunque posiblemente en sectores cuentapropistas.

agrupacion = datos2 |> group_by(clusters,PCIA) |>
  summarise(Frecuencia=n()) |> ungroup() |>
  pivot_wider(id_cols = clusters, names_from = PCIA, values_from = Frecuencia)

colnames(agrupacion) = c("Estrato","Mendoza","San Juan","Catamarca","La Rioja")

agrupacion[is.na(agrupacion)] = 0

kbl(agrupacion, align = "c") |>  
  kable_styling(full_width = FALSE, position = "c", font_size = 9) |>  
  column_spec(1:5, width = "4cm") |>
  row_spec(0, bold = T, background = "#87CEEB") |>
  add_header_above(c("Tabla 7 - Cantidad de departamentos según provincia y estrato de pertenencia"=5), font_size = 12, align = "l")

Tabla 7 - Cantidad de departamentos según provincia y estrato de pertenencia
Estrato	Mendoza	San Juan	Catamarca	La Rioja
1	1	13	0	0
2	0	2	13	18
3	15	3	0	0
4	2	1	3	0

En esta tabla, la distribución de frecuencias nos permite caracterizar las provincias según los estratos ya descritos. En ese caso, los departamentos san juaninos se concentran en el estrato número 1 de baja mortalidad infantil y baja accesibilidad al mercado laboral con elevada informalidad o inestbailidad sociolaboral.

En el estrato número 2, Catamarca y La Rioja presentan las mejores condiciones en términos del ya citado IDECV con baja informalidad, buen acceso al mercado laboral y escasos problemas habitacionales y educacionales.

En el número 3 tienen mayor peso los departamentos de mendoza caracterizados por tener mayores niveles de informalidad y buenos resultados en materia de indicadores educativos y vivienda.

El cuarto estrato de condiciones más vulnerables en el conjunto de los indicadores reflejados en los 5 componentes, se reparte entre 6 departamentos de tres provincias distintas.

A diferencia de la conclusión a la que llegan a los autores, parece haber dependencia entre la provincia de pertenencia y los estratos del IDECV. Es por esta independencia a priori a la que llegan que deciden realizar un último paso con el Análisis de la varianza. Esto es, probar la hipótesis de existencia o no de asociación del porcentaje de hogares con NBI, a su jurisdicción de pertenencia, o a su estrato de pertenencia

Análisis de la varianza.

En este eje final, avanzamos con el test de ANOVA para corroborar si existe diferencia entre las medias del porcentaje de hogares con NBI según estratos, y luego según jurisdicción/provincia.

## ANOVA Y ETA CUADRADO DE CLUSTER GENRADOS

anova_estratos = aov(HOGNBI ~ as.factor(clusters), data = datos2)


summary(anova_estratos)

                    Df Sum Sq Mean Sq F value   Pr(>F)    
as.factor(clusters)  3  18.62   6.208   8.095 0.000112 ***
Residuals           67  51.38   0.767                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

etaSquared(anova_estratos)

                       eta.sq eta.sq.part
as.factor(clusters) 0.2660422   0.2660422

## ANOVA Y ETA CUADRADO DE PCIAS

anova_pcia = aov(HOGNBI ~ as.factor(PCIA), data = datos2)

summary(anova_pcia)

                Df Sum Sq Mean Sq F value    Pr(>F)    
as.factor(PCIA)  3  19.68   6.560   8.734 0.0000572 ***
Residuals       67  50.32   0.751                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

etaSquared(anova_pcia)

                   eta.sq eta.sq.part
as.factor(PCIA) 0.2811369   0.2811369

Según los resultados del eta cuadrado, es la variable de la jurisdicción o provincia la que explica un porcentaje superior de la varianza de la variable dependiente, siendo de 28,1% contra 26,6% de nuestros cluster generados. Esto puede ser posible en tanto los autores no hayan incluido el indicador NBI en el análisis inicial de componentes principales, siendo el resultado de su estadístico, notablemente superior y arrojando 0.87.

Conclusiones finales.

Replicando el abordaje llevado a cabo por los autores, desarrollamos un PCA que permitió reducir 26 variables observables en 6 y luego 5 componentes conformados por los mismos atributos identificados por los autores, con cargas factoriales idénticas, aunque con distintos puntos de corte para la interpretación de los mismos.

Al calcular los puntuaciones con el método de Bartlett, los resultados variaron aunque se mantuvieron bastante idénticos con los que obtuvieron los autores.

Es a partir de la aglomeración no jerárquica que se observan algunas diferencias en la distribución de los atributos predominantes de los factores al interior de cada cluster. El caso del estrato 1, de baja inestabilidad sociolaboral y mortalidad infantil, se caracteriza, además, por ser de baja accesibilidad al mercado de trabajo, siendo este mismo estrato (que es el estrato 2 en el trabajo) de alta accesibilidad. Lo mismo sucede en el estrato 2 (estrato 3 para los autores) con la accesibilidad al mercado laboral.

Finalmente, hemos comprobado la validez del análisis de cluster con que decidimos trabajar. A diferencia del resultado de los autores, nuestro ANOVA no performó tan bien pero sí discriminó correctamente la agrupación de nuestras unidades por provincia a las condiciones de pobreza, más que por nuestros cluster generados.

Restaría, entonces, probar nuevos modelos para la agrupación de variables y observaciones, que a pesar de explicar un porcentaje no tan alto de la variabilidad de una variable en particular, pueden resultar de utilidad para identificar estructuras intrínsecas en los datos que permitan clasificar espacialmente a unidades regionales del territorio argentino que sea de utilidad para la posterior intervención sobre los mismos.

Notas

1. Los autores, de todas formas, trabajan con cargas superiores a 0.5 tal como lo indican en su matriz de factores rotados.

2. Se excluyó una sexta dimensión, identificada como Inserción incompleta en el mercado laboral, debido a su baja carga factorial y a que su varianza explicada se encontraba justo por debajo del umbral establecido. Esta exclusión no compromete la comunalidad, dado que los cinco componentes retenidos explican en conjunto el 78% de la variabilidad total.

3. Una puntuación factorial alta para una unidad indica que la misma tiene una fuerte presencia de las características capturadas por ese factor, o lo que es lo mismo, cuánto representa ese factor a una observación específica.

4. Previo a la elección de tal, calculamos el coeficiente de aglomeración con distintos métodos, y nos quedamos con el más alto.