Introducción a las técnicas multivariantes de clasificación

El objetivo de estas técnicas es agrupar los objetos de interés (_observaciones_) en grupos internamente homogéneos, de forma que los objetos del mismo grupo serán similares respecto a los criterios establecidos.

En algunos casos interesará que el resultado sea una partición de los objetos en \(k\) grupos de forma que cada objeto pertenezca solo a uno de los grupos. Estas técnicas se conocen como no jerárquicas.

En otras ocasiones interesará establecer una jerarquía en la estructura de relaciones entre los objetos, produciendo, a posteriori, distintas particiones en función del grado de similitud entre objetos. Estas son técnicas de clasificación jerárquica.

Puesto que estas técnicas sirven para asignar a los objetos a diferentes grupos, este tipo de técnicas también la puedes encontrar bajo la denominación de análisis de conglomerados o clustering.

Clasificación iterativa

Cómo se crean los grupos

El algoritmos más conocido y utilizado —aunque no el único— es K-medias, Kmeans en inglés.

Su fundamento es el siguiente:

Partimos de una serie de observaciones distribuidas en el espacio en función de su matriz de distancias.

Se eligen de forma aleatoria una serie de puntos, que pueden pertenecer o no al conjunto de datos. Tantos puntos como grupos — clusters —. Estos puntos actúan como centroides. Serán el “centro del universo” de esa clase en esa iteración.

Las observaciones se irán asignando al centroide más cercano…

…creándose una primera clasificación de la primera iteración.

A continuación se calcula nuevamente los centroides, esta vez como el punto medio de cada clase.

Es decir, se calcula la distancia media entre todos los puntos pertenecientes a una clase y se vuelve a asignar las observaciones al centroide más cercano —al nuevo centroide—.

Esto se repetirá hasta que haya convergencia, lo que quiere decir que ya no habrá nuevas reasignaciones o hasta que se alcance el máximo de iteraciones fijadas, generando así una clasificación de las observaciones en grupos homogéneos.

El algoritmo de k-medias converge en un óptimo local, eso significa que el resultado que obtengamos no tiene por que ser el más óptimo. Si cambiamos los centroides de origen muy posiblemente obtendremos diferentes resultados. Por tanto la forma en que se inicializan los centroides es crítica, y aunque no entremos en ello por quedar fuera del alcance de este curso hay métodos como kmeans++ para optimizar esto.

Más información en este libro y en esta página de la que se han sacado estas imágenes.

Kmeans cómo funciona

En el siguiente enlace hay un vídeo muy ilustrativo de cómo funciona Kmeans, por si queda alguna duda.

Vídeo Explicativo Kmeans

Otro ejemplo animado de formación de los grupos

Coge el código de este chunk y ejecútalo en tu Rstudio. Recuerda que tienes que tener cargada la librería animation

# demostración del algoritmo. Ejecutar en Rstudio
mydata <- iris[ , -5 ]
kmeans.ani( mydata[,1:4], centers = 3 )

Ejemplo práctico clustering iterativo

# Preparar datos
data ( iris )
mydata <- iris[ , -5 ]
head( mydata )
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1          5.1         3.5          1.4         0.2
2          4.9         3.0          1.4         0.2
3          4.7         3.2          1.3         0.2
4          4.6         3.1          1.5         0.2
5          5.0         3.6          1.4         0.2
6          5.4         3.9          1.7         0.4
kc <- kmeans( mydata, 3 , nstart = 25, iter.max = 100)
# colores por grupos
kc <- kmeans( mydata, 3 , nstart = 25, iter.max = 100 )
clusplot( mydata, kc$cluster, color = TRUE, 
          shade = FALSE, labels = 4, lines = 1,
          col.p = kc$cluster )

# colores por especies
mycolor <- as.character( factor( iris$Species, 
                                 labels = c( 1:3 ) ) )
clusplot( mydata, kc$cluster, color = TRUE, 
         shade = FALSE, labels = 4, lines = 1,
         col.p = mycolor )

table( iris$Species, kc$cluster, dnn = c("Original", "Kmeans" ) )
            Kmeans
Original      1  2  3
  setosa     50  0  0
  versicolor  0 48  2
  virginica   0 14 36

Datos generados por kmeans:

kmeans(X, k)\(clusters | grupo al cual ha sido clasificado cada individuo| |kmeans(X, k)\)centers centro (vector de medias) de cada grupo
kmeans(X, k)\(withinss| suma de cuadrados dentro de cada grupo| |kmeans(X, k)\)totss suma de cuadrados total
kmeans(X, k)\(tot.withinss| suma de cuadrados total dentro| |kmeans(X, k)\)betweenss suma de cuadrados entre grupos
kmeans(X, k)$size tamaño de los grupos

Determinar número de grupos

# Clustering iterativo kmeans
x <- NULL
for ( i in 1:10 ){
  kc     <- kmeans( mydata, i )
  x[ i ] <- kc$tot.withinss
}
plot( c( 1:10 ), x, type = "b" )

