logo

Introducción

Los métodos de particionamiento son una clase de técnicas de formación de grupos que producen una partición de las observaciones en un número específico de grupos, mediante la minimización o maximización de algún criterio numérico. Veamos a continuación los criterios más conocidos de métodos de particionamiento.

Método de k-medias

A pesar de ser uno de los métodos de particionamiento más antiguos, sigue siendo muy utilizado, puesto que es un algoritmo fácil de implementar, computacionalmente eficiente y necesita de poco espacio de almacenamiento. Este algoritmo de clasificación no supervisada está basado en la partición de un conjunto de \(n\) observaciones en \(k\) grupos en el que cada observación pertenece al grupo cuya distancia es menor.

El algoritmo k-medias (k-means) necesita definir antes de ejecutar el algoritmo un número de grupos o clústers. Cada grupo o clúster está definido por un punto, generalmente identificado como el centro y llamado semilla o centroide del clúster.

A grandes rasgos, el algoritmo funciona en dos fases principales, que corresponden a las dos siguientes tareas:

  1. En primer lugar, la fase de inicialización, identifica \(k\) puntos como centroides iniciales. Aunque no es necesario que sean puntos del conjunto de datos.

  2. La segunda fase es iterativa, y consiste en:

  • Asignar a cada centroide los puntos del conjunto de datos más próximos, formando \(k\) grupos disjuntos.
  • Recalcular los centroides en base a los puntos que forman parte de su grupo o partición.

Estos dos pasos se repiten hasta que se satisface alguna condición de parar.

De los pasos anteriores, se hace necesario definir algunos criterios.

1. ¿Cómo podemos inicializar los centroides?

Un primer enfoque simple y ampliamente usado, consiste en inicializar los centroides con \(k\) puntos aleatorios del conjunto de datos. Una de las principales consecuencias derivadas de este tipo de inicialización de centroides es que implica que el proceso de clustering será no determinista. Es decir, diferentes ejecuciones del algoritmo pueden conducir a diferentes modelos y, lógicamente, a diferentes niveles de calidad del resultado.

Esta característica, que en principio podría parecer un inconveniente importante, puede ser superada gracias a la simplicidad y efectividad del método, que permite ejecutar el algoritmo múltiples veces con distintos centroides y escoger el mejor resultado.

2. ¿Cómo podemos calcular la distancia entre cada punto y los centroides?

El algoritmo k-medias, generalmente, utiliza la distancia euclídeana. Aunque es posible utilizar cualquier otra métrica.

3. ¿Cómo podemos recalcular los centroides a partir de las instancias que forman su grupo o partición?

El valor de cada centroide es calculado como la media de todos los puntos que pertenecen a este segmento (de aquí el nombre del algoritmo, k-medias). Por lo tanto, este algoritmo solo es aplicable a atributos continuos. En caso de tener atributos no continuos en el conjunto de datos, deberíamos aplicar alguna transformación previa.

4. ¿Cómo podemos establecer la condición de parada del algoritmo?

La convergencia del algoritmo es la condición de parada natural. Esta se alcanza cuando no hay cambios en el recálculo de los centroides durante una iteración completa, provocando que no haya alteraciones en la distribución de las particiones o clústeres. Es decir, se llega a una situación de estabilidad en las particiones.

En la práctica, no es necesario que el método converja, y generalmente es suficiente cuando se tiene pequeños cambios en la pertenencia de las instancias en los distintos grupos o particiones, es decir, cuando estamos próximos a una situación de convergencia.

5. ¿Cómo seleccionar un valor \(k\)?

Una de las principales características de k-medias es que se le debe indicar el número de particiones o clústeres a identificar. Aunque esto le permite mantener un buen nivel de simplicidad y eficiencia, puede ser considerado un inconveniente importante. Identificar el número correcto de particiones \(k\) no es una tarea elemental, en muchas ocasiones el investigador se puede proponer un valor natural de \(k\) por las características de la base de datos. Aun así, aunque no se tenga un valor óptimo de \(k\), en muchos casos el rango de valores no es muy grande, con lo cual se pueden probar diferentes valores y seleccionar el que proporcione mejores resultados. A continuación hablaremos de los tres método más usados para esta tarea.

1. Método de silueta

El método de silueta mide la distancia de separación entre los clústers. El número óptimo de clústers \(k\) es aquel que maximiza la media de los coeficientes de silueta para un rango de valores de \(k\).

Esta medida de distancia se encuentra en el rango de \([-1, 1]\). Los coeficientes de silueta cercanos a +1 indican que la observación se encuentra lejos de los clústers vecinos. Un valor del coeficiente de 0 indica que la observación está muy cerca o en la frontera de decisión entre dos clústers. Valores negativos indican que esas muestras quizás estén asignadas al clúster erróneo.

2. Método del codo

