Técnicas de reducción de dimensionalidad robustas y clustering

RWeekend Segunda Edición, Sociedad Ecuatoriana de Estadística

Introducción al análisis no supervisado

El aprendizaje no supervisado es un método de aprendizaje no existe un conocimiento apriori de un fenómeno en los datos. No tenemos una variable dependiente o output. Buscamos entonces encontrar patrones sobre los datos. Las principales técnicas existentes de análisis no supervisado recaen sobre dos categorías:

• Técnicas de reducción de dimensionalidad

• Técnicas de agrupamiento o clustering

A lo largo de este taller visitaremos 3 técnicas, 2 de reducción de dimensionalidad y 1 de agrupamiento. Específicamente, traemos un problema entre manos: Con mira a coadyuvar al diseño de política pública ¿cómo agrupamos provincias, en función a algunas de sus características socioeconómicas, niveles de contagio por COVID-19 y número de detenciones durante el toque de queda?. Para resolverlo compararemos los resultados de usar ACP y ACP Robusto sobre nuestro conjunto de datos para luego agruparlos con k-medias.

Técnicas de reducción de dimensionalidad

Los métodos de reducción de dimensionalidad consisten en resumir y visualizar la información más importante contenida en un dataset. Existen varias técnicas alrededor de esta temática, tanto si asumimos que existen patrones lineales como no lineales en los datos. A continuación se enumeran algunas de estas técnicas, clasificándolas acorde al tipo de datos con el que estemos trabajando y algunos de los paquetes existentes para su uso:

TipoVariables	Tecnica	Librerias
Numéricos	Análisis de componentes principales, t-SNE	base, FactoMineR
Categóricos	Análisis de correspondencias múltiples	FactoMiner, ade4, epMCA
Mixtos	Análisis factorial mixto	FactoMiner

Técnicas de agrupación o clustering

Las técnicas de agrupamiento o clustering nos permiten obtener conocimiento a partir del descubrimiento de patrones existentes en los datos. Específicamente, el objetivo de los métodos de clustering yacen en la identificación de grupos de objetos similares en un conjunto de datos de interés a través de una medida de similaridad entre puntos (e.g la euclideana). A continuación se enumeran algunas de estas técnicas, clasificándolas acorde al tipo de datos con el que estemos trabajando y algunos de los paquetes existentes para su uso:

TipoVariables	Tecnica	Librerias
Numéricos	k-medias, GMM, CLARA	cluster, FactoMineR, mclust
Categóricos	k-modas, otras medidas de distancia	klaR
Mixtos	k-prototypes, otras medidas de distancia	clustMixType

Cabe notar que algunas técnicas podrían ser útiles para varios tipos de datos, dependiendo de la métrica que sea utilizada como distancia.

Resolución del problema

Para iniciar, carguemos las librerías necesarias:

library(tidyverse) # Manejo de datos
library(sp) # Datos geográficos
library(FactoMineR) # Técnicas de red
library(factoextra) # Complemento de visualización
library(rpca) # ACP Robusto
library(highcharter) # Gráficos interactivos
library(stargazer) # Tablas de resultados

Luego cargaremos nuestros datos:

# Carga de variables socioeconómicas y de contagio a nivel de provincia
CoronavirusProvincia_03 = readRDS("Data/CoronavirusProvincia_03.RDS")
# Carga de datos geográficos
provincias = readRDS("Mapas/ec_provincias.RDS")

Hechemos un rápido vistazo a los datos:

CoronavirusProvincia_03 %>% colnames()

 [1] "CODIGO_PROVINCIA"                      "CONFIRMADO_COVID"                     
 [3] "TOQUE_QUEDA_ECU911"                    "AGLOMERACIONES_ECU911"                
 [5] "DEFUNCIONES_RC"                        "DETENIDOS"                            
 [7] "DETENIDOS_SECUNDARIA"                  "DETENIDOS_HOMBRES"                    
 [9] "DETENIDOS_PAREJA"                      "DETENIDOS_SOLTERO"                    
