Conglomerados
Cargamos la data:
folder="data"
fileName="IDE_limpio.xlsx"
fileToRead=file.path(folder,fileName)
library(readxl)
ideProvincias=read_excel(fileToRead)Método jerárquico
Veo estructura de toda la base de datos:
str(ideProvincias)## Classes 'tbl_df', 'tbl' and 'data.frame': 195 obs. of 12 variables:
## $ regionUbigeo : chr "010000" "010000" "010000" "010000" ...
## $ provinciaUbigeo : chr "010200" "010300" "010400" "010500" ...
## $ regionNombre : chr "AMAZONAS" "AMAZONAS" "AMAZONAS" "AMAZONAS" ...
## $ provinciaNombre : chr "Bagua" "Bongará" "Condorcanqui" "Luya" ...
## $ pob2012 : num 77438 32317 51802 52185 30236 ...
## $ ide2012 : num 0.662 0.632 0.46 0.605 0.631 ...
## $ identificacion2012 : num 94.6 97.5 86.2 96.2 97.3 ...
## $ medicos2012 : num 14.61 9.01 8.56 12.42 14.88 ...
## $ escolaridad2012 : num 79.8 76.4 52.2 74.7 79.4 ...
## $ AguaDesague2012 : num 64.5 54.8 37.7 43.3 46.5 ...
## $ electrificacion2012: num 67.9 72.2 39.5 67.4 67.5 ...
## $ mapProv : chr "0102" "0103" "0104" "0105" ...- Hago un subconjunto que contenga solo las variables numéricas:
ideProvincias_sub=ideProvincias[c(7:11)]
head(ideProvincias_sub)- Hacemos que los nombres de las provincias sean los nombres de las filas:
row.names(ideProvincias_sub)=ideProvincias$provinciaNombre## Warning: Setting row names on a tibble is deprecated.head(ideProvincias_sub)- Ahora hacemos la estandarización usando scale:
ideProvincias_sub.scaled <- scale(ideProvincias_sub)
#resultado
head(ideProvincias_sub.scaled)## identificacion2012 medicos2012 escolaridad2012
## Bagua -1.09436397 0.36776299 -0.05610104
## Bongará 0.06856468 -0.39755431 -0.35513210
## Condorcanqui -4.49992185 -0.45950881 -2.50574699
## Luya -0.44997860 0.06825845 -0.50598013
## Rodríguez de Mendoza 0.01775379 0.40460944 -0.08859598
## Utcubamba -0.86397982 -0.24720043 -0.28901303
## AguaDesague2012 electrificacion2012
## Bagua 0.39995304 -0.267919475
## Bongará -0.06962008 0.003631957
## Condorcanqui -0.90310156 -2.082170112
## Luya -0.62880972 -0.301012951
## Rodríguez de Mendoza -0.47528375 -0.291439587
## Utcubamba -0.18230687 -0.574083951Ahora ya tenemos puntuaciones Z para todas las columnas.
Ahora, calculo las ‘distancias’. Como aquí estamos usando el paquete NbClust, aquí las distancias se calculan automáticamente. [Pero OJO, toma conciencia de este paso]
Aquí identificamos el número óptimo de clusters
library(NbClust)
nb <- NbClust(ideProvincias_sub.scaled, method = "complete") #utilizamos el método complete
## *** : The Hubert index is a graphical method of determining the number of clusters.
## In the plot of Hubert index, we seek a significant knee that corresponds to a
## significant increase of the value of the measure i.e the significant peak in Hubert
## index second differences plot.
## 
## *** : The D index is a graphical method of determining the number of clusters.
## In the plot of D index, we seek a significant knee (the significant peak in Dindex
## second differences plot) that corresponds to a significant increase of the value of
## the measure.
##
## *******************************************************************
## * Among all indices:
## * 6 proposed 2 as the best number of clusters
## * 9 proposed 3 as the best number of clusters
## * 5 proposed 4 as the best number of clusters
## * 2 proposed 5 as the best number of clusters
## * 1 proposed 13 as the best number of clusters
## * 1 proposed 15 as the best number of clusters
##
## ***** Conclusion *****
##
## * According to the majority rule, the best number of clusters is 3
##
##
## *******************************************************************El número sugerido de clusters es:
length(table(nb$Best.partition))## [1] 3- Aplicamos el algoritmo deseado, sabiendo cuántos clusters pedir (en este caso, 3):
library(factoextra)## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 3.4.4## Welcome! Related Books: `Practical Guide To Cluster Analysis in R` at https://goo.gl/13EFCZalgoritmo="hclust" #hierarchical method
cuantosClusters=3
solucionJerarquica1 <- eclust(ideProvincias_sub.scaled,
FUNcluster = algoritmo,
k = cuantosClusters,
method = "complete", #linkage
graph = FALSE)
fviz_dend(solucionJerarquica1, rect = TRUE, show_labels = T) #veamos el dendograma
- Guardamos los resultados de la aglomeración en un data frame:
idehc = as.data.frame(solucionJerarquica1$cluster)
colnames(idehc) = c("hc") #le cambio el nombre para reconocer que estoy haciendo hierarchichal clusterHacemos merge:
provinciasClust=merge(ideProvincias,idehc,
by.x = 'provinciaNombre', #la variable de ideProvincias que me permitirá hacer merge
by.y=0) # 'by.y=0' es para usar los row.names ¡que en este caso son los nombres de las provincias!
head(provinciasClust)- Ahora podemos explorar nuestros grupos
Veamos una tabla de frecuencias para ver cuántas provincias hay en cada cluster:
table(provinciasClust$hc)##
## 1 2 3
## 91 56 48Si quiero ver, por ejemplo, quiénes están en el cluster 3, hago esto:
provinciasClust[provinciasClust$hc=="3",]$provinciaNombre## [1] "Abancay" "Andahuaylas" "Arequipa" "Ascope"
## [5] "Barranca" "Cajamarca" "Callao" "Camaná"
## [9] "Canchis" "Cañete" "Casma" "Castilla"
## [13] "Chachapoyas" "Chepén" "Chiclayo" "Chincha"
## [17] "Cusco" "Huamanga" "Huancayo" "Huaral"
## [21] "Huaraz" "Huarmey" "Huarochirí" "Huaura"
## [25] "Ica" "Ilo" "Islay" "Jorge Basadre"
## [29] "Lima" "Mariscal Nieto" "Nazca" "Pacasmayo"
## [33] "Palpa" "Pisco" "Piura" "Puno"
## [37] "Recuay" "San Martín" "San Román" "Santa"
## [41] "Sullana" "Tacna" "Talara" "Tambopata"
## [45] "Tarata" "Trujillo" "Tumbes" "Yauli"Puedo pedir las medias por cluster para todas las variables componentes del IDE:
aggregate(cbind(identificacion2012,medicos2012,escolaridad2012,AguaDesague2012,electrificacion2012) ~ hc, data=provinciasClust,FUN=mean)IMPORTANTE: Si calculásemos clusters con data de otro año, podríamos comparar si una provincia ha mejorado o decaído de uno a otro año
- Evaluamos qué tan buena ha sido nuestra aglomeración:
Lo importante del paquete NbClust es que nos permite validar qué tan buena ha sido nuestra aglomeración. Para ello, calculo siluetas:
fviz_silhouette(solucionJerarquica1)## cluster size ave.sil.width
## 1 1 91 0.27
## 2 2 56 0.16
## 3 3 48 0.41
Un valor cercano a 1 (UNO) indica que el valor está bien aglomerado. Si está cercano a 0 (CERO) indica que está entre dos conglomerados. Si es NEGATIVO, está en un cluster erróneo (es decir, luego de finalizado el algoritmo, no se pudo ubicar a esos casos en el grupo adecuado).
El objeto solucionJerarquica1 tiene guarda esta información de silueta. Accedamos a esa información:
siluetas <- solucionJerarquica1$silinfo$widthsVeo a quienes están mal agrupados:
siluetas[siluetas$sil_width<0,]OJO: Ya sabemos que si a los resultados que acabamos de ver les aplicamos la función as.data.frame podremos manipularlos como cualquier tabla de datos. [FIJA PARA EL EXAMEN]
IMPORTANTE: Cuando usamos scale para estandarizar, usamos los valores por defecto de esta función: la media y la desviación típica. Podríamos usar las medianas y las desviaciones de la mediana en vez de estas [cuando hay muchos outliers, por ejemplo]. En este último caso, pueden ser muy diferentes
Método k-medias
Aquí nos vamos directamente al paso 6: aplicar el algoritmo
algoritmo="kmeans" #k-medias
cuantosClusters=3
solucionKmeans1 <- eclust(ideProvincias_sub.scaled,
FUNcluster = algoritmo,
k = cuantosClusters, #como lo hicimos previamente
graph = F)Si dibujásemos la solución NO tendríamos un dendograma, sino un mapa de provincias que las separa y junta según la cercanía o lejanía entre sus perfiles en el IDE:
fviz_cluster(solucionKmeans1, geom = "point", ellipse = F)
Guardo los clusters como un data.frame:
ideK=as.data.frame(solucionKmeans1$cluster)
colnames(ideK)=c("km")Hago el merge:
provinciasClust=merge(provinciasClust,ideK,
by.x = 'provinciaNombre',
by.y=0) # 'by.y=0' está usando los row.names
head(provinciasClust)Saquemos una tabla de frecuencias para ver cuántas provincias están en cada cluster de k-medias:
table(provinciasClust$km)##
## 1 2 3
## 46 53 96Si quiero ver las provincias que están en el conglomerado de 1 de k-medias, hago lo siguiente:
provinciasClust[provinciasClust$km=="1",]$provinciaNombre## [1] "Abancay" "Arequipa" "Ascope" "Barranca"
## [5] "Cajamarca" "Callao" "Camaná" "Cañete"
## [9] "Casma" "Castilla" "Chachapoyas" "Chepén"
## [13] "Chiclayo" "Chincha" "Cusco" "Huamanga"
## [17] "Huancayo" "Huánuco" "Huaral" "Huaraz"
## [21] "Huarmey" "Huarochirí" "Huaura" "Ica"
## [25] "Ilo" "Islay" "Jorge Basadre" "Lima"
## [29] "Mariscal Nieto" "Nazca" "Pacasmayo" "Palpa"
## [33] "Piura" "Puno" "San Martín" "San Román"
## [37] "Santa" "Sullana" "Tacna" "Tahuamanú"
## [41] "Talara" "Tambopata" "Tarata" "Trujillo"
## [45] "Tumbes" "Yauli"Puedo pedir las medias por cluster de k-medias para todas las variables componentes del IDE:
aggregate(cbind(identificacion2012,medicos2012,escolaridad2012,AguaDesague2012,electrificacion2012) ~ km, data=provinciasClust,FUN=mean)Veamos la performance de la técnica en este data:
fviz_silhouette(solucionKmeans1)## cluster size ave.sil.width
## 1 1 46 0.42
## 2 2 53 0.17
## 3 3 96 0.29
E identificamos quienes han sido pobremente agrupados:
siluetasKmeans <- solucionKmeans1$silinfo$widths
#quedarse con los negativos:
siluetasKmeans[siluetasKmeans$sil_width<0,]