La suma de cuadrados intragrupos disminuye a medida que aumenta \(k\). El número de grupos debería ser aquel en el que la disminución de la suma de cuadrados intragrupos no sea significativa.

En el gráfico de sedimentación equivaldría al punto donde se genera un “codo” en la curva.

La función fviz_bnclust para determinar el número de grupos

Tiene varios métodos, algunos sugieren número de grupos y en otros hay que deducirlo de la gráfica

fviz_nbclust(mydata, kmeans, method="wss")

fviz_nbclust(mydata, kmeans, method="silhouette")

fviz_nbclust(mydata, kmeans, method="gap_stat")

El paquete NbClust

Utiliza 30 índices para decidir el número de grupos. El número de grupos propuesto por más índices será —teóricamente— el óptimo.

nb <- NbClust(mydata, distance = "euclidean", min.nc = 2,
        max.nc = 10, method = "kmeans")

*** : 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:                                                
* 11 proposed 2 as the best number of clusters 
* 11 proposed 3 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 
 
 
******************************************************************* 
fviz_nbclust(nb)
Among all indices: 
===================
* 2 proposed  0 as the best number of clusters
* 11 proposed  2 as the best number of clusters
* 11 proposed  3 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 .

Otras funciones

K-medoids: la funcion pam

Es más robusto, menos sensible a los outliers, puede usar matrices con disimilaridades no euclideas. Los medioides sería algo parecido, aunque no son las medianas.

El medioide es un elemento perteneciente al grupo, a diferencia de k-means que es un valor medio que puede coincidir o no con la posición de un elemento.

kp <- pam( mydata, 3 )
clusplot( mydata, kp$cluster, color = TRUE, 
         shade = FALSE, labels = 4, lines = 1,
         col.p = mycolor )

table( iris$Species, kp$cluster, dnn = c( "Original", "Pam" ) )
            Pam
Original      1  2  3
  setosa     50  0  0
  versicolor  0 48  2
  virginica   0 14 36

¿Y si tengo una matriz de distancias?

la función pam() también puede trabajar sobre matriz de distancias.

K-medoids: la función clara

Puede ejecutarse sobre matriz de datos o distancias igual que pam(), pero según los autores da mejores resultados con grandes conjuntos de datos.

mydist <- as.matrix( dist( mydata ) )
kcl <- clara( mydist, 3 )
clusplot( mydata, kcl$cluster, color = TRUE, 
          shade = TRUE, labels = 4, lines = 1, col.p = mycolor )

table( iris$Species, kcl$cluster, dnn = c( "Original", "Clara" ) )
            Clara
Original      1  2  3
  setosa     50  0  0
  versicolor  0  3 47
  virginica   0 39 11

La función silhouette para evaluar el número \(k\).

Se utiliza la función Silhouette para determina cómo de bueno o de homogéneo es cada cluster formado.

Rango de valores Interpretación
0.71 – 1.00 Se ha encontrado estructura fuerte
0.51 – 0.70 La estructura es razonable
0.26 – 0.50 La estructura es débil y podría ser artificial
< 0.25 No se ha encontrado una estructura sustancial

¿Qué significan esos valores?

Si \(a_i\) es la disimilaridad —distancia— media desde la observación \(i\) al resto de observaciones de su grupo, y \(b_i\) es la menor distancia media desde esa observación a los individuos de los otros grupos, el valor medio de la función será:

\[\text{Silhouette}_i = \dfrac{ b_i - a_i } {\max( a_i, b_i )}\]

skmeans <- silhouette( kp )
plot( skmeans )

Representación 3d

Colores por grupos

mydataPCA <- PCA( mydata[ , 1:4 ], scale.unit = TRUE, graph = FALSE )
plot3d( mydataPCA$ind$coord[ , 1:3 ], type = "s", 
        col = kp$clustering, radius = 0.1,
        xlim = c(-3, 3 ), ylim = c( -3, 3 ) )

You must enable Javascript to view this page properly.

Colores por species

plot3d( mydataPCA$ind$coord[ , 1:3], type = "s", 
        col = factor( iris$Species, labels = c( 1:3 ) ),
        radius = 0.1, box=F )

You must enable Javascript to view this page properly.

Clasificación jerárquica

Cómo se crea la jerarquía

  • Se crea una estructura jerárquica basada en las distancias entre individuos

  • Aglomerativa: cada observación es su propio grupo y los grupos se van mezclando

  • Divisiva: todas las observaciones están en el mismo grupo y en cada iteración se van dividiendo los grupos.

En las técnicas jerárquicas se parte de una matriz de distancias o disimilaridades, estableciendo las distancias entre grupos según diferentes criterios. Algunos de los más conocidos son:

ward
Este criterio de agrupación mezclará aquel par de grupos que lleven a un incremento mínimo de la varianza del cluster después de ser mezclado.
single
Se mezclarán grupos con la mínima distancia entre los vecinos más próximos
complete
Se mezclarán grupos con la mínima distancia entre los vecinos más alejados