[11] "DETENIDOS_CUANTIL_1_24A"               "DETENIDOS_CUANTIL_2_30A"              
[13] "DETENIDOS_CUANTIL_3_33A"               "DETENIDOS_CUANTIL_4_69A"              
[15] "PROVINCIA"                             "ANOS_ESCOLARIDAD_IND"                 
[17] "TASA_ANALFABETISMO_IND"                "NET_PRIMARIA_IND"                     
[19] "NET_BASICA_IND"                        "NET_SECUNDARIA_IND"                   
[21] "NET_BACHILLERATO_IND"                  "GROSS_PRIMARIA_IND"                   
[23] "GROSS_BASICA_IND"                      "GROSS_SECUNDARIA_IND"                 
[25] "GROSS_BACHILLERATO_IND"                "EMPLEO_BRUTO_IND"                     
[27] "EMPLEO_GLOBAL_IND"                     "EMPLEO_PLENO_IND"                     
[29] "SUBEMPLEO_IND"                         "EMPLEO_NO_REMUNERADO_IND"             
[31] "OTRO_EMPLEO_NO_PLENO_IND"              "DESEMPLEO_IND"                        
[33] "PARTICIPACION_GLOBAL_IND"              "PARTICIPACION_BRUTA_IND"              
[35] "SECTOR_INFORMAL_IND"                   "POBREZA_INGRESOS_IND"                 
[37] "POBREZA_EXTREMA_INGRESOS_IND"          "GINI_IND"                             
[39] "IPM_IND"                               "POBREZA_NBI_IND"                      
[41] "HACINAMIENTO_IND"                      "AGUA_RED_PUBLICA_SENAGUA_IND"         
[43] "AGUA_RED_PUBLICA_IND"                  "EXCRETAS_SENAGUA_IND"                 
[45] "EXCRETAS_IND"                          "DEFICIT_HABITACIONAL_CUALITATIVO_IND" 
[47] "DEFICIT_HABITACIONAL_CUANTITATIVO_IND" "ELECTRICIDAD_IND"                     
[49] "ALUMBRADO_PUBLICO_IND"                 "RECOLECCION_DESECHOS_SOLIDOS_IND"     
[51] "SERVICIOS_BASICOS_IND"                 "POBLACION_2020"                       
[53] "AREA_KM2"                             

Preprocesemos los datos a ser utilizados:

CoronavirusProvincia_03 = CoronavirusProvincia_03 %>% 
  select(-c(DETENIDOS_SECUNDARIA, DETENIDOS_HOMBRES, DETENIDOS_PAREJA, DETENIDOS_SOLTERO,
            DETENIDOS_CUANTIL_1_24A, DETENIDOS_CUANTIL_2_30A,
            DETENIDOS_CUANTIL_3_33A, DETENIDOS_CUANTIL_4_69A))%>% 
  left_join(data.frame(CODIGO_PROVINCIA=provincias$CC_1, 
                       LONG=coordinates(provincias)[,1], 
                       LAT=coordinates(provincias)[,2]))

Joining, by = "CODIGO_PROVINCIA"

MapaProvincias03 = provincias

Creamos algunas variables que permitan realizar comparaciones entre provincias:

CoronavirusProvincia_03 = CoronavirusProvincia_03 %>%
  mutate(RATIO_CONFIRMADOS = CONFIRMADO_COVID/(POBLACION_2020/AREA_KM2),
         RATIO_DEFUNCIONES = DEFUNCIONES_RC/(POBLACION_2020/AREA_KM2),
         RATIO_DETENIDOS = DETENIDOS/(POBLACION_2020/AREA_KM2),
         AGLOMERACIONES_PC = AGLOMERACIONES_ECU911/(POBLACION_2020/AREA_KM2),
         TOQUE_QUEDA_PC = TOQUE_QUEDA_ECU911/(POBLACION_2020/AREA_KM2))

Realizamos una selección de variables para realizar la reducción de dimensionalidad:

CoronavirusProvincia_03 = CoronavirusProvincia_03 %>% 
  select(CODIGO_PROVINCIA, PROVINCIA,
         RATIO_CONFIRMADOS, TOQUE_QUEDA_PC, AGLOMERACIONES_PC, RATIO_DEFUNCIONES, RATIO_DETENIDOS,
         CONFIRMADO_COVID, TOQUE_QUEDA_ECU911, AGLOMERACIONES_ECU911, DEFUNCIONES_RC, DETENIDOS,
         DESEMPLEO_IND, SUBEMPLEO_IND, SECTOR_INFORMAL_IND,
         ANOS_ESCOLARIDAD_IND, 
         POBREZA_INGRESOS_IND, POBREZA_EXTREMA_INGRESOS_IND,
         GINI_IND, IPM_IND,
         HACINAMIENTO_IND, SERVICIOS_BASICOS_IND,
         LONG, LAT)

Observamos algunas medidas de tendencia central descriptivas:

stargazer(as.data.frame(CoronavirusProvincia_03), type = "text", summary = T)

==================================================================================
Statistic                    N   Mean   St. Dev.   Min   Pctl(25) Pctl(75)   Max  
----------------------------------------------------------------------------------
RATIO_CONFIRMADOS            23  0.553   1.016    0.000   0.097    0.502    4.839 
TOQUE_QUEDA_PC               23  8.517   10.831   0.000   2.261    9.442   49.072 
AGLOMERACIONES_PC            23  4.458   2.834    0.883   2.486    5.576   12.549 
RATIO_DEFUNCIONES            23  1.242   1.695    0.324   0.575    1.112    8.627 
RATIO_DETENIDOS              23  0.892   0.890    0.056   0.231    1.199    2.824 
CONFIRMADO_COVID             23 82.000  284.390     0       5        25     1,376 
TOQUE_QUEDA_ECU911           23 517.522 570.718     0     121.5     673     1,864 
AGLOMERACIONES_ECU911        23 394.478 513.760    27      73.5    405.5    1,944 
DEFUNCIONES_RC               23 178.913 505.040     5       12       96     2,453 
DETENIDOS                    23 48.217   65.642     3       14       44      307  
DESEMPLEO_IND                23  0.035   0.018    0.013   0.025    0.042    0.086 
SUBEMPLEO_IND                23  0.168   0.046    0.099   0.143    0.181    0.275 
SECTOR_INFORMAL_IND          23  0.529   0.115    0.237   0.461    0.583    0.750 
ANOS_ESCOLARIDAD_IND         23  9.473   0.964    7.174   9.070    10.051  11.979 
POBREZA_INGRESOS_IND         23  0.308   0.116    0.127   0.232    0.396    0.530 
POBREZA_EXTREMA_INGRESOS_IND 23  0.121   0.084    0.035   0.060    0.159    0.327 
GINI_IND                     23  0.485   0.040    0.424   0.449    0.509    0.575 
IPM_IND                      23  0.448   0.141    0.102   0.333    0.537    0.670 
HACINAMIENTO_IND             23  0.094   0.038    0.030   0.067    0.119    0.174 
SERVICIOS_BASICOS_IND        23  0.628   0.144    0.337   0.513    0.760    0.924 
LONG                         23 -78.738  1.083   -80.578 -79.357  -78.217  -76.375
LAT                          23 -1.461   1.427   -4.170   -2.336   -0.471   0.752 
----------------------------------------------------------------------------------

Escalemos las variables de análisis:

CoronavirusProvincia_03_Scaled = CoronavirusProvincia_03 %>% 
  select("CONFIRMADO_COVID", "TOQUE_QUEDA_ECU911", "AGLOMERACIONES_ECU911", 
         "DEFUNCIONES_RC", "DETENIDOS",
         "RATIO_CONFIRMADOS", "TOQUE_QUEDA_PC", "AGLOMERACIONES_PC", 
         "RATIO_DEFUNCIONES", "RATIO_DETENIDOS",
         "ANOS_ESCOLARIDAD_IND", 
         "POBREZA_INGRESOS_IND", "POBREZA_EXTREMA_INGRESOS_IND",
         "GINI_IND", "IPM_IND",
         "HACINAMIENTO_IND") %>% 
  scale()

Utilicemos ACP sobre los datos:

PCA_res3 = PCA(CoronavirusProvincia_03_Scaled %>% 
                   data.frame() %>% 
                   select(RATIO_DETENIDOS,
                          ANOS_ESCOLARIDAD_IND, 
                          POBREZA_INGRESOS_IND) %>% 
                   as.matrix())

