Profesor:Dr. José Manuel MAGALLANES REYES, Ph.D

Profesor Principal del Departamento de Ciencias Sociales, Sección de Ciencia Política y Gobierno.
Oficina 223 - Edificio CISEPA / ECONOMIA / CCSS
Telefono: (51) 1 - 6262000 anexo 4302
Correo Electrónico: jmagallanes@pucp.edu.pe

Análisis de Conglomerados

1 Presentación

La idea intuitiva de conglomerar es poder organizar los casos (filas) en subconjuntos o grupos, tal que la similitud entre casos justifique que pertenezca a un grupo y no a otro.

En la sesión anterior, vimos como crear un índice numérico que resuma indicadores de un concepto. Utilicemos los datos ahí preprocesados:

Traigamos algunos datos de los paises del mundo para el ejemplo de esta sesión:

rm(list = ls())

link_idhdemo="https://docs.google.com/spreadsheets/d/e/2PACX-1vRB4YNe2KdIrTmQUAMScYuWcA2ig8d5fKgJBQIlRPVKcryiurAY3dz4Dy8-fpa_MjqmPeTeYet1ggDR/pub?gid=1870508685&single=true&output=csv"

idhdemo=read.csv(link_idhdemo)
names(idhdemo)

##  [1] "country"          "hdiRanking"       "hdi"              "hdiLife"         
##  [5] "hdiSchoolExpec"   "hdiMeanEduc"      "hdiGni"           "ideRanking"      
##  [9] "ideRegime"        "ide"              "ideElectoral"     "ideFunctioning"  
## [13] "ideParticipation" "ideCulture"       "ideLiberties"

2 Transformación de datos

Para este ejercicio sólo usaremos los componentes del IDH. La distribución de los componentes del IDH podemos verla en la Figura 2.1.

boxplot(idhdemo[,c(4:7)],horizontal = F,las=2,cex.axis = 0.5)

Figure 2.1: Distribución de los componentes del IDH

Como primera estrategia cambiemos sus rangos. Elijamos un rango del 0 al 1, cuyo resultado se ve en la Figura ??

library(BBmisc)

## 
## Attaching package: 'BBmisc'

## The following object is masked from 'package:base':
## 
##     isFALSE

boxplot(normalize(idhdemo[,c(4:7)],method='range',range=c(0,10)))

Figure 2.2: Distribución de los componentes del IDH con nuevo rango (0-1)

Una segunda estrategia sería tipificarla ¹. El resultado se muestra en la Figura 2.3.

boxplot(normalize(idhdemo[,c(4:7)],method='standardize'))

Figure 2.3: Distribución de los componentes del IDH tipificados

Nos quedaremos con la segunda opción.

idhdemo[,c(4:7)]=normalize(idhdemo[,c(4:7)],method='standardize')

3 Correlación

Veamos correlaciones entre estas variables tipificadas:

cor(idhdemo[,c(4:7)])

##                  hdiLife hdiSchoolExpec hdiMeanEduc    hdiGni
## hdiLife        1.0000000      0.8057425   0.7659252 0.7838335
## hdiSchoolExpec 0.8057425      1.0000000   0.8159101 0.7311884
## hdiMeanEduc    0.7659252      0.8159101   1.0000000 0.7139462
## hdiGni         0.7838335      0.7311884   0.7139462 1.0000000

Si hubiera alguna correlación negativa sería bueno invertir el rango, tal que el menor sea el mayor y viceversa. Esto no sucede aquí, por lo que no se hace ningún ajuste.

4 Preparación de los datos para la clusterización

No podemos usar la columna Pais en la clusterización, pero tampoco debemos perderla, por lo que se recomienda usar esos nombres en lugar del nombre de fila.

dataClus=idhdemo[,c(4:7)]
row.names(dataClus)=idhdemo$country

Ya con los datos en el objeto dataClus, calculemos la matriz de distancias entre los casos (paises):

library(cluster)
g.dist = daisy(dataClus, metric="gower")

Hay diversas maneras de calculas matrices de distancias entre casos. Cuando las variables son todas numéricas, es comun usar la distancia Euclideana). Hay otras técnicas útiles como la Mahattan (revisar este debate). En nuestro caso, usaremos la distancia Gower útil cuando las variables (columnas) están de diversos tipos de escalas.

