Análisis de Cluster
LRomero, MRamirez, JRojas, EDarghan
27/6/2020
Introducción
Una de las habilidades que marcan a las criaturas vivas es la capacidad de agrupar elementos similares en categorías que comparten ciertas propiedades. La clasificación hace parte esencial en varias disciplinas como el lenguaje, la biología, la química, la astronomía, etc. La clasificación se desarrolla principalmente como un método para organizar grandes cantidades de datos y describir patrones de similaridad y disimilaridad.
A continuación se presentan una serie de ejercicios que permitiran interiorizar el manejo de R y su aplicación en métodos de clasificación por clusters. El documento presenta una serie de ejercicios con datos generados y disponibles en el paquete Iris, para demostrar el uso principalmente de los métodos K-means y CLARA. Este contenido hace parte del contenido académico desarrollado por el profesor Aquiles Enrique Darghan Contreras.
Previo a la solución de los ejercicios es necesario instalar y cargar las siguientes librerías.
library(ggplot2)
library(car)
library(dplyr)
library(factoextra)
library(cluster)
library(simstudy)
library(data.table)
library(tidyverse)
library(PerformanceAnalytics)
library(corrr)Datos
Inicialmente se construirán una serie de datos ajustados para ejemplificar los métodos de clasificación por Cluster.
x <- c(4, 4.5, 4, 7.5, 7, 6, 5, 5.5, 5, 6)
y <- c(4, 4.5, 4, 7.5, 7, 6, 5, 5.5, 5, 6) + 4
z <- c(5, 5.5, 4.8, 5.4, 4.7, 5.6, 5.3, 5.5, 5.2, 4.8)
datos <- data.frame(punto = rep(c("x", "y", "z"), each = 10),
respuesta = c(x,y,z),predictor = 1:10)
datos## punto respuesta predictor
## 1 x 4.0 1
## 2 x 4.5 2
## 3 x 4.0 3
## 4 x 7.5 4
## 5 x 7.0 5
## 6 x 6.0 6
## 7 x 5.0 7
## 8 x 5.5 8
## 9 x 5.0 9
## 10 x 6.0 10
## 11 y 8.0 1
## 12 y 8.5 2
## 13 y 8.0 3
## 14 y 11.5 4
## 15 y 11.0 5
## 16 y 10.0 6
## 17 y 9.0 7
## 18 y 9.5 8
## 19 y 9.0 9
## 20 y 10.0 10
## 21 z 5.0 1
## 22 z 5.5 2
## 23 z 4.8 3
## 24 z 5.4 4
## 25 z 4.7 5
## 26 z 5.6 6
## 27 z 5.3 7
## 28 z 5.5 8
## 29 z 5.2 9
## 30 z 4.8 10Visualización de datos
ggplot(data = datos, aes(x = as.factor(predictor), y = respuesta,
colour = punto)) + geom_path(aes(group = punto)) +
geom_point() + labs(x = "predictor") + theme_bw() +
theme(legend.position = "bottom")
Una primera aproximación para determinar la similaridad entre obervaciones es calculando la distancia total entre los puntos. Uno de los métodos más empleados es la distancia euclidiana. Hay que tener en cuenta, que entre más distantes sean dos puntos, más disimilares serán o menos correlacionados serán.
x1=matrix(c(x,y))
x2=matrix(c(x,z))
x3=matrix(c(y,z))
sum(dist(x1, method = "euclidean",diag=F,upper=F,p=2),1)## [1] 528sum(dist(x2, method = "euclidean",diag=F,upper=F,p=2),1)## [1] 180.3sum(dist(x3, method = "euclidean",diag=F,upper=F,p=2),1)## [1] 509.3R=matrix(c(x,y,z),ncol=3)
dcor=1-cor(R)
dist(x = rbind(x,y,z), method = "euclidean")## x y
## y 12.64911
## z 3.75100 13.98821dcor <- 1 - cor(x = cbind(x,y,z), method = "pearson")
dcor## x y z
## x 0.0000000 0.0000000 0.9466303
## y 0.0000000 0.0000000 0.9466303
## z 0.9466303 0.9466303 0.0000000Ya que se está utilizando una medida de 1 - Correlación, se está midiendo disimilaridad. Entre más cercano el valor a 0 indica menor disimilaridad, mostrando que x y Y son los menos disimilares, y por ende podrían formar un grupo. Esta es una de las maneras de agrupar (por similaridad o disimilaridad). Ahora se realizará un método de clasificación con datos más complejos.
Análisis exploratorio de los datos
df <- data.frame(iris)
iris.2 <- iris[,-5]
species <- iris[,5]
di=data.frame(df) %>%
mutate(Species=dplyr::recode(Species,
setosa="st",
versicolor="vs",
virginica="vg"))
pairs(di[,1:4], col = species,lower.panel = NULL)
par(xpd = TRUE)
legend(x = 0.05, y = 0.4, cex = 2,
legend=as.character(levels(species)),
fill = unique(species))
par(xpd = NA)En los resultados del análisis exploratorio se evidencia que para varias combinaciones de variables pueden formarse grupos diferenciados, especialmente de la especie setosa. Sabiendo esto, se justifica intentar correr los métodos de clasificación de clusters.
Método de K-means
Este algoritmo de clasificación no supervisada agrupa objetos en k grupos basándose en la mínima suma de distancias entre cada objeto y el centroide de su grupo o cluster.
set.seed(20)
k.means.fit <-kmeans(di[,1:4], 3, nstart = 10)
k.means.fit ## K-means clustering with 3 clusters of sizes 50, 62, 38
##
## Cluster means:
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 5.006000 3.428000 1.462000 0.246000
## 2 5.901613 2.748387 4.393548 1.433871
## 3 6.850000 3.073684 5.742105 2.071053
##
## Clustering vector:
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [71] 2 2 2 2 2 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2 3 3 3
## [106] 3 2 3 3 3 3 3 3 2 2 3 3 3 3 2 3 2 3 2 3 3 2 2 3 3 3 3 3 2 3 3 3 3 2 3
## [141] 3 3 2 3 3 3 2 3 3 2
##
## Within cluster sum of squares by cluster:
## [1] 15.15100 39.82097 23.87947
## (between_SS / total_SS = 88.4 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss"
## [5] "tot.withinss" "betweenss" "size" "iter"
## [9] "ifault"En la línea 103 se eligieron 3 clústers (centros) debido al análisis exploratorio que se realizó. Ya se ha realizado una primera clasificación, ¿pero qué tanta coincidencia tiene esta clasificación?
k.means.fit$centers## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 5.006000 3.428000 1.462000 0.246000
## 2 5.901613 2.748387 4.393548 1.433871
## 3 6.850000 3.073684 5.742105 2.071053k.means.fit$ifault## [1] 0grupos=k.means.fit$cluster
table(di$Species,grupos) #Matriz de confusión## grupos
## 1 2 3
## st 50 0 0
## vs 0 48 2
## vg 0 14 36dif=data.frame(di,grupos)
dif=data.frame(dif) %>%
mutate(grupos=dplyr::recode(grupos,
"3"="st",
"2"="vs",
"1"="vg"))
table(dif$grupos,dif$Species)##
## st vs vg
## st 0 2 36
## vg 50 0 0
## vs 0 48 14Al observar la matriz de confusión, se observa que la especie setosa quedó totalmente clasificada en un grupo, 48 de los datos de versicolor fueron agrupados mientras que dos datos quedaron dentro de otro grupo. Virginica es el que presenta datos más confusos, ya que 16 de los 50 datos se encuentran en otro grupo junto con Versicolor.
Combinación K-means y PCA
A pesar de que la agrupación de variables puede ayudar en el análisis de la totalidad de datos, el proceso se hace cada vez más demandante entre mayor sea el número de variables. Por esta razón el PCA se hace una herramienta importante para identificar las variables que mayor aportan a la variabilidad de datos y trabajar con estas únicamente.
d2 <- scale(di[,1:4])
rownames(d2) <- di$Species
fviz_nbclust(x = d2, FUNcluster = kmeans, method = "wss", k.max = 15,
diss = get_dist(d2, method = "euclidean"), nstart = 50)
Esta figura de números óptimos de clusters indica que, con los datos actuales, puede hacerse una clasificación de máximo 15 clusters, aunque ya se ha evidenciado que tres clusters pueden ser suficientes.
#Aplicación del algoritmo K-means
x|set.seed(123)## logical(0)d2f=data.frame(d2)
km_clusters <- kmeans(x = d2f, centers = 3, nstart = 50)
# Las funciones del paquete factoextra emplean el nombre de las filas del
# dataframe que contiene los datos como identificador de las observaciones.
# Esto permite añadir labels a los gráficos.
fviz_cluster(object = km_clusters, data = d2f, show.clust.cent = TRUE,
ellipse.type = "euclid", star.plot = TRUE, repel = TRUE,
pointsize=0.5,outlier.color="darkred") +
labs(title = "Resultados clustering K-means") +
theme_bw() + theme(legend.position = "none")
En esta imagen se combinan entonces el método K-Means con PCA. Los dos componentes graficados representan 95% de la variabilidad de los datos, por lo que es más que suficiente trabajar con estos dos elementos.
Se puede realizar unos ajustes para mejorar la presentación de los clusters.
require(cluster)
pam.res <- pam(d2f, 3)
# Visualización
fviz_cluster(pam.res, geom = "point", ellipse.type = "norm",
show.clust.cent = TRUE,star.plot = TRUE)+
labs(title = "Resultados clustering K-means")+ theme_bw()
Con esta imagen se evidencia que hay cierto grado de conflicto entre los clusters 2 y 3, ya que tienen cierta área solapada entre ambos. Estas representaciones mezclan cluster y PCA, pero no indican el grado de representación de cada variable en cada uno de los componentes.
Biplot PCA y K-Means para medir representatividad
data(iris)
# PCA
pca <- prcomp(iris[,-5], scale=TRUE)
df.pca <- pca$x
# Cluster over the three first PCA dimensions
kc <- kmeans(df.pca[,1:3], 3)
fviz_pca_biplot(pca, label="var", habillage=as.factor(kc$cluster)) +
labs(color=NULL) + ggtitle("") +
theme(text = element_text(size = 15),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
legend.key = element_rect(fill = "white"))
Gracias a esta gráfica ahora se sabe que la longitud de pétalo es la variable más reresentada por el primer componente, al igual que el ancho de pétalo. En el caso del componente dos, la variable ancho de pétalo es la que está más representada. No obstante, hay algunos datos que presentan inconsistencia, especialmente en la zona de conflicto entre los clusters. Vamos a probar un método más robusto que no se vea afectado por valores extraños.
K-Medioides con algoritmo PAM
fviz_nbclust(x = d2f[,1:4], FUNcluster = pam, method = "wss", k.max = 15,
diss = dist(d2, method = "manhattan"))
set.seed(123)
pam_clusters <- pam(x = d2f[,1:4], k = 3, metric = "manhattan")
pam_clusters$medoids## Sepal.Length Sepal.Width Petal.Length Petal.Width
## st.7 -1.0184372 0.7861738 -1.2791040 -1.3110521
## vg.16 0.7930124 -0.1315388 0.9868021 0.7880307
## vs.44 -0.2938574 -0.8198233 0.2503826 0.1320673fviz_cluster(object = pam_clusters, data = d2f[,1:4],
ellipse.type = "t",repel = TRUE) +
theme_bw() + labs(title = "Resultados clustering PAM") +
theme(legend.position = "none")
medoids <- prcomp(d2f[,1:4])$x
medoids <- medoids[rownames(pam_clusters$medoids), c("PC1", "PC2")]
medoids <- as.data.frame(medoids)
colnames(medoids) <- c("x", "y")
# Creación del gráfico
fviz_cluster(object = pam_clusters, data = d2f[,1:4], ellipse.type = "t",
repel = TRUE) + theme_bw() +
# Se resaltan las observaciones que actúan como medoids
geom_point(data = medoids, color = "firebrick", size = 2) +
labs(title = "Resultados clustering PAM") +
theme(legend.position = "none")
El punto rojo en cada cluster es su medioide, a partir del cual se hace la clasificación.
Método no jerárquico CLARA basado en simulación y muestreo
Cuando los datos a clasificar son demasiados, el algoritmo K-means o K-medoides exige una alta capacidad computacional para poder desarrollarlo. El método CLARA (Clustering Large Applications) combina la idea de K-medoides con el remuestreo y así poder trabajar con grandes volúmenes de datos de manera más eficiente.
Para ver las aplicaciones de este método se usaran nuevos datos relacionados con el color del café registrados en cinco medidas y dos fechas de maduración. El algoritmo que se utilizará es PAM.
set.seed(555)
dt1 <- genCorData(200, mu = c(37.88,16.79, 11.74,21.4,34.04),
sigma = 1, rho = 0.5, corstr = "cs" )
#Con la función rho se generan datos aleatorios CORRELACIONADOS, y no dependientes como se generan con los paquetes básicos.
dt2 <- genCorData(200, mu = c(35.94,17.90, 11.81,21.77,32.86),
sigma = 1, rho = 0.7, corstr = "cs" )
dac <- rbind(dt1,dt2)
DDA=gl(2,200,400,labels=c("DDA224","DDA231"))
dfc=data.frame(DDA,dac)
colnames(dfc) <- c("DDA","id","L", "a","b","C","h")
dfc=dfc[,-2]
head(dfc)## DDA L a b C h
## 1 DDA224 36.71324 16.59757 11.88936 20.80325 34.02769
## 2 DDA224 38.31073 15.76697 13.00676 21.46628 35.65491
## 3 DDA224 38.66201 17.00284 12.14174 22.81913 34.63488
## 4 DDA224 37.53979 17.02302 14.26060 21.62405 34.72869
## 5 DDA224 38.97227 15.98841 12.11759 20.35450 34.15341
## 6 DDA224 38.61491 17.14862 11.28327 22.67647 34.46416options(digits = 3)
dfc %>%
split(.$DDA) %>%
map(select, -c(DDA)) %>%
map(cor) ## $DDA224
## L a b C h
## L 1.000 0.468 0.415 0.535 0.417
## a 0.468 1.000 0.463 0.505 0.454
## b 0.415 0.463 1.000 0.391 0.437
## C 0.535 0.505 0.391 1.000 0.526
## h 0.417 0.454 0.437 0.526 1.000
##
## $DDA231
## L a b C h
## L 1.000 0.664 0.681 0.699 0.680
## a 0.664 1.000 0.714 0.683 0.735
## b 0.681 0.714 1.000 0.686 0.696
## C 0.699 0.683 0.686 1.000 0.707
## h 0.680 0.735 0.696 0.707 1.000clara_clusters <- clara(x = dfc, k = 2, metric = "manhattan", stand = TRUE,
samples = 60, pamLike = TRUE)
clara_clusters$sample## [1] 2 9 16 25 26 33 37 41 61 87 90 94 95 107 110 115 116
## [18] 141 142 146 148 150 159 161 183 193 204 222 223 228 244 249 251 253
## [35] 256 257 282 284 296 302 304 314 325 380clara_clusters$medoids## DDA L a b C h
## [1,] 1 37.5 17 11.7 21.2 34.4
## [2,] 2 36.1 18 11.8 21.9 32.6clara_clusters$i.med## [1] 94 302clara_clusters$clustering## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [71] 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2
## [106] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [141] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [176] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2
## [211] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [246] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [281] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [316] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [351] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [386] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2table(clara_clusters$clustering)##
## 1 2
## 197 203fviz_cluster(object = clara_clusters, ellipse.type = "t", geom = "point",
pointsize = 1.5) + theme_bw() +
labs(title = "Resultados clustering CLARA") +
theme(legend.position = "none")
medoides <- prcomp(dfc[,-1])$x
medoides <- medoides[rownames(pam_clusters$medoides), c("PC1", "PC2")]
medoides <- as.data.frame(medoides)
colnames(medoides) <- c("x", "y")
# Creación del gráfico
fviz_cluster(object = clara_clusters, data = d2f[,-1], ellipse.type = "t",
repel = TRUE,pointsize = 0.5) + theme_bw() +
# Se resaltan las observaciones que actúan como medoids
geom_point(data = medoides, color = "firebrick", size = 1,) +
labs(title = "Resultados clustering PAM") +
theme(legend.position = "none")
La designación de estos grupos se debe a las fechas de maduración en que fueron tomadas las cinco variables. Los 2000 datos fueron categorizados rápidamente por la mecánica del método empleado. ¿Cómo podemos saber cuál de las fechas de maduración del café quedó mejor representado en los dos primeros componentes.
Biplot CLARA con PCA
# PCA
pca <- prcomp(dfc[,-1], scale=TRUE)
df.pca <- pca$x
# Cluster over the three first PCA dimensions
CLA <- clara(df.pca[,1:3], k = 2, metric = "manhattan", stand = TRUE,
samples = 60, pamLike = TRUE)
fviz_pca_biplot(pca, label="var", habillage=as.factor(CLA$cluster)) +
labs(color=NULL) + ggtitle("") +
theme(text = element_text(size = 15),
panel.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black"),
legend.key = element_rect(fill = "white"))
Con este biplot sabemos que la variable “b” es el más representativo para el primer componente y que la variable a es la más representativa para el segundo componente.
Los métodos que se han presentado hasta ahora son no jerárquicos de Reasignación por centroides (K-means) y por medoides (K-medoides y Clara). Estos métodos categorizan los elementos según un número de Cluster dado. Ahora se presentarán una serie de métodos Jerárquicos que van generando grupos en cada una de las fases del proceso, buscando el número de cluster que hacen una agrupación óptima.
Métodos de clasificación Jerárquica
datj <- scale(dfc[,-1])
rownames(datj) <- dfc[,1]
# Matriz de distancias euclideas
mat_dist <- dist(x = datj, method = "euclidean")
# Dendrogramas con linkage complete y average
hc_euclidea_complete <- hclust(d = mat_dist, method = "complete")
hc_euclidea_average <- hclust(d = mat_dist, method = "average")
hc_euclidea_single <- hclust(d = mat_dist, method = "single")
hc_euclidea_ward.D2 <- hclust(d = mat_dist, method = "ward.D2")
hc_euclidea_median <- hclust(d = mat_dist, method = "median")
hc_euclidea_centroid <- hclust(d = mat_dist, method = "centroid")
hc_euclidea_mcquitty <- hclust(d = mat_dist, method = "mcquitty")
cor(x = mat_dist, cophenetic(hc_euclidea_complete))## [1] 0.486cor(x = mat_dist, cophenetic(hc_euclidea_average))## [1] 0.618cor(x = mat_dist, cophenetic(hc_euclidea_single))## [1] 0.435cor(x = mat_dist, cophenetic(hc_euclidea_ward.D2))## [1] 0.502cor(x = mat_dist, cophenetic(hc_euclidea_median))## [1] 0.381cor(x = mat_dist, cophenetic(hc_euclidea_centroid)) ## Tiene la mayor correlación## [1] 0.538cor(x = mat_dist, cophenetic(hc_euclidea_mcquitty))## [1] 0.481De los siete métodos evaluados, el método de centroide es el que tiene mayor correlación (0.538), por lo que es la mejor opción para hacer la clasificación de los datos
set.seed(101)
hc_euclidea_av <- hclust(d = dist(x = datj, method = "euclidean"),
method = "centroid")
fviz_dend(x = hc_euclidea_av, k = 2, cex = 0.5,
k_colors = c("red","green"),color_labels_by_k = T,
lwd = 0.2,type = "c",label_cols = rainbow(400),
rect_lty = "lightblue") +
geom_hline(yintercept = 3.65, linetype = "dashed") +
labs(title = "Herarchical clustering",
subtitle = "Distancia euclídea, Centroide, k=2")
También es posible cambiar la visiualización a árbol de decisión cambiando el argumento type= “c” por “r”.
set.seed(101)
hc_euclidea_av <- hclust(d = dist(x = datj, method = "euclidean"),
method = "centroid")
fviz_dend(x = hc_euclidea_av, k = 2, cex = 0.5,
k_colors = c("red","green"),color_labels_by_k = T,
lwd = 0.2,type = "r",label_cols = rainbow(400),
rect_lty = "lightblue") +
geom_hline(yintercept = 3.65, linetype = "dashed") +
labs(title = "Herarchical clustering",
subtitle = "Distancia euclídea, Centroide, k=2")
Este árbol de decisión, a pesar de tener la mayor correlación, no es un árbol adecuado ya que no se definen claramente los clusters. Se va a probar ahora con el método “Average”.
hc_euclidea_av <- hclust(d = dist(x = datj, method = "euclidean"),
method = "average")
fviz_dend(x = hc_euclidea_av, k = 2, cex = 0.5,
k_colors = c("red","green"),color_labels_by_k = T,
lwd = 0.2,type = "c",label_cols = rainbow(400),
rect_lty = "lightblue") +
geom_hline(yintercept = 3.65, linetype = "dashed") +
labs(title = "Herarchical clustering",
subtitle = "Distancia euclídea, Average, k=2")
####################################### Otra visualización (línea 231 "r")
hc_euclidea_av <- hclust(d = dist(x = datj, method = "euclidean"),
method = "average")
fviz_dend(x = hc_euclidea_av, k = 2, cex = 0.5,
k_colors = c("red","green"),color_labels_by_k = T,
lwd = 0.2,type = "r",label_cols = rainbow(400),
rect_lty = "lightblue") +
geom_hline(yintercept = 3.65, linetype = "dashed") +
labs(title = "Herarchical clustering",
subtitle = "Distancia euclidea, Average, k=2")
Con este método si se evidencia la clasificación en dos grandes clusters, ajustado mejor a estos datos y brindando información de manera más adecuada.
Bibliografía
Everitt Brian, Landau Sabine, Leese Morven & Stahl Daniel (2011). Cluster Analysis. Wiley Series in Probability and Statistics. 337 pp.