Y ACP Robusto:

RPCA_res3 = rpca(CoronavirusProvincia_03_Scaled %>% 
                   data.frame() %>%
                   select(RATIO_DETENIDOS,
                          ANOS_ESCOLARIDAD_IND, 
                          POBREZA_INGRESOS_IND) %>% 
                   as.matrix(), 
                 trace = F)
L_matrix3 = data.frame(RPCA_res3$L)
names(L_matrix3) = paste0("Dim",1:ncol(L_matrix3))
rownames(L_matrix3) = paste0(CoronavirusProvincia_03$CODIGO_PROVINCIA)

Observemos la matriz del ACP:

PCA_res3$ind$coord %>% head()

        Dim.1      Dim.2       Dim.3
1 -1.21558240  0.5066878 -0.63249670
2  0.03386516 -1.4821559  0.17143625
3 -0.40084948 -1.1562322 -0.48386835
4 -0.41619415 -0.7214811  0.18192222
5  0.04916951 -1.3123051 -0.05146349
6  1.64599620 -2.5000927  0.50378592

La matriz L del ACPR:

RPCA_res3$L %>% head()

            [,1]         [,2]        [,3]
[1,] -0.02765675  0.018291207 -0.03245406
[2,] -0.01705788 -0.014471168 -0.00536073
[3,] -0.43594302  0.032869155 -0.36618408
[4,] -0.18762382  0.004805349 -0.15228448
[5,] -0.01705788 -0.014471168 -0.00536073
[6,] -0.01705788 -0.014471168 -0.00536073

Y la matriz S del ACPR:

RPCA_res3$S %>% head()

     RATIO_DETENIDOS ANOS_ESCOLARIDAD_IND POBREZA_INGRESOS_IND
[1,]       0.0000000            0.5764712           -1.2664365
[2,]      -0.9213479           -1.0892828            0.1831394
[3,]      -0.1891067           -0.9964548           -0.2145111
[4,]      -0.5641851           -0.3316678            0.0000000
[5,]      -0.6899042           -1.0584686            0.0388222
[6,]      -0.8035747           -2.3698880            1.5720040

Realicemos ahora el agrupamiento con ambas matrices.

Primero seleccionemos el número de clusters con la métrica del ancho de silueta para la matriz generada por el ACP:

fviz_nbclust(PCA_res3$ind$coord, kmeans)

Y luego para la matriz L:

fviz_nbclust(L_matrix3, kmeans)

Para la matriz del ACP se sugieren entre 2 y 3 grupos, y para la matriz L se sugieren 3 grupos. Por fines de comparación utilizaremos 3 grupos.

Primero realicemos el k-medias sobre la matriz del ACP:

PCA_KM3 = eclust(PCA_res3$ind$coord, "kmeans", 3, nboot = 20, graph = F, seed = 123)
fviz_cluster(PCA_KM3, L_matrix3)

Y a la par sobre la matriz L:

RPCA_KM3 = eclust(L_matrix3, "kmeans", 3, nboot = 20, graph = F, seed = 123)
fviz_cluster(RPCA_KM3, L_matrix3)

Ahora observemos medidas de validación para los clusters.

Primero la varianza explicada con la el algoritmo hecho sobre matriz del ACP:

PCA_KM3$betweenss/PCA_KM3$totss

[1] 0.6475804

Y sobre la matriz L:

RPCA_KM3$betweenss/RPCA_KM3$totss

[1] 0.9356374

¿Y si observamos el ancho de silueta?

Sobre la matriz del ACP:

fviz_silhouette(PCA_KM3)

Y sobre la matriz L:

fviz_silhouette(RPCA_KM3)

Finalmente, ¿qué tal si lo observamos en un mapa?

Para ello, primero coloquemos un nombre a cada cluster acorde a sus medidas de tendencia central:

CoronavirusProvincia_03 = CoronavirusProvincia_03 %>% 
  mutate(ClusterPCA=paste0(PCA_KM3$cluster),
         ClusterRPCA = paste0(RPCA_KM3$cluster))