5 Procesos de clusterización

Hay diversas estrategías de clusterización. Veremos dos de ellas en nuestro curso:

La técnica de Partición
La técnica de Jerarquización
- Jerarquización Aglomerativa
- Jerarquización Divisiva

5.1 Estrategia de Partición

Como su nombre lo indica, la estrategia de partición busca partir los casos en grupos. El algoritmo básico establece puntos que traten de ser el centro de los demás casos, tal que estos se separen. Claro está, que estos centros de atracción van moviéndose conforme los grupos se van formando, hasta que al final se han partido todos los casos.

Hay diversos algoritmos que buscan una implementación de estos principios básicos. El más conocido es el de K-medias, pero para ciencias sociales tiene la desventaja que requiere que todas las variables sean numéricas, no siendo muy adecuado cuando haya presencia de variables sean categóricas. La alternativa a las necesidades en ciencias sociales es la técnica de k-medoides.

5.1.1 Decidir cantidad de clusters:

La Figura 5.1 sirve para determinar la cantidad de clusters a solicitar (usando el estadístico gap).

## para PAM

library(factoextra)
fviz_nbclust(dataClus, pam,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F)

Figure 5.1: Clusters sugeridos para algoritmo PAM.

5.1.2 Clusterizar via PAM:

La técnica de k-medoides se implementa en la función pam. Esta función retorna diversos valores, en este caso crearemos una columna con la etiqueta del cluster. Usemos la sugerencia de la Figura 5.1, y hallamos:

library(kableExtra)
set.seed(123)
res.pam=pam(g.dist,3,cluster.only = F)

#nueva columna
dataClus$pam=res.pam$cluster

# ver

head(dataClus,15)%>%kbl()%>%kable_styling()

	hdiLife	hdiSchoolExpec	hdiMeanEduc	hdiGni	pam
Afghanistan	-1.1720163	-1.1223433	-1.7899045	-0.9036118	1
Albania	0.6655042	0.2825071	0.7113515	-0.2954005	2
Algeria	0.6546035	0.3425723	-0.2580009	-0.4600137	2
Angola	-1.2150350	-0.4816372	-1.0570320	-0.7236515	1
Argentina	0.5293798	1.4330952	0.6694140	0.0403686	2
Armenia	0.1046748	-0.1644630	0.7245624	-0.3434922	2
Australia	1.6888036	2.5007014	1.1453419	1.4396130	3
Austria	1.3148582	0.8062709	1.0036890	1.6560856	3
Azerbaijan	-0.2350715	-0.0368326	0.4873350	-0.2891916	2
Bahrain	0.9571050	0.9035056	0.6390787	0.9582010	3
Bangladesh	0.1475666	-0.3910236	-0.4659696	-0.7233309	2
Belarus	0.1547871	0.5249118	0.9696084	-0.0622427	2
Belgium	1.3528009	2.0137324	1.0395414	1.5905903	3
Benin	-1.4462954	-0.9535138	-1.3922662	-0.8252918	1
Bhutan	0.0757291	-0.1281194	-1.1314933	-0.5273580	1

5.1.3 Evaluando el uso de PAM

Una manera práctica de ver el desempeño del algoritmo es calcular las silhouettes. Para el caso reciente, veamos la Figura 5.2.

fviz_silhouette(res.pam,print.summary = F)

Figure 5.2: Evaluando resultados de PAM

La Figura 5.2 muestra barras, donde cada una es un país (caso). Mientras más alta la barra, la pertenencia a ese cluster es clara. La barra negativa indica un país mal clusterizado. Para este caso, estos serían los mal clusterizados:

silPAM=data.frame(res.pam$silinfo$widths)
silPAM$country=row.names(silPAM)
poorPAM=silPAM[silPAM$sil_width<0,'country']%>%sort()
poorPAM

## [1] "Chile"  "Latvia"

5.1.4 Verificando etiqueta de clusters

Exploremos el promedio de cada cluster:

aggregate(.~ pam, data=dataClus,mean)