otros criterios incluidos en la función hclust son mcquitty, median, centroid

Uno de los que suele dar mejores resultados es el criterio de ward.

Cómo funciona: los datos

Sepal.Length Sepal.Width
1 5.1 3.5
2 4.9 3.0
3 4.7 3.2
4 4.6 3.1

Cómo funciona: iter 1

          1         2         3
2 0.5385165                    
3 0.5000000 0.2828427          
4 0.6403124 0.3162278 0.1414214

Cómo funciona: iter 2

            1         2
2   0.5385165          
3-4 0.5700877 0.2915476

Cómo funciona: iter 3

              1
3-4-2 0.5350234

Resultado

Dendrograma clásico

mydata <- iris[ , -5 ]
mydist <- dist( mydata, method = "euclidean" )
hclu   <- hclust( mydist, method = "ward.D" )
# convertir objeto cluster en dendrograma, útil para representar aunque no necesario
dendro <- as.dendrogram( hclu ) 
#
plot( dendro , horiz = FALSE, type = "r" ) # type = "t"
rect.hclust( hclu, k = 3 )

Dendrograma personalizado

#muestra de tamaño 20 del data.frame
subdata <- iris[ sample( nrow( iris ), 20 ), ]

hc <- as.phylo( hclust( dist(subdata[ , -5 ] ), 
                method = "ward.D" )  )

mycolor <- as.character( factor( subdata$Species, labels = c( 1:3 ) ) )
plot( hc, type = "phylogram", cex = .8, tip.color = mycolor )

plot( hc, type = "cladogram", cex = .8, tip.color = mycolor )

plot( hc, type = "unrooted",  cex = .8, tip.color = mycolor )

plot( hc, type = "fan",       cex = .8, tip.color = mycolor )

plot( hc, type = "radial",    cex = .8, tip.color = mycolor )

Combinando técnicas multivariantes

Clustering de los resultados

También se puede combinar un PCA con un clustering jerárquico para intentar obtener alguna pista que nos ayude a interpretar los resultados.

biom <- read.table( "http://ares.inf.um.es/00Rteam/pub/clas/data/biom2003.dat" )
Xdf <- biom[ , 2:7 ] 
biomPCA <- PCA( Xdf , graph = FALSE )
# k-means 3 grupos sobre datos originales
kmRaw <- kmeans( Xdf, 3 )
clusplot( Xdf, kmRaw$cluster, color = TRUE, 
         shade = FALSE, labels = 4, lines = 1,
         col.p = kmRaw$cluster )

# k-means 3 grupos sobre matriz de PCA
kmPCA <- kmeans( biomPCA$ind$coord, 3 )

clusplot( Xdf, kmPCA$cluster, color = TRUE, 
          shade = FALSE, labels = 4, lines = 1, 
          col.p = kmPCA$cluster )

table( kmRaw$cluster, kmPCA$cluster, dnn = c( "Raw", "PCA" ) )
   PCA
Raw  1  2  3
  1  0 28  0
  2 31  0 15
  3  2  2 20

Clustering jerárquico sobre PCA

biomClust <- hclust( dist( biomPCA$ind$coord ), method = "ward.D" )
plot( biomClust, xlab = "", sub = "" )

clusters <- as.factor( cutree( biomClust, k = 3 ) )
clusplot( Xdf, clusters, color = TRUE, 
          shade = FALSE, labels = 4, lines = 1, 
          col.p = clusters, main = "hclust sobre PCA" )

# plot( biomPCA$x, type = "n" )
# text( biomPCA$x[ , 1 ], biomPCA$x[ , 2 ], 
      # labels = rownames( biomPCA$x ), col = c( clusters ) )
# Clustering jerárquico sobre datos RAW
hcRaw <- hclust( dist ( Xdf ), method = "ward.D" )
clusterRaw <- as.factor( cutree( hcRaw, k = 3 ) )
clusplot( Xdf, clusterRaw, color = TRUE, 
          shade = FALSE, labels = 4, lines = 1, 
          col.p = clusterRaw, main = "hclust sobre datos originales" )

table( clusterRaw, clusters, dnn = c( "Raw", "PCA" ) )
   PCA
Raw  1  2  3
  1 16 27  1
  2  1  0 33
  3 19  0  1

Factominer en acción

La función HCPC

cl_hcpc <- HCPC( Xdf, nb.clust = 3 )

cl_hcpc <- HCPC( Xdf, nb.clust = -1 )

# Modo interactivo, ejecutar en Rstudio
cl_hcpc <- HCPC( Xdf)

Manos a la obra

En esta ocasión, vas a hacer un ejercicio de clasificación iterativa, con los datos biom2003.dat que ya hemos utilizado en otros módulos. Verás si hay diferencia entre normalizar y no normalizar los datos. Finalmente harás una clasificación jerárquica con los mismos datos y comprobarás los diferentes resultados obtenidos al aplicar distintos algoritmos aglomerativos. Vamos a ello.

Referencias