CoronavirusProvincia_03 %>%
  select("CONFIRMADO_COVID", "TOQUE_QUEDA_ECU911", "AGLOMERACIONES_ECU911", 
         "DEFUNCIONES_RC", "DETENIDOS",
         "RATIO_CONFIRMADOS", "TOQUE_QUEDA_PC", "AGLOMERACIONES_PC", 
         "RATIO_DEFUNCIONES", "RATIO_DETENIDOS",
         "ANOS_ESCOLARIDAD_IND", 
         "POBREZA_INGRESOS_IND", "POBREZA_EXTREMA_INGRESOS_IND",
         "GINI_IND", "IPM_IND",
         "HACINAMIENTO_IND",ClusterPCA) %>% 
  group_by(ClusterPCA) %>% 
  summarise_all(mean) %>% 
  arrange(RATIO_DETENIDOS)

CoronavirusProvincia_03 %>%
  select("CONFIRMADO_COVID", "TOQUE_QUEDA_ECU911", "AGLOMERACIONES_ECU911", 
         "DEFUNCIONES_RC", "DETENIDOS",
         "RATIO_CONFIRMADOS", "TOQUE_QUEDA_PC", "AGLOMERACIONES_PC", 
         "RATIO_DEFUNCIONES", "RATIO_DETENIDOS",
         "ANOS_ESCOLARIDAD_IND", 
         "POBREZA_INGRESOS_IND", "POBREZA_EXTREMA_INGRESOS_IND",
         "GINI_IND", "IPM_IND",
         "HACINAMIENTO_IND",ClusterRPCA) %>% 
  group_by(ClusterRPCA) %>% 
  summarise_all(mean) %>% 
  arrange(RATIO_DETENIDOS)

Conociendo las descripciones por cluster, observemos un mapa:

library(rjson)
library(tidyverse)
library(highcharter)
library(stringi)
library(viridisLite)
library(scales)

rename_map = function(x){
  x$properties$name=stri_trans_toupper(stri_trans_general(x$properties$name, "Latin-ASCII"))
  return(x)
}

ecuador = fromJSON(file = "https://raw.githubusercontent.com/Rusersgroup/mapa_ecuador/master/ec-all.geo.json")
ecuador$features = lapply(ecuador$features,rename_map)

CoronavirusProvincia_03$PROVINCIA = stri_trans_toupper(stri_trans_general(CoronavirusProvincia_03$PROVINCIA, "Latin-ASCII"))
CoronavirusProvincia_03$PROVINCIA[CoronavirusProvincia_03$PROVINCIA=="SANTO DOMINGO"] = "SANTO DOMINGO DE LOS TSACHILAS"
CoronavirusProvincia_03$ClusterPCA = factor(CoronavirusProvincia_03$ClusterPCA)
CoronavirusProvincia_03$ClusterRPCA = factor(CoronavirusProvincia_03$ClusterRPCA)

CoronavirusProvincia_03$RATIO_DETENIDOS = round(CoronavirusProvincia_03$RATIO_DETENIDOS, 2)
CoronavirusProvincia_03$ANOS_ESCOLARIDAD_IND = number(CoronavirusProvincia_03$ANOS_ESCOLARIDAD_IND, accuracy = 1)
CoronavirusProvincia_03$POBREZA_INGRESOS_IND = paste0(round(CoronavirusProvincia_03$POBREZA_INGRESOS_IND*100,1)," %")
CoronavirusProvincia_03$SECTOR_INFORMAL_IND = paste0(round(CoronavirusProvincia_03$SECTOR_INFORMAL_IND*100,1)," %")
CoronavirusProvincia_03 = plyr::rbind.fill(CoronavirusProvincia_03, data.frame(PROVINCIA="GALAPAGOS"))