##   pam    hdiLife hdiSchoolExpec hdiMeanEduc     hdiGni
## 1   1 -1.0811743     -1.0681695  -1.1822831 -0.8111138
## 2   2  0.1340138      0.1819924   0.3393855 -0.2301917
## 3   3  1.2520905      1.1502463   1.0317146  1.5181170

El número asignado al cluster no tiene significado necesariamente, por lo que recomiendo hacer cálculos para dárselo. En este caso, las etiquetas ascienden al igual que el promedio, por lo que no es necesario recodificar la etiqueta.

Antes de continuar, guardemos la columna de PAM en la data integrada, y eliminemos la de dataClus.

idhdemo$pamIDHpoor=idhdemo$country%in%poorPAM
idhdemo$pamIDH=as.ordered(dataClus$pam)
dataClus$pam=NULL

5.2 Estrategia Jerárquica

La jerarquización busca clusterizar por etapas, hasta que todas las posibilidades de clusterizacion sean visible. Este enfoque tiene dos familias de algoritmos:

Aglomerativos
Divisivos

5.2.1 Estrategia Aglomerativa

En esta estrategia se parte por considerar cada caso (fila) como un cluster, para de ahi ir creando miniclusters hasta que todos los casos sean un solo cluster. El proceso va mostrando qué tanto esfuerzo toma juntar los elementos cluster tras cluster.

5.2.1.1 Decidir linkages

Aunque se tiene la distancia entre elementos, tenemos que decidir como se irá calculando la distancia entre los clusters que se van formando (ya no son casos individuales). Los tres mas simples metodos:

Linkage tipo SINGLE.
Linkage tipo COMPLETE.
Linkage tipo AVERAGE

Otro metodo adicional, y muy eficiente, es el de Ward. Al final, lo que necesitamos saber cual de ellos nos entregará una mejor propuesta de clusters. Usemos este último para nuestro caso.

5.2.1.2 Decidir cantidad de Clusters

La Figura 5.3 sirve para determinar la cantidad de clusters a solicitar (usando el estadístico gap).

## PARA JERARQUICO

fviz_nbclust(dataClus, hcut,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F,hc_func = "agnes")

Figure 5.3: Clusters sugeridos para algoritmo AGNES.

5.2.1.3 Clusterizar vía AGNES

La función hcut es la que usaremos para el método jerarquico, y el algoritmo aglomerativo se emplea usando agnes. El linkage será ward (aquí ward.D):

set.seed(123)
library(factoextra)

res.agnes<- hcut(g.dist, k = 3,hc_func='agnes',hc_method = "ward.D")

dataClus$agnes=res.agnes$cluster

# ver

head(dataClus,15)%>%kbl()%>%kable_styling()

	hdiLife	hdiSchoolExpec	hdiMeanEduc	hdiGni	agnes
Afghanistan	-1.1720163	-1.1223433	-1.7899045	-0.9036118	1
Albania	0.6655042	0.2825071	0.7113515	-0.2954005	2
Algeria	0.6546035	0.3425723	-0.2580009	-0.4600137	2
Angola	-1.2150350	-0.4816372	-1.0570320	-0.7236515	1
Argentina	0.5293798	1.4330952	0.6694140	0.0403686	2
Armenia	0.1046748	-0.1644630	0.7245624	-0.3434922	2
Australia	1.6888036	2.5007014	1.1453419	1.4396130	3
Austria	1.3148582	0.8062709	1.0036890	1.6560856	3
Azerbaijan	-0.2350715	-0.0368326	0.4873350	-0.2891916	2
Bahrain	0.9571050	0.9035056	0.6390787	0.9582010	2
Bangladesh	0.1475666	-0.3910236	-0.4659696	-0.7233309	2
Belarus	0.1547871	0.5249118	0.9696084	-0.0622427	2
Belgium	1.3528009	2.0137324	1.0395414	1.5905903	3
Benin	-1.4462954	-0.9535138	-1.3922662	-0.8252918	1
Bhutan	0.0757291	-0.1281194	-1.1314933	-0.5273580	2

El dendograma de la Figura 5.4 nos muestra el proceso de conglomeración AGNES:

# Visualize
fviz_dend(res.agnes, cex = 0.7, horiz = T,main = "")