Este método utiliza la distancia media de las observaciones a su centroide, buscando minimizar la suma de cuadrados dentro del clúster. Es decir, la suma de las distancias de cualquier punto a su centroide más cercano. Este criterio busca busca el valor \(k\) que satisfaga que un incremento de \(k\), no mejore sustancialmente la distancia media.

3. Método de brecha

Este método es similar al método del codo, pero ahora la finalidad es encontrar la mayor diferencia o distancia que hay entre los diferentes. Es decir, maximizar la suma de distancias entre sus centroides. En este caso estaríamos priorizando tener segmentos los más alejados posible entre sí, es decir, tener segmentos lo más diferenciados posible.

Cualquiera de los criterios anteriorespuede usarse para establecer un valor adecuado para el parámetro \(k\). Pero es importante resaltar, que dependiendo de la base de datos, estos métodos pueden generar distintos valores de \(k\).

En R se pueden visualizar los tres métodos por medio de la función viz_nbclust() de la librería factoextra . Esta función permite cambiar el método: sulieta (silhouette), codo (wss) y brecha (gap_stat). Además, también se puede cambiar la función para encontrar el número óptimo de clúster. Esto son: hcut, kmeans, cluster::pam, cluster::clara, cluster::fanny.

Por otra parte, es importante resaltar que en R, existe la función NbClust() que determina el número de clúster de manera automatica probando con 30 índices diferentes.

Métodos derivados de k-medias

El uso de la distancia media para el cálculo de los centroides proporciona un modelo simple y eficiente, y además proporciona resultados fácilmente interpretables. Aun así, tiene algunas limitaciones. Por ejemplo, cuando los datos están distribuidos de forma asimétrica y contienen valores extremos (outliers) en algunos grupos, hara que el calculo de la distancia media proporcione muchos errores. Para solventar este problema, en la literatura se han propuesto algunas alternativas.

k-mediana

El algoritmo k-mediana (k-medians) es una variación del método anterior, donde se sustituye el cálculo basado en el valor medio, por el valor de la mediana. Aunque la obtención del valor de la mediana requiere mayor complejidad computacional, es preferible en determinados contextos.

Al igual que en el caso del k-medias, este método puede implementarse mediante distintas métricas de proximidad, aunque existe el riesgo de perder la garantía de convergencia. Es por ello, que una de las métricas de distancia más empleadas en este algoritmo es la distancia de Manhattan.

k-medoides

El algoritmo k-medoides (k-medoids) propone el recálculo de los centros a partir de grupos que presentan un valor de proximidad mínimo respecto a los demás grupos de la partición. La identificación de estos puntos, conocidos como medoides, es computacionalmente más costosa que las aproximaciones vistas anteriormente.

Tambien es conocido por sus siglas en inglés pam (partición alrededor de los medoides) , en el que, cada grupo está representado por uno de los objetos en el grupo. Existe también una variante de pam llamado clara (agrupación de grandes aplicaciones), que se utiliza para analizar grandes conjuntos de datos.

Aunque esta aproximación sea considerablemente más costosa computacionalmente, es preferible en algunos casos donde puede ser interesante que los puntos de centro sean puntos del conjunto de datos. Además, este método es particularmente robusto al ruido y a los valores extremos.

En R, es posible comparar los métodos vistos. La librería clValid() tiene programado una lista extensa de métodos que se podrían evaluar para validar el método que obtiene mejores resultados.

Ejemplos en R

Ejemplo 1:

Para este primer ejemplo, generaremos seis conjuntos de datos normales multivariados.

library(mvtnorm)
# Generamos los datos
set.seed(198)
s1 <- matrix(c(1,.2,0.2,1), ncol=2)
x1 <- rmvnorm(50, mean = c(2, 3),s1)
s2 <- matrix(c(1,0,0,2), ncol=2)
x2 <- rmvnorm(50, mean = c(3, 9),s2)
s3 <- matrix(c(2,0.5,0.5,1.5), ncol=2)
x3 <- rmvnorm(100, mean = c(10, 5),s3)
s4 <- matrix(c(3,0.8,0.8,1.5), ncol=2)
x4 <- rmvnorm(100, mean = c(12, 10),s4)
s5 <- matrix(c(2,0.7,0.7,1.5), ncol=2)
x5 <- rmvnorm(80, mean = c(5, 5),s5)
s6 <- matrix(c(1,0.5,0.5,1), ncol=2)
x6 <- rmvnorm(100, mean = c(2, 12),s6)
x  <- c(x1[,1],x2[,1],x3[,1],x4[,1],x5[,1],x6[,1])
y  <- c(x1[,2],x2[,2],x3[,2],x4[,2],x5[,2],x6[,2])
data <- data.frame(scale(cbind(x,y)))
summary(data)
##        x                 y           
##  Min.   :-1.6848   Min.   :-2.02963  
##  1st Qu.:-0.8983   1st Qu.:-0.88012  
##  Median :-0.3038   Median : 0.03345  
##  Mean   : 0.0000   Mean   : 0.00000  
##  3rd Qu.: 0.9348   3rd Qu.: 0.92221  
##  Max.   : 2.4278   Max.   : 1.90970
# Gráfica de los datos
library(ggplot2)
Conjuntos <- factor(c(rep(1,50),rep(2,50),rep(3,100),rep(4,100),rep(5,80),rep(6,100)))
ggplot() + geom_point(aes(x =x, y = y,color=Conjuntos), data = data,alpha = 0.5) + ggtitle("Datos")

