Clustering y heatmaps: aprendizaje no supervisado
Joaquín Amat Rodrigo
Septiembre, 2017
Más sobre ciencia de datos en cienciadedatos.net
Versión PDF: Github
Introducción
El término clustering hace referencia a un amplio abanico de técnicas no supervisadas cuya finalidad es encontrar patrones o grupos (clusters) dentro de un conjunto de observaciones. Las particiones se establecen de forma que, las observaciones que están dentro de un mismo grupo, son similares entre ellas y distintas a las observaciones de otros grupos. Se trata de un método no supervisado, ya que el proceso ignora la variable respuesta que indica a que grupo pertenece realmente cada observación (si es que existe tal variable). Esta característica diferencia al clustering de las técnicas supervisadas, que emplean un set de entrenamiento en el que se conoce la verdadera clasificación.
Dada la utilidad del clustering en disciplinas muy distintas (genómica, marketing…), se han desarrollado multitud de variantes y adaptaciones de sus métodos y algoritmos. Pueden diferenciarse tres grupos principales:
Partitioning Clustering: Este tipo de algoritmos requieren que el usuario especifique de antemano el número de clusters que se van a crear (K-means, K-medoids, CLARA).
Hierarchical Clustering: Este tipo de algoritmos no requieren que el usuario especifique de antemano el número de clusters. (agglomerative clustering, divisive clusterig).
Métodos que combinan o modifican los anteriores (hierarchical K-means, fuzzy clustering, model based clustering y density based clustering).
En el entorno de programación R existen múltiples paquetes que implementan algoritmos de clustering y funciones para visualizar sus resultados. En este documento se emplean los siguientes:
stats: contiene las funcionesdist()para calcular matrices de distancias,kmeans(),hclust(),cuttree()para crear los clusters yplot.hclust()para visualizar los resultados.cluster,mclust: contienen múltiples algoritmos de clustering y métricas para evaluarlos.factoextra: extensión basada enggplot2para crear visualizaciones de los resultados de clustering y su evaluación.dendextend: extensión para la customización de dendrogramas.
Medidas de distancia
Todos los métodos de clustering tienen una cosa en común, para poder llevar a cabo las agrupaciones necesitan definir y cuantificar la similitud entre las observaciones. El término distancia se emplea dentro del contexto del clustering como cuantificación de la similitud o diferencia entre observaciones. Si se representan las observaciones en un espacio p dimensional, siendo p el número de variables asociadas a cada observación, cuando más se asemejen dos observaciones más próximas estarán, de ahí que se emplee el término distancia. La característica que hace del clustering un método adaptable a escenarios muy diversos es que puede emplear cualquier tipo de distancia, lo que permite al investigador escoger la más adecuada para el estudio en cuestión. A continuación, se describen algunas de las más utilizadas.
Distancia euclídea
La distancia euclídea entre dos puntos p y q se define como la longitud del segmento que une ambos puntos. En coordenadas cartesianas, la distancia euclídea se calcula empleando el teorema de Pitágoras. Por ejemplo, en un espacio de dos dimensiones en el que cada punto está definido por las coordenadas \((x,y)\), la distancia euclídea entre p y q viene dada por la ecuación:
\[d_{euc}(p,q) = \sqrt{(x_p - x_q)^2 + (y_p - y_q)^2}\]
Esta ecuación puede generalizarse para un espacio euclídeo n-dimensional donde cada punto está definido por un vector de n coordenadas: \(p = (p_1,p_2,p_3,...,p_n)\) y \(q = (q_1,q_2,q_3,...,q_n)\).
\[d_{euc}(p,q) = \sqrt{(p_1 - q_1)^2 + (p_2 - q_2)^2 +...+(p_n - q_n)^2} =\]
\[\sqrt{ \sum_{i=1}^n (p_i - q_i)^2}\] Una forma de dar mayor peso a aquellas observaciones que están más alejadas es emplear la distancia euclídea al cuadrado. En el caso del clustering, donde se busca agrupar observaciones que minimicen la distancia, esto se traduce en una mayor influencia de aquellas observaciones que están más distantes.
La siguiente imagen muestra el perfil de dos observaciones definidas por 10 variables (espacio con 10 dimensiones).
library(ggplot2)
observacion_a <- c(4, 4.5, 4, 7.5, 7, 6, 5, 5.5, 5, 6)
observacion_b <- c(4, 4.5, 4, 7.5, 7, 6, 5, 5.5, 5, 6) + 4
datos <- data.frame(observacion = rep(c("a", "b"), each = 10),
valor = c(observacion_a, observacion_b),
predictor = 1:10)
ggplot(data = datos, aes(x = as.factor(predictor), y = valor,
colour = observacion)) +
geom_path(aes(group = observacion)) +
geom_point() +
geom_line(aes(group = predictor), colour = "firebrick", linetype = "dashed") +
labs(x = "predictor") +
theme_bw() +
theme(legend.position = "none")
La distancia euclídea entre las dos observaciones equivale a la raíz cuadrada de la suma de las longitudes de los segmentos rojos que unen cada par de puntos. Tiene en cuenta por lo tanto el desplazamiento individual de cada una de las variables.
Distancia de Manhattan
La distancia de Manhattan, también conocida como taxicab metric, rectilinear distance o L1 distance, define la distancia entre dos puntos p y q como el sumatorio de las diferencias absolutas entre cada dimensión. Esta medida se ve menos afectada por outliers (es más robusta) que la distancia euclídea debido a que no eleva al cuadrado las diferencias.
\[d_{man}(p,q) = \sum_{i=1}^n |(p_i - q_i)|\]
La siguiente imagen muestra una comparación entre la distancia euclídea (segmento azul) y la distancia de manhattan (segmento rojo y verde) en un espacio bidimensional. Existen múltiples caminos para unir dos puntos con el mismo valor de distancia de manhattan, ya que su valor es igual al desplazamiento total en cada una de las dimensiones.
datos <- data.frame(observacion = c("a", "b"), x = c(2,7), y = c(2,7))
manhattan <- data.frame(
x = rep(2:6, each = 2),
y = rep(2:6, each = 2) + rep(c(0,1), 5),
xend = rep(2:6, each = 2) + rep(c(0,1), 5),
yend = rep(3:7, each = 2))
manhattan_2 <- data.frame(
x = c(2, 5, 5, 7),
y = c(2, 2, 4, 4),
xend = c(5, 5, 7, 7),
yend = c(2, 4, 4, 7))
ggplot(data = datos, aes(x = x, y = y)) +
geom_segment(aes(x = 2, y = 2, xend = 7, yend = 7), color = "blue", size = 1.2) +
geom_segment(data = manhattan, aes(x = x, y = y, xend = xend, yend = yend),
color = "red", size = 1.2) +
geom_segment(data = manhattan_2, aes(x = x, y = y, xend = xend, yend = yend),
color = "green3", size = 1.2) +
geom_point(size = 3) +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_line(size = 2),
panel.background = element_rect(fill = "gray",
colour = "white",
size = 0.5, linetype = "solid"))
Correlación
La correlación es una medida de distancia muy útil cuando la definición de similitud se hace en términos de patrón o forma y no de desplazamiento o magnitud. ¿Qué quiere decir esto? En la imagen del apartado de la distancia euclídea, las dos observaciones tienen exactamente el mismo patrón, la única diferencia es que una de ellas está desplazada 4 unidades por encima de la otra. Si se emplea como medida de similitud 1 menos el valor de la correlación, ambas observaciones se consideran idénticas (su distancia es 0).
\[d_{cor}(p,q) = 1 - \text{correlacion}(p,q)\] donde la correlación puede ser de distintos tipos (Pearson, Spearman, Kendall…) Correlación lineal.
En la siguiente imagen se muestra el perfil de 3 observaciones. Acorde a la distancia euclídea, las observaciones b y c son las más similares, mientras que acorde a la correlación de Pearson, las observaciones más similares son a y b.
library(ggplot2)
observacion_a <- c(4, 4.5, 4, 7.5, 7, 6, 5, 5.5, 5, 6)
observacion_b <- c(4, 4.5, 4, 7.5, 7, 6, 5, 5.5, 5, 6) + 4
observacion_c <- c(5, 5.5, 4.8, 5.4, 4.7, 5.6, 5.3, 5.5, 5.2, 4.8)
datos <- data.frame(observacion = rep(c("a", "b", "c"), each = 10),
valor = c(observacion_a, observacion_b, observacion_c),
predictor = 1:10)
ggplot(data = datos, aes(x = as.factor(predictor), y = valor,
colour = observacion)) +
geom_path(aes(group = observacion)) +
geom_point() +
labs(x = "predictor") +
theme_bw() +
theme(legend.position = "bottom")
## observacion_a observacion_b
## observacion_b 12.64911
## observacion_c 3.75100 13.98821## observacion_a observacion_b observacion_c
## observacion_a 0.0000000 0.0000000 0.9466303
## observacion_b 0.0000000 0.0000000 0.9466303
## observacion_c 0.9466303 0.9466303 0.0000000Este ejemplo pone de manifiesto que no existe una única medida de distancia que sea mejor que las demás, sino que, dependiendo del contexto, una será más adecuada que otra.
Jackknife correlation
El coeficiente de correlación de Pearson resulta efectivo en ámbitos muy diversos. Sin embargo, tiene la desventaja de no ser robusto frente a outliers a pesar de que se cumpla la condición de normalidad. Si dos variables tienen un pico o un valle común en una única observación, por ejemplo, por un error de lectura, la correlación va a estar dominada por este registro a pesar de que entre las dos variables no haya correlación real alguna. Lo mismo puede ocurrir en la dirección opuesta. Si dos variables están altamente correlacionadas excepto para una observación, en la que los valores son muy dispares, entonces la correlación existente quedará enmascarada. Una forma de evitarlo es recurrir a la Jackknife correlation, que consiste en calcular todos los posibles coeficientes de correlación entre dos variables si se excluye cada vez una de las observaciones. El promedio de todas las Jackknife correlations calculadas atenua en cierta medida el efecto del outlier.
\[\bar{\theta}_{(A,B)} = \text{Promedio Jackknife correlation (A,B)} = \frac{1}{n} \sum_{i=1}^n \hat r_i\]
donde \(n\) es el número de observaciones y \(\hat r_i\) es el coeficiente de correlación entre las variables \(A\) y \(B\), habiendo excluido la observación \(i\).
Además del promedio, se puede estimar su error estándar (\(SE\)) y así obtener intervalos de confianza para la Jackknife correlation y su correspondiente p-value.
\[SE = \sqrt{\frac{n-1}{n} \sum_{i=1}^n (\hat r_i - \bar{\theta})^2}\]
Intervalo de confianza del 95% \((Z=1.96)\)
\[\text{Promedio Jackknife correlation (A,B)} \pm \ 1.96 *SE\]
\[\bar{\theta} - 1.96 \sqrt{\frac{n-1}{n} \sum_{i=1}^n (\hat r_i - \bar{\theta})^2} \ , \ \ \bar{\theta}+ 1.96 \sqrt{\frac{n-1}{n} \sum_{i=1}^n (\hat r_i - \bar{\theta})^2}\]
P-value para la hipótesis nula de que \(\bar{\theta} = 0\):
\[Z_{calculada} = \frac{\bar{\theta} - H_0}{SE}= \frac{\bar{\theta} - 0}{\sqrt{\frac{n-1}{n} \sum_{i=1}^n (\hat r_i - \bar{\theta})^2}}\]
\[p_{value} = P(Z > Z_{calculada})\]
Cuando se emplea este método, es conveniente calcular la diferencia entre el valor de correlación obtenido por Jackknife correlation \((\bar{\theta})\) y el que se obtiene si se emplean todas las observaciones \((\bar{r})\). A esta diferencia se le conoce como Bias. Su magnitud es un indicativo de cuanto está influenciada la estimación de la correlación entre dos variables debido a un valor atípico u outlier.
\[Bias = (n-1)*(\bar{\theta} - \hat{r})\]
Si se calcula la diferencia entre cada correlación \((\hat r_i)\) estimada en el proceso de Jackknife y el valor de correlación \((\hat r)\) obtenido si se emplean todas las observaciones, se puede identificar que observaciones son más influyentes.
Cuando el estudio requiere minimizar al máximo la presencia de falsos positivos, a pesar de que se incremente la de falsos negativos, se puede seleccionar como valor de correlación el menor de entre todos los calculados en el proceso de Jackknife.
\[\text{Correlacion} = \text{min} \{ \hat r_1, \hat r_2,..., \hat r_n \}\]
A pesar de que el método de Jackknife permite aumentar la robustez de la correlación de Pearson, si los outliers son muy extremos su influencia seguirá siendo notable.
Simple matching coefficient
Cuando las variables con las que se pretende determinar la similitud entre observaciones son de tipo binario, a pesar de que es posible codificarlas de forma numérica como \(1\) o \(0\), no tiene sentido aplicar operaciones aritméticas sobre ellas (media, suma…). Por ejemplo, si la variable sexo se codifica como \(1\) para mujer y \(0\) para hombre, carece de significado decir que la media de la variable sexo en un determinado set de datos es \(0.5\). En situaciones como esta, no se pueden emplear medidas de similitud basadas en distancia euclídea, manhattan, correlación…
Dado dos objetos A y B, cada uno con \(n\) atributos binarios, el simple matching coefficient (SMC) define la similitud entre ellos como:
\[SMC= \frac{\text{número coincidencias}}{\text{número total de atributos}} = \frac{M_{00} + M{11}}{M_{00} + M_{01}+ M_{10}+M_{11}}\]
donde \(M_{00}\) y \(M_{11}\) son el número de variables para las que ambas observaciones tienen el mismo valor (ambas \(0\) o ambas \(1\)), y \(M_{01}\) y \(M_{10}\) el número de variables que no coinciden. El valor de distancia simple matching distance (SMD) se corresponde con \(1 - SMC\).
Índice Jaccard
El índice Jaccard o coeficiente de correlación Jaccard es similar al simple matching coefficient (SMC). La diferencia radica en que el SMC tiene el término \(M_{00}\) en el numerador y denominador, mientras que el índice de Jaccard no. Esto significa que SMC considera como coincidencias tanto si el atributo está presente en ambos sets como si el atributo no está en ninguno de los sets, mientras que Jaccard solo cuenta como coincidencias cuando el atributo está presente en ambos sets.
\[J= \frac{M_{11}}{M_{01}+ M_{10}+M_{11}}\] o en términos matemáticos de sets:
\[J(A,B) = \frac{A \cap B}{A \cup B}\]
La distancia de Jaccard (\(1-J\)) supera a la simple matching distance en aquellas situaciones en las que la coincidencia de ausencia no aporta información. Para ilustrar este hecho, supóngase que se quiere cuantificar la similitud entre dos clientes de un supermercado en base a los artículos comprados. Es de esperar que cada cliente solo adquiera unos pocos artículos de los muchos disponibles, por lo que el número de artículos no comprados por ninguno (\(M_{00}\)) será muy alto. Como la distancia de Jaccard ignora las coincidencias de tipo \(M_{00}\), el grado de similitud dependerá únicamente de las coincidencias entre los artículos comprados.
Distancia coseno
El coseno del ángulo que forman dos vectores puede interpretarse como una medida de similitud de sus orientaciones, independientemente de sus magnitudes. Si dos vectores tienen exactamente la misma orientación (el ángulo que forman es 0º) su coseno toma el valor de 1, si son perpendiculares (forman un ángulo de 90º) su coseno es 0 y si tienen orientaciones opuestas (ángulo de 180º) su coseno es de -1.
\[\text{cos}(\alpha) = \frac{\textbf{x} \cdot \textbf{y}}{||\textbf{x}|| \ ||\textbf{y}||}\]
Si los vectores se normalizan restándoles la media (\(X − \bar{X}\)), la medida recibe el nombre de coseno centrado y es equivalente a la correlación de Pearson.
a <- c(4, 4.5, 4, 7, 7, 6, 5, 5.5, 5, 6)
b <- c(5, 5.5, 4.8, 5.4, 4.7, 5.6, 5.3, 5.5, 5.2, 4.8)
# Correlación de Pearson
cor(x = a, y = b, method = "pearson")## [1] 0.02427991# Coseno
coseno <- function(x, y){
resultado <- x%*%y / (sqrt(x %*% x) * sqrt(y %*%y ))
return(as.numeric(resultado))
}
coseno(a,b)## [1] 0.9802813## [1] 0.02427991Escala de las variables
Al igual que en otros métodos estadísticos (PCA, ridge regression, lasso…), la escala en la que se miden las variables y la magnitud de su varianza pueden afectar en gran medida a los resultados obtenidos por clustering. Si una variable tiene una escala mucho mayor que el resto, determinará en gran medida el valor de distancia/similitud obtenido al comparar las observaciones, dirigiendo así la agrupación final. Escalar y centrar las variables de forma que todas ellas tengan media 0 y desviación estándar 1 antes de calcular la matriz de distancias asegura que todas las variables tengan el mismo peso cuando se realice el clustering. Sebastian Raschka hace un análisis muy explicativo de los distintos tipos de escalado y normalización.
\[\frac{x_i - mean(x)}{sd(x)}\]
Para ilustrar este hecho, supóngase que una tienda online quiere clasificar a los compradores en función de los artículos que adquieren, por ejemplo, calcetines y ordenadores. La siguiente imagen muestra el número de artículos comprados por 8 clientes a lo largo de un año, junto con el gasto total de cada uno.
Si se intenta agrupar a los clientes por el número de artículos comprados, dado que los calcetines se compran con mucha más frecuencia que los ordenadores, van a tener más peso al crear los clusters. Por el contrario, si la agrupación se hace en base al gasto total de los clientes, como los ordenadores son mucho más caros, van a determinar en gran medida la clasificación. Escalando y centrando las variables se consigue igualar la influencia de calcetines y ordenadores.
Cabe destacar que, si se aplica la estandarización descrita, existe una relación entre la distancia euclídea y la correlación de Pearson que hace que los resultados obtenidos por clustering en ambos casos sean equivalentes.
\[d_{euc}(p,q\ estandarizados) = \sqrt{2n(1-cor(p,q ))}\]
Ejemplo
El set de datos USArrests contiene información sobre el número de delitos (asaltos, asesinatos y secuestros) junto con el porcentaje de población urbana para cada uno de los 50 estados de USA. Empleando estas variables se pretende calcular una matriz de distancias que permita identificar los estados más similares
Dos de las funciones en R que permiten calcular matrices de distancia empleando variables numéricas son dist() y get_dist(). Esta última incluye más tipos de distancias.
data(USArrests)
# Se escalan las variables
datos <- scale(USArrests, center = TRUE, scale = TRUE)
# Distancia euclídea
mat_dist <- dist(x = datos, method = "euclidean")
round(as.matrix(mat_dist)[1:5, 1:5], 2)## Alabama Alaska Arizona Arkansas California
## Alabama 0.00 2.70 2.29 1.29 3.26
## Alaska 2.70 0.00 2.70 2.83 3.01
## Arizona 2.29 2.70 0.00 2.72 1.31
## Arkansas 1.29 2.83 2.72 0.00 3.76
## California 3.26 3.01 1.31 3.76 0.00# Distancia basada en la correlación de pearson
library(factoextra)
mat_dist <- get_dist(x = datos, method = "pearson")
round(as.matrix(mat_dist)[1:5, 1:5], 2)## Alabama Alaska Arizona Arkansas California
## Alabama 0.00 0.71 1.45 0.09 1.87
## Alaska 0.71 0.00 0.83 0.37 0.81
## Arizona 1.45 0.83 0.00 1.18 0.29
## Arkansas 0.09 0.37 1.18 0.00 1.59
## California 1.87 0.81 0.29 1.59 0.00# Esto es equivalente a 1 - correlación pearson
round(1 - cor(x = t(datos), method = "pearson"), 2)[1:5, 1:5]## Alabama Alaska Arizona Arkansas California
## Alabama 0.00 0.71 1.45 0.09 1.87
## Alaska 0.71 0.00 0.83 0.37 0.81
## Arizona 1.45 0.83 0.00 1.18 0.29
## Arkansas 0.09 0.37 1.18 0.00 1.59
## California 1.87 0.81 0.29 1.59 0.00La función fviz_dist del paquete factoextra resulta muy útil para generar representaciones gráficas (heatmap) de matrices de distancia.
K-means clustering
Idea intuitiva
El método K-means clustering (MacQueen, 1967) agrupa las observaciones en K clusters distintos, donde el número K lo determina el analista antes de ejecutar del algoritmo. K-means clustering encuentra los K mejores clusters, entendiendo como mejor cluster aquel cuya varianza interna (intra-cluster variation) sea lo más pequeña posible. Se trata por lo tanto de un problema de optimización, en el que se reparten las observaciones en K clusters de forma que la suma de las varianzas internas de todos ellos sea lo menor posible. Para poder solucionar este problema es necesario definir un modo de cuantificar la varianza interna.
Considérense \(C_1\),…, \(C_K\) como los sets formados por los índices de las observaciones de cada uno de los clusters. Por ejemplo, el set \(C_1\) contiene los índices de las observaciones agrupadas en el cluster 1. La nomenclatura empleada para indicar que la observación \(i\) pertenece al cluster \(k\) es: \(i \in C_k\). Todos los sets satisfacen dos propiedades:
\(C_1 \cup C_2 \cup ... \cup C_K = \{1,...,n\}\). Significa que toda observación pertenece al menos a uno de los K clusters.
\(C_k \cap C_{k'} = \emptyset\) para todo \(k \neq k'\). Implica que los clusters no solapan, ninguna observación pertenece a más de un cluster a la vez.
Dos de las medidas más comúnmente empleadas definen la varianza interna de un cluster (\(W(C_k)\)) como:
- La suma de las distancias euclídeas al cuadrado entre cada observación (\(x_i\)) y el centroide (\(\mu\)) de su cluster. Esto equivale a la suma de cuadrados internos del cluster.
\[W(C_k) = \sum_{x_i,\in C_k} (x_{i} - \mu_k)^2\]
- La suma de las distancias euclídeas al cuadrado entre todos los pares de observaciones que forman el cluster, dividida entre el número de observaciones del cluster.
\[W(C_k) = \frac{1}{|C_k|} \sum_{i,i' \in C_k} \sum^p_{j=1}(x_{ij} - x_{i'j})^2\]
Minimizar la suma total de varianza interna \(\sum^k_{k=1}W(C_k)\) de forma exacta (encontrar el mínimo global) es un proceso muy complejo debido a la inmensa cantidad de formas en las que n observaciones se pueden dividir en K grupos. Sin embargo, es posible obtener una solución que, aun no siendo la mejor de entre todas las posibles, es muy buena (óptimo local). El algoritmo empleado para ello es:
Asignar aleatoriamente un número entre 1 y \(K\) a cada observación. Esto sirve como asignación inicial aleatoria de las observaciones a los clusters.
Iterar los siguientes pasos hasta que la asignación de las observaciones a los clusters no cambie o se alcance un número máximo de iteraciones establecido por el usuario.
2.1 Para cada uno de los clusters calcular su centroide. Entendiendo por centroide la posición definida por la media de cada una de las dimensiones (variables) de las observaciones que forman el cluster. Aunque no es siempre equivalente, puede entenderse como el centro de gravedad.
2.2 Asignar cada observación al cluster cuyo centroide está más próximo.
Este algoritmo garantiza que, en cada paso, se reduzca la intra-varianza total de los clusters hasta alcanzar un óptimo local. La siguiente imagen muestra cómo van cambiando las asignaciones de las observaciones a medida que se ejecuta cada paso del algoritmo.
Otra forma de implementar el algoritmo de K-means clustering es la siguiente:
Especificar el número K de clusters que se quieren crear.
Seleccionar de forma aleatoria k observaciones del set de datos como centroides iniciales.
Asignar cada una de las observaciones al centroide más cercano.
Para cada uno de los K clusters recalcular su centroide.
Repetir los pasos 3 y 4 hasta que las asignaciones no cambien o se alcance el número máximo de iteraciones establecido.
Debido a que el algoritmo de K-means no evalúa todas las posibles distribuciones de las observaciones sino solo parte de ellas, los resultados obtenidos dependen de la asignación aleatoria inicial (paso 1). Por esta razón, es importante ejecutar el algoritmo varias veces (25-50), cada una con una asignación aleatoria inicial distinta, y seleccionar aquella que haya conseguido un menor valor de varianza total.
Ventajas y desventajas
K-means es uno de los métodos de clustering más utilizados. Destaca por la sencillez y velocidad de su algoritmo, sin embargo, presenta una serie de limitaciones que se deben tener en cuenta.
Requiere que se indique de antemano el número de clusters que se van a crear. Esto puede ser complicado si no se dispone de información adicional sobre los datos con los que se trabaja. Se han desarrollado varias estrategias para ayudar a identificar potenciales valores óptimos de K (ver más adelante), aunque todas ellas son orientativas.
Las agrupaciones resultantes pueden variar dependiendo de la asignación aleatoria inicial de los centroides. Para minimizar este problema se recomienda repetir el proceso de clustering entre 25-50 veces y seleccionar como resultado definitivo el que tenga menor suma total de varianza interna. Aun así, solo se puede garantizar la reproducibilidad de los resultados si se emplean semillas.
Presenta problemas de robustez frente a outliers. La única solución es excluirlos o recurrir a otros métodos de clustering más robustos como K-medoids (PAM).
Ejemplo 1
Los siguientes datos simulados contienen observaciones que pertenecen a cuatro grupos distintos. Se pretende aplicar K-means-clustering con el fin de identificarlos.
library(tidyverse)
library(ggpubr)
set.seed(101)
# Se simulan datos aleatorios con dos dimensiones
datos <- matrix(rnorm(n = 100*2), nrow = 100, ncol = 2,
dimnames = list(NULL,c("x", "y")))
datos <- as.data.frame(datos)
# Se determina la media que va a tener cada grupo en cada una de las dos
# dimensiones. En total 2*4 medias. Este valor se utiliza para separar
# cada grupo de los demás.
media_grupos <- matrix(rnorm(n = 8, mean = 0, sd = 4), nrow = 4, ncol = 2,
dimnames = list(NULL, c("media_x", "media_y")))
media_grupos <- as.data.frame(media_grupos)
media_grupos <- media_grupos %>% mutate(grupo = c("a","b","c","d"))
# Se genera un vector que asigne aleatoriamente cada observación a uno de
# los 4 grupos
datos <- datos %>% mutate(grupo = sample(x = c("a","b","c","d"),
size = 100,
replace = TRUE))
# Se incrementa el valor de cada observación con la media correspondiente al
# grupo asignado.
datos <- left_join(datos, media_grupos, by = "grupo")
datos <- datos %>% mutate(x = x + media_x,
y = y + media_y)
ggplot(data = datos, aes(x = x, y = y, color = grupo)) +
geom_point(size = 2.5) +
theme_bw()
La función kmeans() del paquete stats realiza K-mean-clustering. Entre sus argumentos destacan: centers, que determina el número K de clusters que se van a generar y nstart, que determina el número de veces que se va a repetir el proceso, cada vez con una asignación aleatoria inicial distinta. Es recomendable que este último valor sea alto, entre 25-50, para no obtener resultados malos debido a una iniciación poco afortunada del proceso.
Como los datos se han simulado considerando que todas las dimensiones tienen aproximadamente la misma magnitud, no es necesario escalarlos ni centrarlos. De no ser así, sí que habría que hacerlo.
## K-means clustering with 4 clusters of sizes 32, 20, 20, 28
##
## Cluster means:
## x y
## 1 -0.5787702 4.7639233
## 2 -3.1104142 1.2535711
## 3 1.4989983 -0.2412154
## 4 -5.6518323 3.3513316
##
## Clustering vector:
## [1] 4 2 1 4 2 1 4 2 1 1 3 1 1 3 2 3 4 3 4 4 4 4 4 3 1 1 2 4 2 1 4 3 4 2 2 3 3
## [38] 2 3 3 4 2 2 4 4 3 4 1 4 2 4 1 1 3 3 2 3 1 1 1 2 4 4 4 2 2 1 1 3 4 4 1 1 3
## [75] 1 3 4 1 1 1 2 1 2 1 4 3 1 4 4 1 1 2 4 2 1 1 3 3 1 1
##
## Within cluster sum of squares by cluster:
## [1] 53.04203 48.52107 34.95921 42.40322
## (between_SS / total_SS = 85.7 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"El objeto devuelto por la función kmeans() contiene entre otros datos: la media de cada una de las variables para cada cluster (centers), un vector indicando a que cluster se ha asignado cada observación (cluster), la suma de cuadrados interna de cada cluster (withinss) y la suma total de cuadrados internos de todos los clusters (tot.withinss). Al imprimir el resultado también se muestra y el ratio de la suma de cuadrados entre-clusters y la suma de cuadrados totales. Este ratio es equivalente al \(R^2\) de los modelos de regresión, indica el porcentaje de varianza explicada por el modelo respecto al total de varianza observada. Puede utilizarse para evaluar el clustering obtenido pero, al igual que ocurre con \(R^2\) al incrementar el número de predictores, el ratio between_SS / total_SS aumenta con el número de clusters creados. Esto último se debe de tener en cuenta para evitar problemas de overfitting.
Al tratarse de una simulación, se conoce el número real de grupos (4) y a cuál de ellos pertenece cada observación. Esto no sucede en la mayoría de casos prácticos, pero es útil para evaluar cómo de bueno es el método de K-means-clustering clasificando las observaciones.
# Se representa el número de cluster al que se ha asignado cada observación y
# se muestra con un código de color el grupo real al que pertenece.
datos <- datos %>% mutate(cluster = km_clusters$cluster)
datos <- datos %>% mutate(cluster = as.factor(cluster),
grupo = as.factor(grupo))
ggplot(data = datos, aes(x = x, y = y, color = grupo)) +
geom_text(aes(label = cluster), size = 5) +
theme_bw() +
theme(legend.position = "none")
El gráfico muestra que solo dos observaciones se han agrupado incorrectamente. Este tipo de visualización es muy útil e informativa, sin embargo, solo es posible cuando se trabaja con dos dimensiones y se conocen los verdaderos grupos a los que pertenecen las observaciones. Si los datos contienen más de dos variables (dimensiones), una posible solución es utilizar las dos primeras componentes principales obtenidas en un PCA previo.
El número de aciertos y errores puede representarse también en modo de matriz de confusión. A la hora de interpretar estas matrices, es importante recordar que el clustering asigna las observaciones a clusters cuyo identificador no tiene que por qué coincidir con la nomenclatura empleada para los grupos reales. Por ejemplo, el grupo b podría haberse llamado en su lugar grupo 2 y haberse asignado al cluster 1. Así pues, por cada fila de la matriz cabe esperar un valor alto (coincidencias) para una de las posiciones y valores bajos en las otras (errores de clasificación), pero no tienen por qué coincidir los nombres.
## grupo real
## cluster a b c d
## 1 31 0 0 1
## 2 1 0 0 19
## 3 0 0 20 0
## 4 0 28 0 0En este análisis, solo 2 de las 100 observaciones se han agrupado erróneamente. De nuevo repetir que, en la realidad, no se suelen conocer los verdaderos grupos en los que se dividen las observaciones, de lo contrario no se necesitaría aplicar clustering.
Supóngase ahora que se trata de un caso real, en el que se desconoce el número de grupos en los que se subdividen las observaciones. El investigador tendría que probar con diferentes valores de K y decidir cuál parece más razonable (en los siguientes apartados se describen métodos más sofisticados). A continuación, se muestran los resultados para \(K = 2\) y \(K = 6\).
# Resultados para K = 2
datos <- datos %>% select(x, y)
set.seed(101)
km_clusters_2 <- kmeans(x = datos, centers = 2, nstart = 50)
datos <- datos %>% mutate(cluster = km_clusters_2$cluster)
p1 <- ggplot(data = datos, aes(x = x, y = y, color = as.factor(cluster))) +
geom_point(size = 3) +
labs(title = "Kmenas con k=2") +
theme_bw() +
theme(legend.position = "none")
# Resultados para K = 6
datos <- datos %>% select(x, y)
set.seed(101)
km_clusters_6 <- kmeans(x = datos, centers = 6, nstart = 50)
datos <- datos %>% mutate(cluster = km_clusters_6$cluster)
p2 <- ggplot(data = datos, aes(x = x, y = y, color = as.factor(cluster))) +
geom_point(size = 3) +
labs(title = "Kmenas con k=6") +
theme_bw() +
theme(legend.position = "none")
ggarrange(p1, p2)
Al observar los resultados obtenidos para K = 2, es intuitivo pensar que el grupo que se encuentra entorno a las coordenadas \(x = -1, y = 6\) (mayoritariamente considerado como azul) debería ser un grupo separado. Para K = 6 no parece muy razonable la separación de los grupos rojo y verde. Este ejemplo muestra la principal limitación del método de K-means, el hecho de tener que escoger de antemano el número de clusters que se generan. El siguiente ejemplo, se muestra una de las estrategias para identificar posibles valores de K óptimos.
Ejemplo 2
El set de datos USArrests contiene información sobre el número de delitos (asaltos, asesinatos y secuestros) junto con el porcentaje de población urbana para cada uno de los 50 estados de USA. Se pretende estudiar si existe una agrupación subyacente de los estados empleando K-means-clustering.
El paquete factoextra creado por Alboukadel Kassambara contiene funciones que facilitan en gran medida la visualización y evaluación de los resultados de clustering.
Si se emplea K-means-clustering con distancia euclídea hay que asegurarse de que las variables empleadas son de tipo continuo, ya que trabaja con la media de cada una de ellas.
## 'data.frame': 50 obs. of 4 variables:
## $ Murder : num 13.2 10 8.1 8.8 9 7.9 3.3 5.9 15.4 17.4 ...
## $ Assault : int 236 263 294 190 276 204 110 238 335 211 ...
## $ UrbanPop: int 58 48 80 50 91 78 77 72 80 60 ...
## $ Rape : num 21.2 44.5 31 19.5 40.6 38.7 11.1 15.8 31.9 25.8 ...Como la magnitud de los valores difiere notablemente entre variables, se procede a escalarlas antes de aplicar el clustering.
Una forma sencilla de estimar el número K óptimo de clusters cuando no se dispone de información adicional en la que basarse, es aplicar el algoritmo de K-means para un rango de valores de K e identificar aquel valor a partir del cual la reducción en la suma total de varianza intra-cluster deja de ser sustancial. A esta estrategia se la conoce como método del codo o elbow method (en los siguientes apartados se detallan otras opciones).
La función fviz_nbclust() automatiza este proceso y genera una representación de los resultados.
library(factoextra)
fviz_nbclust(x = datos, FUNcluster = kmeans, method = "wss", k.max = 15,
diss = get_dist(datos, method = "euclidean"), nstart = 50)
Este mismo análisis también puede realizarse sin recurrir a la función fviz_nbclust().
calcular_totwithinss <- function(n_clusters, datos, iter.max=1000, nstart=50){
# Esta función aplica el algoritmo kmeans y devuelve la suma total de
# cuadrados internos.
cluster_kmeans <- kmeans(centers = n_clusters, x = datos, iter.max = iter.max,
nstart = nstart)
return(cluster_kmeans$tot.withinss)
}
# Se aplica esta función con para diferentes valores de k
total_withinss <- map_dbl(.x = 1:15,
.f = calcular_totwithinss,
datos = datos)
total_withinss## [1] 196.00000 102.86240 78.32327 56.40317 48.94420 42.83303 38.25764
## [8] 33.85843 29.86789 26.18348 24.05222 21.47090 20.15762 18.04643
## [15] 16.81152data.frame(n_clusters = 1:15, suma_cuadrados_internos = total_withinss) %>%
ggplot(aes(x = n_clusters, y = suma_cuadrados_internos)) +
geom_line() +
geom_point() +
scale_x_continuous(breaks = 1:15) +
labs(title = "Evolución de la suma total de cuadrados intra-cluster") +
theme_bw()
En este análisis, a partir de 4 clusters la reducción en la suma total de cuadrados internos parece estabilizarse, indicando que K = 4 es una buena opción.
El paquete factoextra también permite obtener visualizaciones de las agrupaciones resultantes. Si el número de variables (dimensionalidad) es mayor de 2, automáticamente realiza un PCA y representa las dos primeras componentes principales.
set.seed(123)
km_clusters <- kmeans(x = datos, centers = 4, 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 = datos, show.clust.cent = TRUE,
ellipse.type = "euclid", star.plot = TRUE, repel = TRUE) +
labs(title = "Resultados clustering K-means") +
theme_bw() +
theme(legend.position = "none")
Otra forma de representar los resultados obtenidos con un clustering k-means es mediante una network que muestre las observaciones conectadas por* clusters*. Es importante tener en cuenta que, a diferencia de la proyección PCA, en esta representación las distancias son arbitrarias, solo importan las conexiones.
set.seed(123)
km_clusters <- kmeans(x = datos, centers = 4, nstart = 50)
resultados <- data.frame(ciudad = names(km_clusters$cluster),
cluster = as.factor(km_clusters$cluster)) %>%
arrange(cluster)
library(igraph)
library(tidygraph)
library(ggraph)
datos_graph <- graph_from_data_frame(d = resultados, directed = TRUE)
datos_graph <- as_tbl_graph(datos_graph)
# Se añade información sobre a que cluster pertenece cada observacion
datos_graph <- datos_graph %>%
activate(nodes) %>%
left_join(resultados, by = c("name" = "ciudad"))
ggraph(graph = datos_graph) +
geom_edge_link(alpha = 0.5) +
geom_node_point(aes(color = cluster)) +
geom_node_text(aes(label = name), repel = TRUE, alpha = 0.5, size = 3) +
labs(title = "Resultados clustering K-means") +
theme_graph()
K-medoids clustering (PAM)
Idea intuitiva
K-medoids es un método de clustering muy similar a K-means en cuanto a que ambos agrupan las observaciones en K clusters, donde K es un valor preestablecido por el analista. La diferencia es que, en K-medoids, cada cluster está representado por una observación presente en el cluster (medoid), mientras que en K-means cada cluster está representado por su centroide, que se corresponde con el promedio de todas las observaciones del cluster pero con ninguna en particular.
Una definición más exacta del término medoid es: elemento dentro de un cluster cuya distancia (diferencia) promedio entre él y todos los demás elementos del mismo cluster es lo menor posible. Se corresponde con el elemento más central del cluster y por lo tanto puede considerarse como el más representativo. El hecho de utilizar medoids en lugar de centroides hace de K-medoids un método más robusto que K-means, viéndose menos afectado por outliers o ruido. A modo de idea intuitiva puede considerarse como la analogía entre media y mediana.
El algoritmo más empleado para aplicar K-medoids se conoce como PAM (Partitioning Around Medoids) y sigue los siguientes pasos:
Seleccionar K observaciones aleatorias como medoids iniciales. También es posible identificarlas de forma específica.
Calcular la matriz de distancia entre todas las observaciones si esta no se ha calculado anteriormente.
Asignar cada observación a su medoid más cercano.
Para cada uno de los clusters creados, comprobar si seleccionando otra observación como medoid se consigue reducir la distancia promedio del cluster, si esto ocurre, seleccionar la observación que consigue una mayor reducción como nuevo medoid.
Si al menos un medoid ha cambiado en el paso 4, volver al paso 3, de lo contrario, se termina el proceso.
A diferencia del algoritmo K-means, en el que se minimiza la suma total de cuadrados intra-cluster (suma de las distancias al cuadrado de cada observación respecto a su centroide), el algoritmo PAM minimiza la suma de las diferencias de cada observación respecto a su medoid.
Por lo general, el método de K-medoids se utiliza cuando se conoce o se sospecha de la presencia de outliers. Si esto ocurre, es recomendable utilizar como medida de similitud la distancia de Manhattan, ya que es menos sensible a outliers que la euclídea.
Ventajas y desventajas
K-medoids es un método de clustering más robusto que K-means, por lo es más adecuado cuando el set de datos contiene outliers o ruido.
Al igual que K-means, necesita que se especifique de antemano el número de clusters que se van a crear. Esto puede ser complicado de determinar si no se dispone de información adicional sobre los datos. Muchas de las estrategias empleadas en K-means para identificar el numero óptimo, puden aplicarse en K-medoids.
Para sets de datos grandes necesita muchos recursos computacionales. En tal situación se recomienda aplicar el método CLARA.
Ejemplo
El set de datos USArrests contiene información sobre el número de delitos (asaltos, asesinatos y secuestros) junto con el porcentaje de población urbana para cada uno de los 50 estados de USA. Se pretende estudiar si existe una agrupación subyacente de los estados mediante clustering. Dado que se sospecha de la presencia de outliers se recurre a K-medoids.
El proceso a seguir en R para aplicar el método de K-medoids es igual al seguido en K-means (ejemplo 2), pero en este caso empleando la función pam() del paquete cluster.
## 'data.frame': 50 obs. of 4 variables:
## $ Murder : num 13.2 10 8.1 8.8 9 7.9 3.3 5.9 15.4 17.4 ...
## $ Assault : int 236 263 294 190 276 204 110 238 335 211 ...
## $ UrbanPop: int 58 48 80 50 91 78 77 72 80 60 ...
## $ Rape : num 21.2 44.5 31 19.5 40.6 38.7 11.1 15.8 31.9 25.8 ...Como la magnitud de los valores difiere notablemente entre variables, se procede a escalarlas antes de aplicar el clustering.
Se evalúa la reducción de varianza total intra-cluster para un rango de valores K con el objetivo de identificar el número óptimo de clusters (elbow method). En este caso, dado que se sospecha de la presencia de outliers, se emplea la distancia de Manhattan como medida de similitud.
library(cluster)
library(factoextra)
fviz_nbclust(x = datos, FUNcluster = pam, method = "wss", k.max = 15,
diss = dist(datos, method = "manhattan"))
# Mismo análisis pero sin recurrir a factoextra
# ==============================================================================
calcular_suma_dif_interna <- function(n_clusters, datos, distancia = "manhattan"){
# Esta función aplica el algoritmo pam y devuelve la suma total de las
# diferencias internas
cluster_pam <- cluster::pam(x = datos, k = n_clusters, metric = distancia)
# El objeto pam almacena la suma de las diferencias respecto a los medoides en
# $objective["swap"]
return(cluster_pam$objective["swap"])
}
# Se aplica esta función con para diferentes valores de k
suma_dif_interna <- map_dbl(.x = 1:15,
.f = calcular_suma_dif_interna,
datos = datos)
data.frame(n_clusters = 1:15, suma_dif_interna = suma_dif_interna) %>%
ggplot(aes(x = n_clusters, y = suma_dif_interna)) +
geom_line() +
geom_point() +
scale_x_continuous(breaks = 1:15) +
labs(title = "Evolución de la suma total de diferencias intra-cluster") +
theme_bw()
Al igual que ocurría al aplicar K-means a estos datos, a partir de 4 clusters la reducción en la suma total de diferencias internas parece estabilizarse, indicando que K = 4 es una buena opción.
## Medoids:
## ID Murder Assault UrbanPop Rape
## Alabama 1 1.2425641 0.7828393 -0.5209066 -0.003416473
## Michigan 22 0.9900104 1.0108275 0.5844655 1.480613993
## Oklahoma 36 -0.2727580 -0.2371077 0.1699510 -0.131534211
## Iowa 15 -1.2829727 -1.3770485 -0.5899924 -1.060387812
## Clustering vector:
## Alabama Alaska Arizona Arkansas California
## 1 2 2 3 2
## Colorado Connecticut Delaware Florida Georgia
## 2 4 3 2 1
## Hawaii Idaho Illinois Indiana Iowa
## 3 4 2 3 4
## Kansas Kentucky Louisiana Maine Maryland
## 3 3 1 4 2
## Massachusetts Michigan Minnesota Mississippi Missouri
## 3 2 4 1 3
## Montana Nebraska Nevada New Hampshire New Jersey
## 3 3 2 4 3
## New Mexico New York North Carolina North Dakota Ohio
## 2 2 1 4 3
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 3 3 3 3 1
## South Dakota Tennessee Texas Utah Vermont
## 4 1 2 3 4
## Virginia Washington West Virginia Wisconsin Wyoming
## 3 3 4 4 3
## Objective function:
## build swap
## 1.730682 1.712075
##
## Available components:
## [1] "medoids" "id.med" "clustering" "objective" "isolation"
## [6] "clusinfo" "silinfo" "diss" "call" "data"El objeto devuelto por pam() contiene entre otra información: las observaciones que finalmente se han seleccionado como medoids ($medoids) y el cluster al que se ha asignado cada observación ($clustering).
fviz_cluster(object = pam_clusters, data = datos, ellipse.type = "t",
repel = TRUE) +
theme_bw() +
labs(title = "Resultados clustering PAM") +
theme(legend.position = "none")
# Como en k-medoids no hay centroides, no se muestran en la representación ni
# tampoco las distancias desde este al resto de observacionesLa función fviz_cluster() no permite resaltar las observaciones que actúan como medoids, sin embargo, al tratarse de un objeto ggplot2, es sencillo conseguirlo.
# Como hay más de 2 variables, se están representando las 2 primeras componentes
# de un PCA. Se tienen que calcular el PCA y extraer las proyecciones almacenadas
# en el elemento x.
medoids <- prcomp(datos)$x
# Se seleccionan únicamente las proyecciones de las observaciones que son medoids
medoids <- medoids[rownames(pam_clusters$medoids), c("PC1", "PC2")]
medoids <- as.data.frame(medoids)
# Se emplean los mismos nombres que en el objeto ggplot
colnames(medoids) <- c("x", "y")
# Creación del gráfico
fviz_cluster(object = pam_clusters, data = datos, 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")
CLARA
Idea intuitiva
Una de las limitaciones del método K-medoids-clustering es que su algoritmo requiere mucha memoria RAM, lo que impide que se pueda aplicar cuando el set de datos contiene varios miles de observaciones. CLARA (Clustering Large Applications) es un método que combina la idea de K-medoids con el resampling para que pueda aplicarse a grandes volúmenes de datos.
En lugar de intentar encontrar los medoids empleando todos los datos a la vez, CLARA selecciona una muestra aleatoria de un tamaño determinado y le aplica el algoritmo de PAM (K-medoids) para encontrar los clusters óptimos acorde a esa muestra. Utilizando esos medoids se agrupan las observaciones de todo el set de datos. La calidad de los medoids resultantes se cuantifica con la suma total de las distancias entre cada observación del set de datos y su correspondiente medoid (suma total de distancias intra-clusters). CLARA repite este proceso un número predeterminado de veces con el objetivo de reducir el bias de muestreo. Por último, se seleccionan como clusters finales los obtenidos con aquellos medoids que han conseguido menor suma total de distancias. A continuación, se describen los pasos del algoritmo CLARA.
Se divide aleatoriamente el set de datos en n partes de igual tamaño, donde n es un valor que determina el analista.
Para cada una de las n partes:
2.1 Aplicar el algoritmo PAM e identificar cuáles son los k medoids.
2.2 Utilizando los medoids del paso anterior agrupar todas las observaciones del set de datos.
2.3 Calcular la suma total de las distancias entre cada observación del set de datos y su correspondiente medoid (suma total de distancias intra-clusters).
Seleccionar como clustering final aquel que ha conseguido menor suma total de distancias intra-clusters en el paso 2.3.
Ejemplo
Para ilustrar la aplicación del método CLARA se simula un set de datos bidimensional (dos variables) con 500 observaciones, de las cuales 200 pertenecen a un grupo y 300 a otro (número de grupos reales = 2).
set.seed(1234)
grupo_1 <- cbind(rnorm(n = 200, mean = 0, sd = 8),
rnorm(n = 200, mean = 0, sd = 8))
grupo_2 <- cbind(rnorm(n = 300, mean = 30, sd = 8),
rnorm(n = 300, mean = 30, sd = 8))
datos <- rbind(grupo_1, grupo_2)
colnames(datos) <- c("x", "y")
head(datos)## x y
## [1,] -9.656526 3.881815
## [2,] 2.219434 5.574150
## [3,] 8.675529 1.484111
## [4,] -18.765582 5.605868
## [5,] 3.432998 2.493448
## [6,] 4.048447 6.083699La función clara() del paquete cluster permite aplicar el algoritmo CLARA. Entre sus argumentos destaca: una matriz numérica x donde cada fila es una observación, el número de clusters k, la medida de distancia empleada metric (euclídea o manhattan), si los datos se tienen que estandarizar stand, el número de partes samples en las que se divide el set de datos (recomendable 50) y si se utiliza el algoritmo PAM pamLike.
library(cluster)
library(factoextra)
clara_clusters <- clara(x = datos, k = 2, metric = "manhattan", stand = TRUE,
samples = 50, pamLike = TRUE)
clara_clusters## Call: clara(x = datos, k = 2, metric = "manhattan", stand = TRUE, samples = 50, pamLike = TRUE)
## Medoids:
## x y
## [1,] -0.2780831 1.269004
## [2,] 30.3298450 29.233511
## Objective function: 0.8629488
## Clustering vector: int [1:500] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
## Cluster sizes: 202 298
## Best sample:
## [1] 17 30 33 34 63 71 80 84 85 100 115 141 149 162 165 167 168 171 188
## [20] 214 223 249 259 261 291 302 318 320 333 386 391 402 412 417 419 422 437 439
## [39] 440 450 469 472 498 499
##
## Available components:
## [1] "sample" "medoids" "i.med" "clustering" "objective"
## [6] "clusinfo" "diss" "call" "silinfo" "data"fviz_cluster(object = clara_clusters, ellipse.type = "t", geom = "point",
pointsize = 2.5) +
theme_bw() +
labs(title = "Resultados clustering CLARA") +
theme(legend.position = "none")
Hierarchical clustering
Hierarchical clustering es una alternativa a los métodos de partitioning clustering que no requiere que se pre-especifique el número de clusters. Los métodos que engloba el hierarchical clustering se subdividen en dos tipos dependiendo de la estrategia seguida para crear los grupos:
Agglomerative clustering (bottom-up): el agrupamiento se inicia en la base del árbol, donde cada observación forma un cluster individual. Los clusters se van combinado a medida que la estructura crece hasta converger en una única “rama” central.
Divisive clustering (top-down): es la estrategia opuesta al agglomerative clustering, se inicia con todas las observaciones contenidas en un mismo cluster y se suceden divisiones hasta que cada observación forma un cluster individual.
En ambos casos, los resultados pueden representarse de forma muy intuitiva en una estructura de árbol llamada dendrograma.
Agglomerative
La estructura resultante de un agglomerative hierarchical clustering se obtiene mediante un algoritmo sencillo.
El proceso se inicia considerando cada una de las observaciones como un cluster individual, formando así la base del dendrograma (hojas).
Se inicia un proceso iterativo hasta que todas las observaciones pertenecen a un único cluster:
2.1 Se calcula la distancia entre cada posible par de los n clusters. El investigador debe determinar el tipo de medida emplea para cuantificar la similitud entre observaciones o grupos (distancia y linkage).
2.2 Los dos clusters más similares se fusionan, de forma que quedan n-1 clusters.
Determinar dónde cortar la estructura de árbol generada (dendrograma).
La siguiente imagen muestra cómo se van sucediendo las agrupaciones a medida que avanzan las primeras iteraciones del algoritmo.
Imagen obtenida del libro ISLR
Para que el proceso de agrupamiento pueda llevarse a cabo tal como indica el algoritmo anterior, es necesario definir cómo se cuantifica la similitud entre dos clusters. Es decir, se tiene que extender el concepto de distancia entre pares de observaciones para que sea aplicable a pares de grupos, cada uno formado por varias observaciones. A este proceso se le conoce como linkage. A continuación, se describen los 5 tipos de linkage más empleados y sus definiciones.
Complete or Maximum: Se calcula la distancia entre todos los posibles pares formados por una observación del cluster A y una del cluster B. La mayor de todas ellas se selecciona como la distancia entre los dos clusters. Se trata de la medida más conservadora (maximal intercluster dissimilarity).
Single or Minimum: Se calcula la distancia entre todos los posibles pares formados por una observación del cluster A y una del cluster B. La menor de todas ellas se selecciona como la distancia entre los dos clusters. Se trata de la medida menos conservadora (minimal intercluster dissimilarity).
Average: Se calcula la distancia entre todos los posibles pares formados por una observación del cluster A y una del cluster B. El valor promedio de todas ellas se selecciona como la distancia entre los dos clusters (mean intercluster dissimilarity).
Centroid: Se calcula el centroide de cada uno de los clusters y se selecciona la distancia entre ellos como la distancia entre los dos clusters.
Ward: Se trata de un método general. La selección del par de clusters que se combinan en cada paso del agglomerative hierarchical clustering se basa en el valor óptimo de una función objetivo, pudiendo ser esta última cualquier función definida por el analista. El conocido método Ward’s minimum variance es un caso particular en el que el objetivo es minimizar la suma total de varianza intra-cluster. En cada paso, se identifican aquellos 2 clusters cuya fusión conlleva menor incremento de la varianza total intra-cluster. La implementación en
Rde este método ha sido causa de confusiones (Ward’s Hierarchical Agglomerative Clustering Method: Which Algorithms Implement Ward’s Criterion? by Fionn Murtagh y Pierre Legendre). Si se emplea la funciónhclust()se tiene que especificarmethod = "ward.D2", mientras que en la funciónagnes()esmethod = ward".
Los métodos de linkage complete, average y Ward’s minimum variance suelen ser los preferidos por los analistas debido a que generan dendrogramas más compensados. Sin embargo, no se puede determinar que uno sea mejor que otro, ya que depende del caso de estudio en cuestión. Por ejemplo, en genómica, se emplea con frecuencia el método de centroides. Junto con los resultados de un proceso de hierarchical clustering siempre hay que indicar qué distancia se ha empleado, así como el tipo de linkage, ya que, dependiendo de estos, los resultados pueden variar en gran medida.
Divisive
El algoritmo más conocido de divisive hierarchical clustering es DIANA (DIvisive ANAlysis Clustering). Este algoritmo se inicia con un único cluster que contiene todas las observaciones, a continuación, se van sucediendo divisiones hasta que cada observación forma un cluster independiente. En cada iteración, se selecciona el cluster con mayor diametro, entendiendo por diámetro de un cluster la mayor de las diferencias entre dos de sus observaciones. Una vez seleccionado el cluster, se identifica la observación más dispar, que es aquella con mayor distancia promedio respecto al resto de observaciones que forman el cluster, esta observación inicia el nuevo cluster. Se reasignan las observaciones en función de si están más próximas al nuevo cluster o al resto de la partición, dividiendo así el cluster seleccionado en dos nuevos clusters.
Todas las n observaciones forman un único cluster.
Repetir hasta que hayan n clusters:
2.1 Calcular para cada cluster la mayor de las distancias entre pares de observaciones (diámetro del cluster).
2.2 Seleccionar el cluster con mayor diámetro.
2.3 Calcular la distancia media de cada observación respecto a las demas.
2.4 La observación más distante inicia un nuevo cluster
2.5 Se reasignan las observaciones restantes al nuevo cluster o al viejo dependiendo de cual está más próximo.
A diferencia del clustering aglomerativo, en el que hay que elegir un tipo de distancia y un método de linkage, en el clustering divisivo solo hay que elegir la distancia, no hay linkage.
Dendrograma
Para ilustrar cómo se interpreta un dendograma, se simula un set de datos y se somete a un proceso de hierarchical clustering.
Supóngase que se dispone de 45 observaciones en un espacio de dos dimensiones, que pertenecen a 3 grupos. Aunque se ha coloreado de forma distinta cada uno de los grupos para facilitar la comprensión de la idea, se va a suponer que se desconoce esta información y que se desea aplicar el método de hierarchical clustering para intentar reconocer los grupos.
Al aplicar hierarchical clustering, empleando como medida de similitud la distancia euclídea y linkage complete, se obtiene el siguiente dendrograma. Como los datos se han simulado en aproximamdamente la misma escala, no es necesario estandarizarlos, de no ser así, sí se tendrían que estandarizar.
En la base del dendrograma, cada observación forma una terminación individual conocida como hoja o leaf del árbol. A medida que se asciende por la estructura, pares de hojas se fusionan formando las primeras ramas. Estas uniones (nodos) se corresponden con los pares de observaciones más similares. También ocurre que ramas se fusionan con otras ramas o con hojas. Cuanto más temprana (más próxima a la base del dendrograma) ocurre una fusión, mayor es la similitud. Esto significa que, para cualquier par de observaciones, se puede identificar el punto del árbol en el que las ramas que contienen dichas observaciones se fusionan. La altura a la que esto ocurre (eje vertical) indica cómo de similares/diferentes son las dos observaciones. Los dendrogramas, por lo tanto, se deben interpretar únicamente en base al eje vertical y no por las posiciones que ocupan las observaciones en el eje horizontal, esto último es simplemente por estética y puede variar de un programa a otro. Por ejemplo, la observación 10 es la más similar a la 8 ya que es la primera fusión que recibe la observación 8 (y viceversa). Podría resultar tentador decir que la observación 26, situada inmediatamente a la izquierda de la 8, es la siguiente más similar, sin embargo, las observaciones 28 y 37 son más similares a la 8 a pesar de que se encuentran más alejadas en el eje horizontal. Del mismo modo, no es correcto decir que la observación 28 es más similar a la observación 8 de lo que lo es la 37 por el hecho de que está más próxima en el eje horizontal. Prestando atención a la altura en que las respectivas ramas se unen, la única conclusión válida es que la similitud entre los pares 8-28 y 8-37 es la misma.
Verificar el árbol resultante
Una vez creado el dendrograma, hay que evaluar hasta qué punto su estructura refleja las distancias originales entre observaciones. Una forma de hacerlo es empleando el coeficiente de correlación entre las distancias cophenetic del dendrograma (altura de los nodos) y la matriz de distancias original. Cuanto más cercano es el valor a 1, mejor refleja el dendrograma la verdadera similitud entre las observaciones. Valores superiores a 0.75 suelen considerarse como buenos. Esta medida puede emplearse como criterio de ayuda para escoger entre los distintos métodos de linkage. En R, la función cophenetic() calcula las distancias cophenetic de un hierarchical clustering.
data(USArrests)
datos <- scale(USArrests)
# Matriz de distancias euclídeas
mat_dist <- dist(x = datos, 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")
cor(x = mat_dist, cophenetic(hc_euclidea_complete))## [1] 0.6979437## [1] 0.7180382Para estos datos, el método de linkage average consigue representar ligeramente mejor la similitud entre observaciones.
Cortar el árbol para generar los clusters
Además de representar en un dendrograma la similitud entre observaciones, se tiene que poder identificar el número de clusters creados y qué observaciones forman parte de cada uno. Si se realiza un corte horizontal a una determinada altura del dendrograma, el número de ramas que sobrepasan (en sentido ascendente) dicho corte se corresponde con el número de clusters. La siguiente imagen muestra dos veces el mismo dendrograma. Si se realiza el corte a la altura de 5, se obtienen dos clusters, mientras que si se hace a la de 3.5 se obtienen 4. La altura de corte tiene por lo tanto la misma función que el valor K en K-means-clustering: controla el número de clusters obtenidos.
library(factoextra)
datos <- USArrests
datos <- scale(datos)
set.seed(101)
hc_euclidea_completo <- hclust(d = dist(x = datos, method = "euclidean"),
method = "complete")
fviz_dend(x = hc_euclidea_completo, k = 2, cex = 0.6) +
geom_hline(yintercept = 5.5, linetype = "dashed") +
labs(title = "Herarchical clustering",
subtitle = "Distancia euclídea, Lincage complete, K=2")
fviz_dend(x = hc_euclidea_completo, k = 4, cex = 0.6) +
geom_hline(yintercept = 3.5, linetype = "dashed") +
labs(title = "Herarchical clustering",
subtitle = "Distancia euclídea, Lincage complete, K=4")
Dos propiedades adicionales se derivan de la forma en que se generan los clusters en el método de hierarchical clustering:
Dada la longitud variable de las ramas, siempre existe un intervalo de altura para el que cualquier corte da lugar al mismo número de clusters. En el ejemplo anterior, todos los cortes entre las alturas 5 y 6 tienen como resultado los mismos 2 clusters.
Con un solo dendrograma se dispone de la flexibilidad para generar cualquier número de clusters desde 1 a n. La selección del número óptimo puede valorarse de forma visual, tratando de identificar las ramas principales en base a la altura a la que ocurren las uniones. En el ejemplo expuesto es razonable elegir entre 2 o 4 clusters.
Una forma menos frecuente de representar los resultados de un hierarchical clustering es combinándolos con una reducción de dimensionalidad por PCA. Primero, se calculan las componentes principales y se representan las observaciones en un scatterplot empleando las dos primeras componentes, finalmente se colorean los clusters mediante elipses.
fviz_cluster(object = list(data=datos, cluster=cutree(hc_euclidea_completo, k=4)),
ellipse.type = "convex", repel = TRUE, show.clust.cent = FALSE,
labelsize = 8) +
labs(title = "Hierarchical clustering + Proyección PCA",
subtitle = "Distancia euclídea, Lincage complete, K=4") +
theme_bw() +
theme(legend.position = "bottom")
Ejemplo
Los siguientes datos simulados contienen observaciones que pertenecen a cuatro grupos distintos. Se pretende aplicar hierarchical clustering aglomerativo con el fin de identificarlos.
set.seed(101)
# Se simulan datos aleatorios con dos dimensiones
datos <- matrix(rnorm(n = 100*2), nrow = 100, ncol = 2,
dimnames = list(NULL,c("x", "y")))
datos <- as.data.frame(datos)
# Se determina la media que va a tener cada grupo en cada una de las dos
# dimensiones. En total 2*4 medias. Este valor se utiliza para separar
# cada grupo de los demás.
media_grupos <- matrix(rnorm(n = 8, mean = 0, sd = 4), nrow = 4, ncol = 2,
dimnames = list(NULL, c("media_x", "media_y")))
media_grupos <- as.data.frame(media_grupos)
media_grupos <- media_grupos %>% mutate(grupo = c("a","b","c","d"))
# Se genera un vector que asigne aleatoriamente cada observación a uno de
# los 4 grupos
datos <- datos %>% mutate(grupo = sample(x = c("a","b","c","d"),
size = 100,
replace = TRUE))
# Se incrementa el valor de cada observación con la media correspondiente al
# grupo asignado.
datos <- left_join(datos, media_grupos, by = "grupo")
datos <- datos %>% mutate(x = x + media_x,
y = y + media_y)
datos <- datos %>% select(grupo, x, y)
ggplot(data = datos, aes(x = x, y = y, color = grupo)) +
geom_point(size = 2.5) +
theme_bw()
Al aplicar un hierarchical clustering aglomerativo se tiene que escoger una medida de distancia (1-similitud) y un tipo de linkage. En este caso, se emplea la función hclust(), a la que se pasa como argumento una matriz de distancia euclidea y el tipo de linkages. Se comparan los resultados con los linkages complete, single y average. Dado que los datos se han simulado considerando que las dos dimensiones tienen aproximadamente la misma magnitud, no es necesario escalarlos ni centrarlos.
# Se calculan las distancias
matriz_distancias <- dist(x = datos[, c("x", "y")], method = "euclidean")
set.seed(567)
hc_euclidea_completo <- hclust(d = matriz_distancias, method = "complete")
hc_euclidea_single <- hclust(d = matriz_distancias, method = "single")
hc_euclidea_average <- hclust(d = matriz_distancias, method = "average")Los objetos devueltos por hclust() pueden representarse en forma de dendrograma con la función plot() o con la función fviz_dend() del paquete factoextra.
par(mfrow = c(3,1))
plot(x = hc_euclidea_completo, cex = 0.6, xlab = "", ylab = "", sub = "",
main = "Distancia euclídea, Linkage complete")
plot(x = hc_euclidea_single, cex = 0.6, xlab = "", ylab = "", sub = "",
main = "Distancia euclídea, Linkage single")
plot(x = hc_euclidea_average, cex = 0.6, xlab = "", ylab = "", sub = "",
main = "Distancia euclídea, Linkage average")
El conocer que existen 4 grupos en la población permite evaluar qué linkage consigue los mejores resultados. En este caso, los tres tipos identifican claramente 4 clusters, si bien esto no significa que en los 3 dendrogramas los clusters estén formados por exactamente las mismas observaciones.
Una vez creado el dendrograma, se tiene que decidir a qué altura se corta para generar los clusters. La función cutree() recibe como input un dendrograma y devuelve el cluster al que se ha asignado cada observación dependiendo del número de clusters especificado (argumento k) o la altura de corte indicada (argumento h).
## [1] 1 2 3 1 2 3 1 2 3 3 4 3 3 4 2 4 1 4 1 1 1 1 1 4 3 3 2 1 2 3 1 4 1 2 2 4 4
## [38] 2 4 4 1 2 2 1 1 4 1 3 1 2 1 3 3 4 4 2 4 2 3 3 2 1 1 1 2 2 3 3 4 1 1 3 3 4
## [75] 3 4 1 3 3 3 2 3 2 3 1 4 3 1 1 3 3 2 1 2 3 3 4 4 3 3## [1] 1 2 3 1 2 3 1 2 3 3 4 3 3 4 2 4 1 4 1 1 1 1 1 4 3 3 2 1 2 3 1 4 1 2 2 4 4
## [38] 2 4 4 1 2 2 1 1 4 1 3 1 2 1 3 3 4 4 2 4 2 3 3 2 1 1 1 2 2 3 3 4 1 1 3 3 4
## [75] 3 4 1 3 3 3 2 3 2 3 1 4 3 1 1 3 3 2 1 2 3 3 4 4 3 3Una forma visual de comprobar los errores en las asignaciones es indicando en el argumento labels el grupo real al que pertenece cada observación. Si la agrupación resultante coincide con los grupos reales, entonces, dentro de cada clusters las labels serán las mismas.
plot(x = hc_euclidea_completo, cex = 0.6, sub = "",
main = "Distancia euclídea, Linkage complete, k=4",
xlab = "", ylab = "", labels = datos[, "grupo"])
abline(h = 6, lty = 2)
##
## a b c d
## 1 0 28 0 0
## 2 1 0 0 20
## 3 31 0 0 0
## 4 0 0 20 0El método de hierarchical clustering aglomerativo con linkage completo y k=4 ha sido capaz de agrupar correctamente todas las observaciones.
Ejemplo divisivo
Con los mismos datos que en el ejemplo anterior, se emplea la función diana() de paquete cluster para crear un heirarchical clustering divisivo con diatancia euclídea.
matriz_distancias <- dist(x = datos[, c("x", "y")], method = "euclidean")
hc_diana <- diana(x = matriz_distancias, diss = TRUE, stand = FALSE)
fviz_dend(x = hc_diana, cex = 0.5) +
labs(title = "Hierarchical clustering divisivo",
subtitle = "Distancia euclídea")
Ejemplo Clasificar tumores por su perfil genético
El set de datos NCI60 contiene información genética de 64 líneas celulares cancerígenas. Para cada una de ellas, se ha cuantificado la expresión de 6830 genes mediante tecnología microarray. Los investigadores conocen el tipo de cáncer (histopatología) al que pertenece cada línea celular y quieren utilizar esta información para evaluar si el método de clustering (agglomerative hierarchical clustering) es capaz de agrupar correctamente las líneas empleando los niveles de expresión génica.
Los métodos de clustering son unsupervised, lo que significa que al aplicarlos no se hace uso de la variable respuesta, en este caso el tipo de cáncer. Una vez obtenidos los resultados se añade esta información para determinar si es posible agrupar a las líneas celulares empleando su perfil de expresión.
## List of 2
## $ data: num [1:64, 1:6830] 0.3 0.68 0.94 0.28 0.485 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:64] "V1" "V2" "V3" "V4" ...
## .. ..$ : chr [1:6830] "1" "2" "3" "4" ...
## $ labs: chr [1:64] "CNS" "CNS" "CNS" "RENAL" ...# Los datos están almacenados en forma de lista, un elemento contiene los niveles
# de expresión y otro el tipo de cáncer
expresion <- NCI60$data
tipo_cancer <- NCI60$labsUna exploración inicial de los datos permite conocer el número de líneas celulares que hay de cada tipo de cáncer.
## tipo_cancer
## BREAST CNS COLON K562A-repro K562B-repro LEUKEMIA
## 7 5 7 1 1 6
## MCF7A-repro MCF7D-repro MELANOMA NSCLC OVARIAN PROSTATE
## 1 1 8 9 6 2
## RENAL UNKNOWN
## 9 1El siguiente paso antes de aplicar el método de clustering es decidir cómo se va a cuantificar la similitud. Los perfiles de expresión génica son un claro ejemplo de que es necesario comprender el problema en cuestión para hacer la elección adecuada. Para ilustrarlo, se simula el perfil de 16 genes en 4 tumores (1 de referencia contra el que se comparan los otros 3).
Debajo de cada gráfico se indica la distancia euclídea entre cada par de tumores. Acorde a esta medida, el tumor menos parecido al de referencia es el A y el más parecido el B. Sin embargo, analizando los perfiles con detenimiento, puede observarse que el tumor A tiene exactamente el mismo perfil que el tumor de referencia, pero desplazado unas unidades; mientras que el tumor B tiene un perfil totalmente distinto. Para evitar que las diferencias en magnitud determinen la similitud, se normalizan los valores de expresión de forma que tengan media 0 y desviación estándar 1.
Una vez aplicada la estandarización, las distancias obtenidas en las comparaciones tienen más sentido dentro del contexto de los perfiles de expresión génica. La similitud entre el tumor A y el de referencia es total (distancia 0), y el tumor B pasa a ser el más distinto. El mismo resultado se hubiese obtenido empleando como medida de similitud la correlación lineal de Pearson.
Se procede a escalar la matriz de expresión NCI60 y se aplica hierarchical clustering con linkage completo, average y single.
expresion <- scale(expresion, center = TRUE, scale = TRUE)
matriz_distancias <- dist(x = expresion, method = "euclidean")
hc_completo <- hclust(d = matriz_distancias, method = "complete")
hc_average <- hclust(d = matriz_distancias, method = "average")
hc_single <- hclust(d = matriz_distancias, method = "single")
par(mfrow = c(3, 1))
plot(hc_completo, labels = tipo_cancer, ylab = "", xlab = "", sub = "",
main = "Linkage completo", cex = 0.8)
plot(hc_average, labels = tipo_cancer, ylab = "", xlab = "", sub = "",
main = "Linkage average", cex = 0.8)
plot(hc_single, labels = tipo_cancer, ylab = "", xlab = "", sub = "",
main = "Linkage single", cex = 0.8)
La elección del tipo de linkage influye de forma notable en los dendrogramas obtenidos. Por lo general, single linkage tiende a formar clusters muy grandes en los que las observaciones individuales se unen una a una. El completo y average suele generar dendrogramas más compensados con clusters más definidos, tal como ocurre en este ejemplo.
A pesar de que la agrupación no es perfecta, los clusters tienden a segregar bastante bien las líneas celulares procedentes de leukemia, melanoma y renal. Véase los aciertos cuando el dendrograma se corta a una altura tal que genera 4 clusters.
clusters <- cutree(tree = hc_completo, k = 4)
table(clusters, tipo_cancer, dnn = list("clusters", "tipo de cáncer"))## tipo de cáncer
## clusters BREAST CNS COLON K562A-repro K562B-repro LEUKEMIA MCF7A-repro
## 1 2 3 2 0 0 0 0
## 2 3 2 0 0 0 0 0
## 3 0 0 0 1 1 6 0
## 4 2 0 5 0 0 0 1
## tipo de cáncer
## clusters MCF7D-repro MELANOMA NSCLC OVARIAN PROSTATE RENAL UNKNOWN
## 1 0 8 8 6 2 8 1
## 2 0 0 1 0 0 1 0
## 3 0 0 0 0 0 0 0
## 4 1 0 0 0 0 0 0Todas las líneas celulares de leukemia caen en el cluster 3, todas las de melanoma y ovarian en el 1. Las líneas celulares de breast son las más distribuidas (heterogéneas en su perfil genético), ya que están presentes en los clusters 1, 2 y 4.
Nota: de los 6830 genes medidos, muchos pueden estar aportando información redundante o puede que no varíen lo suficiente como para contribuir al modelo. Con el fin de eliminar todo ese ruido y mejorar los resultados del clustering, es aconsejable filtrar aquellos genes cuya expresión no supere un mínimo de varianza. Del mismo modo, puede ser útil evaluar si aplicando un Principal Component Analysis se consigue capturar la mayor parte de la información en unas pocas componentes y utilizarlas como variables de clustering. Este es un problema muy importante en la aplicación del clustering a perfiles genéticos. En la siguiente sección se describe con más detalle cómo solucionarlo.
Hierarchical K-means clustering
K-means es uno de los métodos de clustering más utilizados y cuyos resultados son satisfactorios en muchos escenarios, sin embargo, como se ha explicado en apartados anteriores, sufre las limitaciones de necesitar que se especifique el número de clusters de antemano y de que sus resultados puedan variar en función de la iniciación aleatoria. Una forma de contrarrestar estos dos problemas es combinando el K-means con el hierarchical clustering. Los pasos a seguir son los siguientes:
Aplicar hierarchical clustering a los datos y cortar el árbol en k clusters. El número óptimo puede elegirse de forma visual o con cualquiera de los métodos explicados en la sección Número óptimo de clusters.
Calcular el centro (por ejemplo, la media) de cada cluster.
Aplicar k-means clustering empleando como centroides iniciales los centros calculados en el paso
El algoritmo de K-means tratará de mejorar la agrupación hecha por el hierarchical clustering en el paso 1, de ahí que las agrupaciones finales puedan variar respecto a las iniciales.
Ejemplo
El set de datos USArrests contiene información sobre el número de delitos (asaltos, asesinatos y secuestros) junto con el porcentaje de población urbana para cada uno de los 50 estados de USA. Se pretende estudiar si existe una agrupación subyacente de los estados empleando Hierarchical K-means clustering.
Como la magnitud de los valores difiere notablemente entre variables, se procede a escalarlas antes de aplicar el clustering.
La función hkmeans() del paquete factoextra permite aplicar el método hierarchical K-means clustering de forma muy similar a la función estándar kmeans().
library(factoextra)
# Se obtiene el dendrograma de hierarchical clustering para elegir el número de
# clusters.
set.seed(101)
hc_euclidea_completo <- hclust(d = dist(x = datos, method = "euclidean"),
method = "complete")
fviz_dend(x = hc_euclidea_completo, cex = 0.5, main = "Linkage completo",
sub = "Distancia euclídea") +
theme(plot.title = element_text(hjust = 0.5, size = 15))
Empleando la representación del dendrograma se considera que existen 4 grupos.
hkmeans_cluster <- hkmeans(x = datos, hc.metric = "euclidean",
hc.method = "complete", k = 4)
hkmeans_cluster## Hierarchical K-means clustering with 4 clusters of sizes 8, 13, 16, 13
##
## Cluster means:
## Murder Assault UrbanPop Rape
## 1 1.4118898 0.8743346 -0.8145211 0.01927104
## 2 0.6950701 1.0394414 0.7226370 1.27693964
## 3 -0.4894375 -0.3826001 0.5758298 -0.26165379
## 4 -0.9615407 -1.1066010 -0.9301069 -0.96676331
##
## Clustering vector:
## Alabama Alaska Arizona Arkansas California
## 1 2 2 1 2
## Colorado Connecticut Delaware Florida Georgia
## 2 3 3 2 1
## Hawaii Idaho Illinois Indiana Iowa
## 3 4 2 3 4
## Kansas Kentucky Louisiana Maine Maryland
## 3 4 1 4 2
## Massachusetts Michigan Minnesota Mississippi Missouri
## 3 2 4 1 2
## Montana Nebraska Nevada New Hampshire New Jersey
## 4 4 2 4 3
## New Mexico New York North Carolina North Dakota Ohio
## 2 2 1 4 3
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 3 3 3 3 1
## South Dakota Tennessee Texas Utah Vermont
## 4 1 2 3 4
## Virginia Washington West Virginia Wisconsin Wyoming
## 3 3 4 4 3
##
## Within cluster sum of squares by cluster:
## [1] 8.316061 19.922437 16.212213 11.952463
## (between_SS / total_SS = 71.2 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault" "data"
## [11] "hclust"fviz_cluster(object = hkmeans_cluster, pallete = "jco", repel = TRUE) +
theme_bw() + labs(title = "Hierarchical k-means Clustering")
Fuzzy clustering
Idea intuitiva
Los métodos de clustering descritos hasta ahora (K-means, hierarchical, K-medoids, CLARA…) asignan cada observación únicamente a un cluster, de ahí que también se conozcan como hard clustering. Los métodos de fuzzy clustering o soft clustering se caracterizan porque, cada observación, puede pertenecer potencialmente a varios clusters, en concreto, cada observación tiene asignado un grado de pertenencia a cada uno de los cluster.
Fuzzy c-means (FCM) es uno de los algoritmos más empleado para generar fuzzy clustering. Se asemeja en gran medida al algoritmo de k-means pero con dos diferencias:
El cálculo de los centroides de los clusters. La definición de centroide empleada por c-means es: la media de todas las observaciones del set de datos ponderada por la probabilidad de pertenecer a al cluster.
Devuelve para cada observación la probabilidad de pertenecer a cada cluster.
Ejemplo
El set de datos USArrests contiene información sobre el número de delitos (asaltos, asesinatos y secuestros) junto con el porcentaje de población urbana para cada uno de los 50 estados de USA. Se pretende estudiar si existe una agrupación subyacente de los estados empleando fuzzy clustering.
Como la magnitud de los valores difiere notablemente entre variables, se procede a escalarlas antes de aplicar el clustering.
La función fanny() (fuzzy analysis) del paquete cluster permite aplicar el algoritmo c-means clustering.
library(cluster)
fuzzy_cluster <- fanny(x = datos, diss = FALSE, k = 3, metric = "euclidean",
stand = FALSE)El objeto devuelto fanny() incluye entre sus elementos: una matriz con el grado de pertenencia de cada observación a cada cluster (las columnas son los clusters y las filas las observaciones).
## [,1] [,2] [,3]
## Alabama 0.4676004 0.3144516 0.2179480
## Alaska 0.4278809 0.3178707 0.2542484
## Arizona 0.5092197 0.2945668 0.1962135
## Arkansas 0.2934077 0.3787718 0.3278205
## California 0.4668527 0.3084149 0.2247324
## Colorado 0.4542018 0.3236683 0.2221299El coeficiente de partición Dunn normalizado y sin normalizar. Valores normalizados próximos a 0 indican que la estructura tiene un alto nivel fuzzy y valores próximos a 1 lo contrario.
## dunn_coeff normalized
## 0.37371071 0.06056606El cluster al que se ha asignado mayoritariamente cada observación.
## Alabama Alaska Arizona Arkansas California Colorado
## 1 1 1 2 1 1Para obtener una representación gráfica del clustering se puede emplear la función fviz_cluster().
library(factoextra)
fviz_cluster(object = fuzzy_cluster, repel = TRUE, ellipse.type = "norm",
pallete = "jco") + theme_bw() + labs(title = "Fuzzy Cluster plot")
Model based clustering
Idea intuitiva
El clustering basado en modelos considera que las observaciones proceden de una distribución que es a su vez una combinación de dos o más componentes (clusters), cada uno con una distribución propia. En principio, cada cluster puede estar descrito por cualquier función de densidad, pero normalmente se asume que siguen una distribución multivariante normal. Para estimar los parámetros que definen la función de distribución de cada cluster (media y matriz de covarianza si se asume que son de tipo normal) se recurre al algoritmo de Expectation-Maximization (EM). Este resuelve distintos modelos en los que el volumen, forma y orientación de las distribuciones pueden considerarse iguales para todos los clusters o distintas para cada uno. Por ejemplo, un posible modelo es: volumen constante, forma variable, orientación variable.
El paquete mclust emplea maximum likelihood para ajustar todos estos modelos con distinto número k de clusters y selecciona el mejor en base al Bayesian Information Criterion (BIC).
Ejemplo
El set de datos diabetes del paquete mclust contiene 3 parámetros sanguíneos medidos en 145 pacientes con 3 tipos distintos de diabetes. Se pretende emplear model-based-clustering para encontrar las agrupaciones.
# Estandarización de variables
datos <- scale(diabetes[, -1])
# Model-based-clustering
model_clustering <- Mclust(data = datos, G = 1:10)
summary(model_clustering)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 3
## components:
##
## log-likelihood n df BIC ICL
## -169.0908 145 29 -482.5069 -501.4662
##
## Clustering table:
## 1 2 3
## 81 36 28El algoritmo de ajuste selecciona como mejor modelo el formado por 3 clusters, cada uno con forma elipsoidal y con volume, shape y orientation propias.
El clustering basado en modelos es de tipo fuzzy, es decir, para cada observación se calcula un grado de pertenencia a cada cluster y se asigna finalmente al que mayor valor tiene.
## [,1] [,2] [,3]
## 1 0.9906745 0.008991332 3.341728e-04
## 2 0.9822128 0.017783229 3.974744e-06
## 3 0.9777871 0.022157665 5.527579e-05
## 4 0.9774763 0.022312280 2.113743e-04
## 5 0.9208978 0.079034264 6.789759e-05
## 6 0.9863472 0.012977950 6.748263e-04## 1 2 3 4 5 6
## 1 1 1 1 1 1La visualización del clustering puede hacerse mediante la función plot.Mclust() o mediante fviz_mclust().