Figure 5.4: Dendograma de AGNES

El eje ‘Height’ nos muestra el “costo” de conglomerar: mientras más corta la distancia mayor similitud y la conglomeracion es más rápida.

5.2.1.4 Evaluando el uso de AGNES

La Figura 5.5 nos muestra las silhouettes para AGNES.

fviz_silhouette(res.agnes,print.summary = F)

Figure 5.5: Evaluando resultados de AGNES

Nótese que también se presentan valores mal clusterizados. Los identificados son estos:

silAGNES=data.frame(res.agnes$silinfo$widths)
silAGNES$country=row.names(silAGNES)
poorAGNES=silAGNES[silAGNES$sil_width<0,'country']%>%sort()
poorAGNES

##  [1] "Bahrain"        "Bhutan"         "Cape Verde"     "Czech Republic"
##  [5] "Estonia"        "Greece"         "Kuwait"         "Lithuania"     
##  [9] "Poland"         "Portugal"       "Saudi Arabia"   "Spain"

5.2.1.5 Verificando etiqueta de clusters

Exploremos el promedio de cada cluster:

aggregate(.~ agnes, data=dataClus,mean)

##   agnes    hdiLife hdiSchoolExpec hdiMeanEduc      hdiGni
## 1     1 -1.1025984     -1.0855778  -1.1832237 -0.81636850
## 2     2  0.2356679      0.2867675   0.3801487 -0.08496801
## 3     3  1.3867428      1.2105608   1.1283409  1.76038883

Estas etiquetas no necesitan recodificación tampoco. Guardemos la columna de AGNES en la data integrada, y eliminemosla de dataClus.

idhdemo$agnesIDHpoor=idhdemo$country%in%poorAGNES
idhdemo$agnesIDH=as.ordered(dataClus$agnes)
dataClus$agnes=NULL

5.2.1.6 Comparando

Veamos qué tanto se parece a la clasificación jerarquica a la de partición:

# verificar recodificacion
table(idhdemo$pamIDH,idhdemo$agnesIDH,dnn = c('Particion','Aglomeracion'))

##          Aglomeracion
## Particion  1  2  3
##         1 54  1  0
##         2  0 70  0
##         3  0 11 29

5.2.2 Estrategia Divisiva

Esta estrategia comienza con todos los casos como un gran cluster; para de ahi dividir en clusters más pequeños.

5.2.2.1 Decidir Cantidad de Clusters

La Figura 5.6 sirve para determinar la cantidad de clusters a solicitar (usando el estadístico gap).

## PARA JERARQUICO

fviz_nbclust(dataClus, hcut,diss=g.dist,method = "gap_stat",k.max = 10,verbose = F,hc_func = "diana")

Figure 5.6: Clusters sugeridos para algoritmo DIANA

5.2.2.2 Clusterizar vía DIANA

La función hcut es la que usaremos para el método jerarquico, y el algoritmo divisivo se emplea usando diana. Aquí una muestra del resultado:

set.seed(123)
res.diana <- hcut(g.dist, k = 4,hc_func='diana')
dataClus$diana=res.diana$cluster
# veamos
head(dataClus,15)%>%kbl%>%kable_styling()

	hdiLife	hdiSchoolExpec	hdiMeanEduc	hdiGni	diana
Afghanistan	-1.1720163	-1.1223433	-1.7899045	-0.9036118	1
Albania	0.6655042	0.2825071	0.7113515	-0.2954005	2
Algeria	0.6546035	0.3425723	-0.2580009	-0.4600137	2
Angola	-1.2150350	-0.4816372	-1.0570320	-0.7236515	1
Argentina	0.5293798	1.4330952	0.6694140	0.0403686	2
Armenia	0.1046748	-0.1644630	0.7245624	-0.3434922	2
Australia	1.6888036	2.5007014	1.1453419	1.4396130	3
Austria	1.3148582	0.8062709	1.0036890	1.6560856	3
Azerbaijan	-0.2350715	-0.0368326	0.4873350	-0.2891916	2
Bahrain	0.9571050	0.9035056	0.6390787	0.9582010	3
Bangladesh	0.1475666	-0.3910236	-0.4659696	-0.7233309	4
Belarus	0.1547871	0.5249118	0.9696084	-0.0622427	2
Belgium	1.3528009	2.0137324	1.0395414	1.5905903	3
Benin	-1.4462954	-0.9535138	-1.3922662	-0.8252918	1
Bhutan	0.0757291	-0.1281194	-1.1314933	-0.5273580	4