Se puede observar que el conjunto de datos se podría partir entre 2 a 6 grupos. El objetivo de los modelos de particionamiento será entonces identificar de manera precisa tales grupos y, sobre todo, las fronteras de cada uno, de manera que cada observación pueda ser asociada a una clase específica.

Número óptimo de particiones

Primero, encontremos el número optimo de clúster. Trabajaremos con la función hcut() , pero como ya lo hemos mencionado, podríamos aplicar cualquier otra función.

1. Método de silueta 

library(factoextra)
fviz_nbclust(data,kmeans, method = "silhouette",k.max = 10)+
  labs(title= "Número optimo de cluster") + 
  xlab("Valor de k ") +
  ylab("Promedio de silueta")

2. Método del codo 

fviz_nbclust(data, hcut, method = "wss",k.max = 10)+
  labs(title= "Número optimo de cluster") + 
  xlab("Valor de k ") +
  ylab("Suma de cuadrados entre grupos")

3. Método de brecha 

fviz_nbclust(data,hcut, method = "gap_stat",k.max = 10)+
  labs(title= "Número optimo de cluster") + 
  xlab("Valor de k ") +
  ylab("Suma de cuadrados entre grupos")

Al parecer, para este conjunto de datos se puede encontrar por los métodos anteriores que el número óptimo de clúster es \(k=4\).

Métodos de particionamiento

Ahora, apliquemos los diferentes algoritmos de particionamiento.

1. k-medias 

library(stats)
set.seed(1234)
kmedias <- kmeans(data,  centers = 4)

data$cluster <- kmedias$cluster
ggplot() + geom_point(aes(x = x, y = y, color = cluster), data = data, size = 2) +
  scale_colour_gradientn(colours=rainbow(4)) +
  geom_point(aes(x = kmedias$centers[, 1], y = kmedias$centers[, 2]), color = 'black', size = 3) + 
  ggtitle("Clasificación óptima") + 
  xlab('X') + ylab('Y')

2. k-mediana 

library(Kmedians)
kmediana <- Kmedians(data,nclust = 4)
data$cluster <- kmediana$allresults$cluster
ggplot() + geom_point(aes(x = x, y = y, color = cluster), data = data, size = 2) +
  scale_colour_gradientn(colours=rainbow(4)) +
  geom_point(aes(x = kmediana$bestresult$centers[, 1], y = kmediana$bestresult$centers[, 2]), color = 'black', size = 3) + 
  ggtitle("Clasificacion optima") + 
  xlab('X') + ylab('Y')

3. k-mediodes (pam) 

library("cluster")
kpam <- pam(data,k=4)
data$cluster <- kpam$cluster
ggplot() + geom_point(aes(x = x, y = y, color = cluster), data = data, size = 2) +
scale_colour_gradientn(colours=rainbow(4)) +
geom_point(aes(x = kpam$medoids[,1], y = kpam$medoids[,2]), color = 'black', size = 3) + 
ggtitle("Clasificación óptima") + 
xlab('X') + ylab('Y')

3. k-mediodes (clara) 

kclara <- clara(data,k=4)
data$cluster <- kclara$cluster
ggplot() + geom_point(aes(x = x, y = y, color = cluster), data = data, size = 2) +
scale_colour_gradientn(colours=rainbow(4)) +
geom_point(aes(x = kclara$medoids[,1], y = kclara$medoids[,2]), color = 'black', size = 3) + 
ggtitle("Clasificación óptima") + 
xlab('X') + ylab('Y')

Ejemplo 2:

Para este segundo ejemplo, trabajaremos la base de datos USArrests , la cual contiene estadísticas de arrestos en los diferentes estados.

data.scaled <- scale(USArrests,
                     center = TRUE,
                     scale = TRUE)
head(data.scaled)
##                Murder   Assault   UrbanPop         Rape
## Alabama    1.24256408 0.7828393 -0.5209066 -0.003416473
## Alaska     0.50786248 1.1068225 -1.2117642  2.484202941
## Arizona    0.07163341 1.4788032  0.9989801  1.042878388
## Arkansas   0.23234938 0.2308680 -1.0735927 -0.184916602
## California 0.27826823 1.2628144  1.7589234  2.067820292
## Colorado   0.02571456 0.3988593  0.8608085  1.864967207