hc_map_ec =highchart(type = "map") %>% 
  hc_plotOptions(map = list(
    borderColor="white",
    allAreas = FALSE,
    joinBy = c("name", "PROVINCIA"),
    mapData = ecuador,
    dataLabels = list(enabled = TRUE,
                      format = '{point.name}'))) %>% 
  hc_title(text = "<b>COVID-19: Grupos de pobreza, años de escolaridad y número de detenciones por violación del toque de queda, matriz ACP</b>",
           margin = 20, align = "center",
           style = list(color = "black", useHTML = TRUE)) %>%
  hc_add_series(name = "Baja", data = CoronavirusProvincia_03 %>% filter(ClusterPCA=="2"), color = "yellow") %>% 
  hc_add_series(name = "Media", data = CoronavirusProvincia_03 %>% filter(ClusterPCA=="1"), color = "orange") %>%
  hc_add_series(name = "Alta", data = CoronavirusProvincia_03 %>% filter(ClusterPCA=="3"), color = "brown") %>%
  hc_add_series(name = "Sin dato", data = CoronavirusProvincia_03 %>% filter(is.na(ClusterPCA)), color = "gray") %>%
  hc_tooltip(followPointer =  FALSE, useHTML=TRUE,
             pointFormat="{point.name} <br>
                          <b> {point.RATIO_DETENIDOS} </b> detenidos / personas por km2.  <br>
                          <b> {point.ANOS_ESCOLARIDAD_IND} </b> años de escolaridad promedio. <br>
                          <b> {point.POBREZA_INGRESOS_IND} </b> de personas bajo la línea de pobreza. <br>
                          <b> {point.SECTOR_INFORMAL_IND} </b> de personas con empleo informal.") %>%
  hc_add_theme(hc_theme_economist()) %>%
  hc_add_annotation(labels = list(list(point = list(x = 9000, y = -8500, xAxis = 0, yAxis = 0), text = "COLOMBIA"),
                                  list(point = list(x = 7000, y = 0, xAxis = 0, yAxis = 0), text = "PERÚ"),
                                  list(point = list(x = 0, y = -5000, xAxis = 0, yAxis = 0), text = "OCEANO PACÍFICO")))
hc_map_ec

hc_map_ec =highchart(type = "map") %>% 
  hc_plotOptions(map = list(
    borderColor="white",
    allAreas = FALSE,
    joinBy = c("name", "PROVINCIA"),
    mapData = ecuador,
    dataLabels = list(enabled = TRUE,
                      format = '{point.name}'))) %>% 
  hc_title(text = "<b>COVID-19: Grupos de pobreza, años de escolaridad y número de detenciones por violación del toque de queda, matriz L</b>",
           margin = 20, align = "center",
           style = list(color = "black", useHTML = TRUE)) %>%
  hc_add_series(name = "Baja", data = CoronavirusProvincia_03 %>% filter(ClusterRPCA=="2"), color = "yellow") %>% 
  hc_add_series(name = "Media", data = CoronavirusProvincia_03 %>% filter(ClusterRPCA=="1"), color = "orange") %>%
  hc_add_series(name = "Alta", data = CoronavirusProvincia_03 %>% filter(ClusterRPCA=="3"), color = "brown") %>%
  hc_add_series(name = "Sin dato", data = CoronavirusProvincia_03 %>% filter(is.na(ClusterRPCA)), color = "gray") %>%
  hc_tooltip(followPointer =  FALSE, useHTML=TRUE,
             pointFormat="{point.name} <br>
                          <b> {point.RATIO_DETENIDOS} </b> detenidos / personas por km2.  <br>
                          <b> {point.ANOS_ESCOLARIDAD_IND} </b> años de escolaridad promedio. <br>
                          <b> {point.POBREZA_INGRESOS_IND} </b> de personas bajo la línea de pobreza. <br>
                          <b> {point.SECTOR_INFORMAL_IND} </b> de personas con empleo informal.") %>%
  hc_add_theme(hc_theme_economist()) %>%
  hc_add_annotation(labels = list(list(point = list(x = 9000, y = -8500, xAxis = 0, yAxis = 0), text = "COLOMBIA"),
                                  list(point = list(x = 7000, y = 0, xAxis = 0, yAxis = 0), text = "PERÚ"),
                                  list(point = list(x = 0, y = -5000, xAxis = 0, yAxis = 0), text = "OCEANO PACÍFICO")))
hc_map_ec

Conclusiones

• El ACP normal es sensible a datos atípicos mientras que el ACPR separa los datos atípicos en una matriz dispersa.
• El realizar k-medias sobre la matriz robusta lleva a mejores indicadores de validación que hacerlo sobre la matriz del ACP.
• ¿Respecto a los resultados económicos?