El dendograma de la Figura 5.7 nos muestra el proceso de conglomeración AGNES:

# Visualize
fviz_dend(res.diana, cex = 0.7, horiz = T, main = "")

Figure 5.7: Dendograma de DIANA

5.2.2.3 Evaluando el uso de DIANA

La Figura 5.8 nos muestra las silhouettes para DIANA.

fviz_silhouette(res.diana,print.summary = F)

Figure 5.8: Evaluando resultados de DIANA

Nótese que también se presentan valores mal clusterizados. Los identificados son estos:

silDIANA=data.frame(res.diana$silinfo$widths)
silDIANA$country=row.names(silDIANA)
poorDIANA=silDIANA[silDIANA$sil_width<0,'country']%>%sort()
poorDIANA

## [1] "Azerbaijan"   "Ecuador"      "Mongolia"     "Turkmenistan"

5.2.2.4 Verificando Etiqueta

Exploremos el promedio de cada cluster:

aggregate(.~ diana, data=dataClus,mean)

##   diana    hdiLife hdiSchoolExpec hdiMeanEduc       hdiGni
## 1     1 -1.1287230     -1.1169561 -1.24824598 -0.823370585
## 2     2  0.3137105      0.4546312  0.57266780 -0.006070452
## 3     3  1.3052423      1.1835769  1.01701723  1.586748700
## 4     4 -0.1965156     -0.2767992 -0.04875179 -0.539435372

Aquí vemos que las etiquetas no muestran un orden. Este sería el orden:

original=aggregate(.~ diana, data=dataClus,mean)
original[order(original$hdiLife),]

##   diana    hdiLife hdiSchoolExpec hdiMeanEduc       hdiGni
## 1     1 -1.1287230     -1.1169561 -1.24824598 -0.823370585
## 4     4 -0.1965156     -0.2767992 -0.04875179 -0.539435372
## 2     2  0.3137105      0.4546312  0.57266780 -0.006070452
## 3     3  1.3052423      1.1835769  1.01701723  1.586748700

Esas posiciones hay que usarlas para recodificar:

dataClus$diana=dplyr::recode(dataClus$diana, `1` = 1, `4`=2,`2`=3,`3`=4)

Guardemos la columna de DIANA en la data integrada, y eliminemosla de dataClus.

idhdemo$dianaIDHpoor=idhdemo$country%in%poorDIANA
idhdemo$dianaIDH=as.ordered(dataClus$diana)
dataClus$diana=NULL

6 Visualización comparativa

Vamos a usar la matriz de distancia para darle a cada país una coordenada, tal que la distancia entre esos paises se refleje en sus posiciones. Eso requiere una técnica que proyecte las dimensiones originales en un plano bidimensional. Para ello usaremos la técnica llamada escalamiento multidimensional. Veams algunas coordenadas.

# k es la cantidad de dimensiones
proyeccion = cmdscale(g.dist, k=2,add = T) 
head(proyeccion$points,20)

##                                [,1]        [,2]
## Afghanistan             0.444400179  0.13113731
## Albania                -0.147175274 -0.13701510
## Algeria                -0.023300723 -0.09169449
## Angola                  0.331282576  0.03588831
## Argentina              -0.236409012 -0.08313538
## Armenia                -0.048973926 -0.16503890
## Australia              -0.529876811  0.16845260
## Austria                -0.429300483  0.11885998
## Azerbaijan             -0.006636829 -0.16267282
## Bahrain                -0.325596910  0.02662569
## Bangladesh              0.138226916 -0.08817999
## Belarus                -0.165689524 -0.13349112
## Belgium                -0.494630520  0.16277289
## Benin                   0.419354114  0.10895568
## Bhutan                  0.172062052 -0.04766684
## Bolivia                 0.067447666 -0.09021494
## Bosnia and Herzegovina -0.094332269 -0.15218399
## Botswana                0.109290126 -0.06847570
## Brazil                 -0.025209702 -0.10957750
## Bulgaria               -0.118824558 -0.14195524