Número óptimo de particiones

Para este ejemplo, trabajaremos con la función NbClust() que prueba 30 índices diferentes para encontrar el número óptimo de particiones

library(NbClust)
set.seed(1234)
res.nbclust <- NbClust(data.scaled, distance = "euclidean",
                  min.nc = 2, max.nc = 10, 
                  method = "complete", index ="all")

## *** : The Hubert index is a graphical method of determining the number of clusters.
##                 In the plot of Hubert index, we seek a significant knee that corresponds to a 
##                 significant increase of the value of the measure i.e the significant peak in Hubert
##                 index second differences plot. 
## 

## *** : The D index is a graphical method of determining the number of clusters. 
##                 In the plot of D index, we seek a significant knee (the significant peak in Dindex
##                 second differences plot) that corresponds to a significant increase of the value of
##                 the measure. 
##  
## ******************************************************************* 
## * Among all indices:                                                
## * 9 proposed 2 as the best number of clusters 
## * 4 proposed 3 as the best number of clusters 
## * 6 proposed 4 as the best number of clusters 
## * 2 proposed 5 as the best number of clusters 
## * 1 proposed 8 as the best number of clusters 
## * 1 proposed 10 as the best number of clusters 
## 
##                    ***** Conclusion *****                            
##  
## * According to the majority rule, the best number of clusters is  2 
##  
##  
## *******************************************************************

Por lo tanto, siguiendo la recomendación dada, tomaremos \(k=2\).

Número óptimo de particiones

Probaremos ahora con la función clValid() , la cual tiene programada 10 métodos diferentes y además de evaluar el método, te proporciona otra evaluación sobre el número óptimo de particiones.

library(clValid)
validclus  <- clValid(data.scaled, nClust = 2:5, 
                clMethods = c("hierarchical","kmeans","diana","fanny","pam","clara","agnes"),
                validation = "internal")
summary(validclus)
## 
## Clustering Methods:
##  hierarchical kmeans diana fanny pam clara agnes 
## 
## Cluster sizes:
##  2 3 4 5 
## 
## Validation Measures:
##                                  2       3       4       5
##                                                           
## hierarchical Connectivity   6.6437  9.5615 13.9563 22.5782
##              Dunn           0.2214  0.2214  0.2224  0.2046
##              Silhouette     0.4085  0.3486  0.3637  0.3213
## kmeans       Connectivity   6.6437 13.6484 16.2413 24.6639
##              Dunn           0.2214  0.2224  0.2224  0.1983
##              Silhouette     0.4085  0.3668  0.3573  0.3377
## diana        Connectivity   6.6437 11.3635 20.9694 23.5623
##              Dunn           0.2214  0.2224  0.1604  0.1709
##              Silhouette     0.4085  0.3734  0.3318  0.3301
## fanny        Connectivity  12.9004 25.1730 22.8071 39.5016
##              Dunn           0.1805  0.0970  0.1233  0.1604
##              Silhouette     0.3870  0.2597  0.2156  0.2286
## pam          Connectivity   6.6437 13.8302 20.4421 29.5726
##              Dunn           0.2214  0.1376  0.1849  0.1849
##              Silhouette     0.4085  0.3144  0.3390  0.3105
## clara        Connectivity   6.6437 13.8302 20.4421 29.5726
##              Dunn           0.2214  0.1376  0.1849  0.1849
##              Silhouette     0.4085  0.3144  0.3390  0.3105
## agnes        Connectivity   6.6437  9.5615 13.9563 22.5782
##              Dunn           0.2214  0.2214  0.2224  0.2046
##              Silhouette     0.4085  0.3486  0.3637  0.3213
## 
## Optimal Scores:
## 
##              Score  Method       Clusters
## Connectivity 6.6437 hierarchical 2       
## Dunn         0.2224 hierarchical 4       
## Silhouette   0.4085 hierarchical 2

Para entender los resultados, debemos buscar el método que minimiza la conectividad, y maximice tanto el índice dunn como el ancho de la silueta.

Por lo tanto, parece que el clúster jerárquico supera a los otros algoritmos de agrupación en cada medida de validación, para casi todos los números de clusters evaluados. Además, se puede corroborar la conclusión de que \(k=2\) es el mejor valor para el número de grupos.

Por lo tanto, si deseamos quedarnos con una clasificación óptima sería la siguiente.

dist      <- dist(data.scaled,method = "euclidean")
modelo    <- hclust(dist, method = "complete")

fviz_dend(modelo, cex = 0.5, k=2, 
          rect = TRUE,  
          k_colors = "jco",
          rect_border = "jco", 
          rect_fill = TRUE,
          horiz = TRUE,
          ggtheme = theme_bw())