Habiendo calculado la proyeccción, recuperemos las coordenadas del mapa del mundo basado en nuestras dimensiones nuevas:

# data frame prep:
idhdemo$dim1 <- proyeccion$points[,1]
idhdemo$dim2 <- proyeccion$points[,2]

Aquí puedes ver el mapa:

library(ggrepel)
base= ggplot(idhdemo,aes(x=dim1, y=dim2,label=row.names(dataClus))) 
base + geom_text_repel(size=3, max.overlaps = 50,min.segment.length = unit(0, 'lines'))

Coloreemos el mapa anterior segun el cluster al que corresponden.

6.1 Gráfica de PAM

# solo paises mal clusterizados
PAMlabels=ifelse(idhdemo$pamIDHpoor,idhdemo$country,'')

#base
base= ggplot(idhdemo,aes(x=dim1, y=dim2))  +
    scale_color_brewer(type = 'qual',palette ='Dark2'  ) + labs(subtitle = "Se destacan los países mal clusterizados")

pamPlot=base + geom_point(size=3, 
                          aes(color=pamIDH))  + 
        labs(title = "PAM") 
# hacer notorios los paises mal clusterizados
pamPlot + geom_text_repel(size=4,
                          aes(label=PAMlabels),
                          max.overlaps = 50,
                          min.segment.length = unit(0, 'lines'))

Figure 6.1: Conglomerados PAM en Mapa Bidimensonal de países

6.2 Gráfica de AGNES

# solo paises mal clusterizados
AGNESlabels=ifelse(idhdemo$agnesIDHpoor,idhdemo$country,'')

agnesPlot=base + geom_point(size=3, 
                            aes(color=as.factor(agnesIDH))) +
          labs(title = "AGNES") 
# hacer notorios los paises mal clusterizados
agnesPlot + geom_text_repel(size=4,
                            aes(label=AGNESlabels),
                            max.overlaps = 50,
                            min.segment.length = unit(0, 'lines'))

Figure 6.2: Conglomerados AGNES en Mapa Bidimensonal de países

6.3 Gráfica de DIANA

# solo paises mal clusterizados
DIANAlabels=ifelse(idhdemo$dianaIDHpoor,idhdemo$country,'')

dianaPlot=base + geom_point(size=3,
                            aes(color=dianaIDH)) + 
          labs(title = "DIANA")

# hacer notorios los paises mal clusterizados
dianaPlot + geom_text_repel(size=4,
                            aes(label=DIANAlabels), 
                            max.overlaps = 50,
                            min.segment.length = unit(0, 'lines'))

Figure 6.3: Conglomerados DIANA en Mapa Bidimensonal de países

Nota que en estas técnicas (partición y jerarquica) todo elemento termina siendo parte de un cluster.

Recuerda que la tipificación producirá variables con media igual a cero y desviación típica igual a uno.↩︎

Sesión 7

1 Presentación

2 Transformación de datos

3 Correlación

4 Preparación de los datos para la clusterización

5 Procesos de clusterización

5.1 Estrategia de Partición

5.1.1 Decidir cantidad de clusters:

5.1.2 Clusterizar via PAM:

5.1.3 Evaluando el uso de PAM

5.1.4 Verificando etiqueta de clusters

5.2 Estrategia Jerárquica

5.2.1 Estrategia Aglomerativa

5.2.1.1 Decidir linkages

5.2.1.2 Decidir cantidad de Clusters

5.2.1.3 Clusterizar vía AGNES

5.2.1.4 Evaluando el uso de AGNES

5.2.1.5 Verificando etiqueta de clusters

5.2.1.6 Comparando

5.2.2 Estrategia Divisiva

5.2.2.1 Decidir Cantidad de Clusters

5.2.2.2 Clusterizar vía DIANA

5.2.2.3 Evaluando el uso de DIANA

5.2.2.4 Verificando Etiqueta

6 Visualización comparativa

6.1 Gráfica de PAM

6.2 Gráfica de AGNES

6.3 Gráfica de DIANA