Capítulo 4: K-Means Clustering

Computar K-means clustering en R

Datos

Los datos solo deben contener variables continuas ya que el algoritmo utiliza la media de las variables.

Primeramente se debe escalar los datos utilizando la función scale ().

data("USArrests")
df <- scale(USArrests)

head(df, n=3)

##             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

Tabla 1. Primeras tres líneas de datos para los arrestos en Estados Unidos.

En el código anterior se utiliza la función head () para definir los datos que serán usados. En este caso se usan las primeras 3 filas.

Paquetes y funciones requeridos

La función principal es kmeans () cuyos elementos son:

X: matriz numérica, data frame o vector numérico
centers: posibles valores de los grupos (k)
iter.max: Número máximo de iteraciones o repeticiones permitidas
nstart: El número de posiciones aleatorias cuando centers es un número

Para crear un gráfico con la función kmeans() es necesario instalar y cargar el paquete factoextra

install.packages("factoextra")

## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.3'
## (as 'lib' is unspecified)

library(factoextra)

## Loading required package: ggplot2

## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

Estimación del numéro optimo de grupos

La agrupación de kmeans requiere a los usuarios especificar el número de grupos para ser generados.

Lo principal para esta función es determinar la cantidad adecuada de grupos, esto será respondido a través del capítulo.

library(factoextra)
fviz_nbclust(df,kmeans, method = "wss") +
  geom_vline(xintercept = 4, linetype = 2)

Figura 1. Número de grupos a traves de la función fviz_nbclust ()

La función fviz_nbclust () permite una solución a la determinación del número óptimo de grupos para la función kmeans. En la figura se muestra la varianza entre los grupos y se observa como decrece con el aumento de k. La línea en k = 4 indica que los grupos después de el cuarto tienen muy poco valor y por lo tanto son innecesarios.

Computando k-means clustering

El algoritmo k-means clustering siempre empieza con una cantidad k aleatoria de valores. Debido a esto, se recomienda usar la función set.seed () para determinar una semilla para el generador de números aleatorios.

Esto se da ya que se busca tener resultados reproducibles por lo que, teniendo en cuenta el resultado anterior (k = 4), se utiliza:

set.seed(123)
km.res <- kmeans(df,4, nstart = 25)
print(km.res)

## 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.9615407 -1.1066010 -0.9301069 -0.96676331
## 3 -0.4894375 -0.3826001  0.5758298 -0.26165379
## 4  0.6950701  1.0394414  0.7226370  1.27693964
## 
## Clustering vector:
##        Alabama         Alaska        Arizona       Arkansas     California 
##              1              4              4              1              4 
##       Colorado    Connecticut       Delaware        Florida        Georgia 
##              4              3              3              4              1 
##         Hawaii          Idaho       Illinois        Indiana           Iowa 
##              3              2              4              3              2 
##         Kansas       Kentucky      Louisiana          Maine       Maryland 
##              3              2              1              2              4 
##  Massachusetts       Michigan      Minnesota    Mississippi       Missouri 
##              3              4              2              1              4 
##        Montana       Nebraska         Nevada  New Hampshire     New Jersey 
##              2              2              4              2              3 
##     New Mexico       New York North Carolina   North Dakota           Ohio 
##              4              4              1              2              3 
##       Oklahoma         Oregon   Pennsylvania   Rhode Island South Carolina 
##              3              3              3              3              1 
##   South Dakota      Tennessee          Texas           Utah        Vermont 
##              2              1              4              3              2 
##       Virginia     Washington  West Virginia      Wisconsin        Wyoming 
##              3              3              2              2              3 
## 
## Within cluster sum of squares by cluster:
## [1]  8.316061 11.952463 16.212213 19.922437
##  (between_SS / total_SS =  71.2 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"

Tabla 2. Agrupaciones k-mean de 4 tamaños

En el algoritmo fue necesario especificar que se debe comenzar con 25 asignaciones aleatorias para posteriormente escoger el mejor resultado. La tabla muestra losvectores de agrupación en donde se asignan valores del 1 al 4

aggregate(USArrests, by=list(cluster = km.res$cluster), mean)

##   cluster   Murder   Assault UrbanPop     Rape
## 1       1 13.93750 243.62500 53.75000 21.41250
## 2       2  3.60000  78.53846 52.07692 12.17692
## 3       3  5.65625 138.87500 73.87500 18.78125
## 4       4 10.81538 257.38462 76.00000 33.19231

dd <- cbind(USArrests, cluster = km.res$cluster)
head(dd)

##            Murder Assault UrbanPop Rape cluster
## Alabama      13.2     236       58 21.2       1
## Alaska       10.0     263       48 44.5       4
## Arizona       8.1     294       80 31.0       4
## Arkansas      8.8     190       50 19.5       1
## California    9.0     276       91 40.6       4
## Colorado      7.9     204       78 38.7       4

Tabla 3 & 4. Media de variables con datos originales y sus clasificaciones

En esta figura se observa la agrupación de los datos utilizando el grupo de datos original o inicial.

Acceso a los resultados de la función kmeans ()

La función devuelve una lista de componentes los cuales son:

Cluster: Un vector de números enteros que indican el grupo en el que cada punto fue asignado
Centers: Una matriz de medias del grupo
totss: La suma total de cuadrados. La suma de la varianza total en los datos
Withinss: vector de suma de cuadrados dentro del grupo. Un componente por grupo
Tot.withinss: Total de la suma de cuadrados dentro del grupo
betweenss: La suma de cuadrados entre grupos
size: numero de observaciones en cada grupo

Para acceder a estos componentesse utiliza:

km.res$cluster

##        Alabama         Alaska        Arizona       Arkansas     California 
##              1              4              4              1              4 
##       Colorado    Connecticut       Delaware        Florida        Georgia 
##              4              3              3              4              1 
##         Hawaii          Idaho       Illinois        Indiana           Iowa 
##              3              2              4              3              2 
##         Kansas       Kentucky      Louisiana          Maine       Maryland 
##              3              2              1              2              4 
##  Massachusetts       Michigan      Minnesota    Mississippi       Missouri 
##              3              4              2              1              4 
##        Montana       Nebraska         Nevada  New Hampshire     New Jersey 
##              2              2              4              2              3 
##     New Mexico       New York North Carolina   North Dakota           Ohio 
##              4              4              1              2              3 
##       Oklahoma         Oregon   Pennsylvania   Rhode Island South Carolina 
##              3              3              3              3              1 
##   South Dakota      Tennessee          Texas           Utah        Vermont 
##              2              1              4              3              2 
##       Virginia     Washington  West Virginia      Wisconsin        Wyoming 
##              3              3              2              2              3

head(km.res$cluster, 4)

##  Alabama   Alaska  Arizona Arkansas 
##        1        4        4        1

Tabla 5. número de grupo para cada observación

km.res$size

## [1]  8 13 16 13

Tabla 6. Tamaño de grupos

km.res$centers

##       Murder    Assault   UrbanPop        Rape
## 1  1.4118898  0.8743346 -0.8145211  0.01927104
## 2 -0.9615407 -1.1066010 -0.9301069 -0.96676331
## 3 -0.4894375 -0.3826001  0.5758298 -0.26165379
## 4  0.6950701  1.0394414  0.7226370  1.27693964

Tabla 7. Media de los grupos

Visualización de grupos k-means

Para visualizar un conjunto de datos con más de 2 variables se realiza un Análisis de componentes principales, o PCA según sus siglas en inglés.

La función fviz_cluster() se puede utilizar para visualizar las agrupaciones. En el plot resultante las observaciones se representan por medio de puntos .

fviz_cluster(km.res, data = df,
             palette = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"), 
             ellipse.type = "euclid",
             star.plot = TRUE, 
             repel = TRUE, 
             ggtheme = theme_minimal()
             )

Figura ??. Visualización de grupos K-means

En la figura anterior se muestra la representación visual de los datos por medio de puntos y elipses. Aquellos datos que poseen número de grupos iguales están del mismo color y son agrupados por medio de una elipse de concentración, como es el caso de Mississippi, Carolina del Norte, Carolina del Sur, Tenesse, Alabama, Arkansas, Louisiana y Georgia con 1.

Capítulo 5: K-medoids

El algoritmo de K-medoids está relacionado a k-means al dividir los datos en cierta cantidad de agrupaciones. Para K-medoidscada agrupación se ve representada por puntos.

Las principales características de este algoritmo es que las diferencias entre los datos dentro de cada grupo son mpinimas. Se toman valores centrales de cada conjunto los cuales pueden ser representativos de este debido a lo establecido anteriormente.

Computar PAM en R

PAM o Partition Around Medoids se basa en la búsqueda de objetos representativos de los K grupos en el conjunto de datos. Por medio de PAM se busca encontrar una función objetiva la cuál determina la suma total de desigualdades que presentan los datos con su medoide más cercano

Datos

data("USArrests")
df <- scale(USArrests)
head(df, n = 3)

##             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

Tabla ??. Conjunto de datos a utilizar para K-medoids

Paquetes y funciones necesarios

Para k-medoids se utilizarán las funciones pam() y pamk() para computar PAM

El formato de la función es:

pam (x, k, metric = “euclidean”, stand = FALSE)

En donde:

X: Valores posibles
k: Número de agrupaciones
metric: Las métricas de distancia a usar
stand: valor lógico

library(cluster)
library(factoextra)

Estimación del número óptimo de agrupaciones

Para la estimación del número óptimo de agrupaciones será utilizado el método Silhouette para comptar el algoritmo PAM usando diferentes valores de grupo.

library(cluster)
library(factoextra)
fviz_nbclust(df, pam, method = "silhouette") + 
  theme_classic()

Figura ??. Número ideal de agrupaciones para K-Medoids usando el método Silhouette

Computar agrupaciones PAM

Para computar el algoritmo con k = 2:

pam.res <- pam(df, 2)
print(pam.res)

## Medoids:
##            ID     Murder    Assault   UrbanPop       Rape
## New Mexico 31  0.8292944  1.3708088  0.3081225  1.1603196
## Nebraska   27 -0.8008247 -0.8250772 -0.2445636 -0.5052109
## Clustering vector:
##        Alabama         Alaska        Arizona       Arkansas     California 
##              1              1              1              2              1 
##       Colorado    Connecticut       Delaware        Florida        Georgia 
##              1              2              2              1              1 
##         Hawaii          Idaho       Illinois        Indiana           Iowa 
##              2              2              1              2              2 
##         Kansas       Kentucky      Louisiana          Maine       Maryland 
##              2              2              1              2              1 
##  Massachusetts       Michigan      Minnesota    Mississippi       Missouri 
##              2              1              2              1              1 
##        Montana       Nebraska         Nevada  New Hampshire     New Jersey 
##              2              2              1              2              2 
##     New Mexico       New York North Carolina   North Dakota           Ohio 
##              1              1              1              2              2 
##       Oklahoma         Oregon   Pennsylvania   Rhode Island South Carolina 
##              2              2              2              2              1 
##   South Dakota      Tennessee          Texas           Utah        Vermont 
##              2              1              1              2              2 
##       Virginia     Washington  West Virginia      Wisconsin        Wyoming 
##              2              2              2              2              2 
## Objective function:
##    build     swap 
## 1.441358 1.368969 
## 
## Available components:
##  [1] "medoids"    "id.med"     "clustering" "objective"  "isolation" 
##  [6] "clusinfo"   "silinfo"    "diss"       "call"       "data"

Tabla ?.Algoritmo PAM con k = 2

dd <- cbind(USArrests, cluster = pam.res$cluster)
head(dd, n = 3)

##         Murder Assault UrbanPop Rape cluster
## Alabama   13.2     236       58 21.2       1
## Alaska    10.0     263       48 44.5       1
## Arizona    8.1     294       80 31.0       1

Tabla ?. Datos originales con algoritmo PAM con k =2

Acceso a los resultados de la función pam()

Luego de aplicar la función pam(), los componentes que esta devuelve son:

Medoids: Objetos que representan las agrupaciones
Clustering: Un vector que contiene el número de grupo de cada objeto

Para acceder a estos componentes es necesario:

pam.res$medoids

##                Murder    Assault   UrbanPop       Rape
## New Mexico  0.8292944  1.3708088  0.3081225  1.1603196
## Nebraska   -0.8008247 -0.8250772 -0.2445636 -0.5052109

Tabla ?. medoides del grupo

head(pam.res$clustering)

##    Alabama     Alaska    Arizona   Arkansas California   Colorado 
##          1          1          1          2          1          1

Tabla ?. Número de agrupaciones

Visualización de las agrupaciones PAM

Para visualizar las particiones de los resultados se utiliza la función fviz_cluster (). Esta se encarga de dibujar diagramas de dispersión y sus datos por números de grupos. En caso de que contengas más de 2 variables, el algoritmo de Principal Component Analysis (PCA) se encarga de reducir las dimensiones.

fviz_cluster(pam.res,
             palette = c("#00AFBB", "#FC4E07"), 
             ellipse.type = "t",
             repel = TRUE,
             ggtheme = theme_classic())

Figura ?. Visualización de las agrupaciones PAM

A partir del código anterior surge este diagrama como representación visual de el conjunto de datos. En el código se establece la paleta de colores, el tipo de elipses (que en este caso son de concentración) y se establece que las marcas no se sobrepongan entre sí.

Capítulo 6: CLARA- Agrupando grandes aplicaciones

CLARA, o Clustering Large Appplications es una extención de el método K-medoid. Surge con el proósito de tratar conjuntos de datos con una gran cantidad de objetos para poder reducir con el tiempo de computación y el problema de RAM.

CLARA considera pequeñas muestras de datos con tamaños reducidos y aplica el algoritmo PAM para generar un set de medoides óptimos.

Computando CLARA en R

Datos y preparación

En el siguiente código se generará un conjunto de datos aleatorio usando la función set.seed().

set.seed(1234)
df <- rbind(cbind(rnorm(200,0,8), rnorm(200, 0, 8)),
            cbind(rnorm(300,50,8), rnorm(300,50,8)))
colnames(df) <- c("x", "y")
rownames(df) <- paste0("S", 1:nrow(df))
head(df, nrow = 6)

##             x        y
## S1  -9.656526 3.881815
## S2   2.219434 5.574150
## S3   8.675529 1.484111
## S4 -18.765582 5.605868
## S5   3.432998 2.493448
## S6   4.048447 6.083699

En el algoritmo anterior se generó un conjunto de datos de 500 objetos dividido en 2 agrupaciones, posteriormente se especificó un nombre de columna y de fila y se hizo una vista previa o un avance de los datos.

Paquetes y funciones requeridos

Para computar CLARA, es necesario hacer uso de la función clara(). A partir de sta función se obtienen:

x: datos numéricos de la matriz o data frame
k: número de agrupaciones
metric: las métricas de distacia que serán usadas
stand: valores lógicos
samples: numero de muestras que se extraen del conjunto de datos
pamLike: lógico. Indica si el algoritmo en la función pam() debe ser usado

Computando CLARA

clara.res <- clara(df, 2, samples = 50, pamLike = TRUE)
print(clara.res)

## Call:     clara(x = df, k = 2, samples = 50, pamLike = TRUE) 
## Medoids:
##              x         y
## S121 -1.531137  1.145057
## S455 48.357304 50.233499
## Objective function:   9.87862
## Clustering vector:    Named int [1:500] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ...
##  - attr(*, "names")= chr [1:500] "S1" "S2" "S3" "S4" "S5" "S6" "S7" ...
## Cluster sizes:            200 300 
## Best sample:
##  [1] S37  S49  S54  S63  S68  S71  S76  S80  S82  S101 S103 S108 S109 S118 S121
## [16] S128 S132 S138 S144 S162 S203 S210 S216 S231 S234 S249 S260 S261 S286 S299
## [31] S304 S305 S312 S315 S322 S350 S403 S450 S454 S455 S456 S465 S488 S497
## 
## Available components:
##  [1] "sample"     "medoids"    "i.med"      "clustering" "objective" 
##  [6] "clusinfo"   "diss"       "call"       "silinfo"    "data"

dd <- cbind(df, cluster = clara.res$cluster)
head(dd, n = 4)

##             x        y cluster
## S1  -9.656526 3.881815       1
## S2   2.219434 5.574150       1
## S3   8.675529 1.484111       1
## S4 -18.765582 5.605868       1

La tabla anterior muestra el algoritmo PAM con k = 2. En este se muestran los componentes de la función y su aplicación en los datos originales

clara.res$medoids

##              x         y
## S121 -1.531137  1.145057
## S455 48.357304 50.233499

head(clara.res$clustering, 10)

##  S1  S2  S3  S4  S5  S6  S7  S8  S9 S10 
##   1   1   1   1   1   1   1   1   1   1

Visualizando agrupaciones CLARA

fviz_cluster(clara.res,
             palette = c("#00AFBB", "#FC4E07"),
             ellipse.type = "t",
             geom = "point", pointsize = 1,
             ggtheme = theme_classic())

Capítulo 7: agrupación aglomerativa

La agrupación aglomerativa es una de los tipos más comunes de agrupaciones jerarquicas usadas para juntar con base a la similaridad. Este algoritmo comienza considerando a cada elemento único, a partir de esto, en cada paso del algoritmo los clusters que son más similares se agrupan generando uno más grande.

Pasos para una agrupación aglomerativa herarquica

Preparar los datos
Computar la información de similaritud en cada par de objetos
Usar la función de enlace para agrupar los objetos
Determinar la partición de los datos

Estructura y preparación de los datos

Los datos deben de ser una matriz numérica con filas representando las observaciones individuales y columnas representando variables.

data("USArrests")
df <- scale(USArrests)
head(df, nrow = 6)

##                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

Medidas de similaridad

Para decidir cuáles objetos o agrupaciones deben ser combinadas o divididas se necesitan métodos de similitud.

La función dist() busca computar la distancia entre cada pareja de objetos, esto es conocido como distancia o matriz de disimilaridad

res.dist <- dist(df, method = "euclidean")
as.matrix(res.dist)[1:6, 1:6]

##             Alabama   Alaska  Arizona Arkansas California Colorado
## Alabama    0.000000 2.703754 2.293520 1.289810   3.263110 2.651067
## Alaska     2.703754 0.000000 2.700643 2.826039   3.012541 2.326519
## Arizona    2.293520 2.700643 0.000000 2.717758   1.310484 1.365031
## Arkansas   1.289810 2.826039 2.717758 0.000000   3.763641 2.831051
## California 3.263110 3.012541 1.310484 3.763641   0.000000 1.287619
## Colorado   2.651067 2.326519 1.365031 2.831051   1.287619 0.000000

En esta tabla se muestra la distancia que posee cada dato con los demás

Dendograma

Los dendogramas corresponden a la representación gráfica de los árboles jerárquicos generados por la función hclust(). Estos son producidos por la función plot(res.hc)

res.hc <- hclust(d = res.dist, method = "ward.D2")
library("factoextra")
fviz_dend(res.hc, cex = 0.5)

## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

En la figura se muestra un dendograma correspondiente a un objeto, aquí se muestra cómo los datos similares entre sí están agrupados

Verificando el arbol de agrupaciones

Luego de enlazar los objetos en un árbol jerárquico es necesario evaluar las distancias con respecto a las originales. Para medir la correlación entre ambas distancias se debe tener en cuenta que entr más cercano sea esta valor a 1, es más precisa la solución de las agrupaciones. Para este proceso se utiliza la función cophenetic().

res.coph <- cophenetic(res.hc)
cor(res.dist, res.coph)

## [1] 0.6975266

El algoritmo anteior busca computar la distancia cophenetica y correlacionarla con la distancia original, el valor dado es esa relación.

res.hc2 <- hclust(res.dist, method = "average")
cor(res.dist, cophenetic(res.hc2))

## [1] 0.7180382

En este algoritmo se ejecutó la función hclust() con el método de enlace promedio. El coeficiente de correlación muestra que usando diferentes métodos de enlace se genera un árbol con una representación mejor de las distancias originales.

Cortar el dendrograma en diferentes grupos

Se puede cortar el árbol jerárquico a cualquier altura para dividir las agrupaciones. Para esto se utiliza la función cutree() . Esta se encarga dedevolver un vectos el cual contiene el número de agrupación de cada observación

grp <- cutree(res.hc, k = 4)
head(grp, n = 4)

##  Alabama   Alaska  Arizona Arkansas 
##        1        2        2        3

En este algoritmo se “cortó” el arbol en 4 grupos

table(grp)

## grp
##  1  2  3  4 
##  7 12 19 12

En este código se muestra el número de miembros en cada agrupación

rownames(df)[grp == 1]

## [1] "Alabama"        "Georgia"        "Louisiana"      "Mississippi"   
## [5] "North Carolina" "South Carolina" "Tennessee"

En este código se muestran los nombres de los miembros de la agrupación 1

El resultado de los cortessepede ver por medio de:

fviz_dend(res.hc, k = 4,
          cex = 0.5,
          k_colors = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
          color_labels_by_k = TRUE,
          rect = TRUE)

Además de esto, usando la función fviz_cluster() se puede visualizar el resultado en un gráfico de dispersión.

fviz_cluster(list(data = df, cluster = grp),
             palette= c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
             ellipse.type = "convex",
             repel = TRUE,
             show.clust.cent = FALSE)

Agrupación en paquetes de R

Las funciones agnes() y diana() para computar aglomeraciones en agrupaciones permiten analizar estos datos.

library("cluster")
res.agnes <- agnes(x = USArrests,
                   stand = TRUE,
                   metric = "euclidean",
                   method = "ward")

res.diana <- diana(x = USArrests,
                   stand = TRUE,
                   metric = "euclidean")
fviz_dend(res.agnes, cex = 0.6, k = 4)

Capítulo 8: Comparando dendrogramas

Para la comparación de dendogramas será usado el paquete dendextend(). Este paquete provee la función tanglegram () y cor.dendlist () los cuales se encargan de comparar visualmente los dendogramas y computar una matriz de correlación, respectivamente.

Preparación de datos

set.seed(123)
ss <- sample(1:50, 10)
df <- df[ss,]

Comparando dendrogramas

Para comparar las dendogramas se debe iniciar creando dos de estos. Esto s ehaace al computar dos agrupaciones jerárquicas usando métodos de enlace. Posteriormente se transforman los resultados en dendogramas y se genera una lista.

library(dendextend)

## 
## ---------------------
## Welcome to dendextend version 1.17.1
## Type citation('dendextend') for how to cite the package.
## 
## Type browseVignettes(package = 'dendextend') for the package vignette.
## The github page is: https://github.com/talgalili/dendextend/
## 
## Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
## You may ask questions at stackoverflow, use the r and dendextend tags: 
##   https://stackoverflow.com/questions/tagged/dendextend
## 
##  To suppress this message use:  suppressPackageStartupMessages(library(dendextend))
## ---------------------

## 
## Attaching package: 'dendextend'

## The following object is masked from 'package:stats':
## 
##     cutree

res.dist <- dist(df, method = "euclidean")

hc1 <- hclust(res.dist, method = "average")
hc2 <- hclust(res.dist, method = "ward.D2")

dend1 <- as.dendrogram (hc1)
dend2 <- as.dendrogram(hc2)

dend_list <-dendlist(dend1, dend2)

Comparación visual de dos dendrogramas

Para comparar visualmente ambos dendrogramas se utiliza la función tanglegram() el cual muestra los dos dendrogramas uno al lado del otro.

tanglegram(dend1, dend2,
           highlight_distinct_edges= FALSE,
           common_subtrees_color_lines = FALSE,
           common_subtrees_color_branches= TRUE,
           main = paste("entanglement =", round(entanglement(dend_list), 2)))

Matriz de correlación entre una lista de dendrogramas

Para realizar una matriz de correlación se utiliza la función cor.dendlist(). Esta función se encarga de correlacionar las listas por medio de un valor entre -1 y 1, siendo 0 un valor representante de la no disimilitud entre los árboles.

cor.dendlist(dend_list, method = "cophenetic")

##           [,1]      [,2]
## [1,] 1.0000000 0.9925544
## [2,] 0.9925544 1.0000000

cor.dendlist(dend_list, method = "baker")

##           [,1]      [,2]
## [1,] 1.0000000 0.9895528
## [2,] 0.9895528 1.0000000

cor_cophenetic(dend1, dend2)

## [1] 0.9925544

cor_bakers_gamma(dend1, dend2)

## [1] 0.9895528

dend1 <- df %>% dist %>% hclust("complete") %>% as.dendrogram()

dend2 <- df %>% dist %>% hclust("single") %>% as.dendrogram()

dend3 <- df %>% dist %>% hclust("average") %>% as.dendrogram()

dend4 <- df %>% dist %>% hclust("centroid") %>% as.dendrogram()

dend_list <- dendlist("Complete" = dend1, "Single" = dend2, "Average" = dend3, "Centroid" = dend4)

cors <- cor.dendlist(dend_list)

round(cors, 2)

##          Complete Single Average Centroid
## Complete     1.00   0.46    0.45     0.30
## Single       0.46   1.00    0.23     0.17
## Average      0.45   0.23    1.00     0.31
## Centroid     0.30   0.17    0.31     1.00

Capítulo 9: visualizando dendrogramas

En la siguiente sección se utilizarán las siguientes funciones:

fviz_dend() parageenrar un dendrograma basado en ggplot2
el paquete dendextend el cuál se encarga de manipular los dendrogramas

data("USArrests")
dd <- dist(scale(USArrests), method = "euclidean")
hc <- hclust(dd, method = "ward.D2")

Visualizando dendrogramas

Para crear un dendrograma básico:

library(factoextra)
fviz_dend(hc, cex = 0.5)

Para cambiar los títulos:

fviz_dend(hc, cex = 0.5,
          main = "Dendrogram - ward.D2",
          xlab = "Objects", ylab = "Distance", sub = "")

Para generar un dendrograma horizontal:

fviz_dend(hc, cex = 0.5, horiz = TRUE)

Para cambiar el tema:

fviz_dend(hc, k = 4,
          cex = 0.5,
          k_colors = c ("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
          color_labels_by_k = TRUE,
          ggtheme = theme_gray())

Para cambiar los colores de los grupos usando una paleta:

fviz_dend(hc, cex = 0.5, k = 4,
          k_colors = "jco")

Para hacer un dendrograma horizontal con los rectangulos alrededor de las agrupaciones:

fviz_dend(hc, k = 4, cex = 0.4, horiz = TRUE, k_colors = "jco", rect = TRUE, rect_border = "jco", rect_fill = TRUE)

Para hacer un dendrograma circular:

fviz_dend(hc, cex = 0.5, k = 4, k_colors = "jco", type = "circular")

Para hacer un dendrograma con forma de “arbol filogenetico”

require("igraph")

## Loading required package: igraph

## 
## Attaching package: 'igraph'

## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum

## The following object is masked from 'package:base':
## 
##     union

fviz_dend(hc, k = 4, k_colors = "jco", type = "phylogenic", repel = TRUE)

require("igraph")
fviz_dend(hc, k = 4, k_colors = "jco", type = "phylogenic", repel = TRUE, phylo_layout = "layout.gem")

Caso de dendrogramas con gran cantidad de datos

Cuando se trabaja con un dendrograma de gran tamaño, muchas veces se quiere realizar un zoom en alguna parte específica de este, para ello son las siguientes erramientas.

Zoom en un dendrograma

fviz_dend(hc, xlim = c(1, 20), ylim = c(1, 8))

Trazar arboles a partir de dendrogramas

Para realizar esto se debe:

Crear el dendograma entero y guardar el resultado en un dend_plot
Usar la función base cut.dendrogram() para recortar la altura del dendrograma en múltiples arboles
Visualizar los árboles usando fviz_dend()

dend_plot <- fviz_dend(hc, k = 4, cex = 0.5, k_colors = "jco")
dend_data <- attr(dend_plot, "dendrogram")

dend_cuts <- cut(dend_data, h = 10)

fviz_dend(dend_cuts$upper)

## Warning in min(-diff(our_dend_heights)): no non-missing arguments to min;
## returning Inf

print(dend_plot)

fviz_dend(dend_cuts$lower[[1]], main = "SUbtree 1")

fviz_dend(dend_cuts$lower[[2]], main = "SUbtree 2")

fviz_dend(dend_cuts$lower[[2]], type = "circular")

Guardar dendrogramas en una página de PDF

Se puede guardar un dendrograma grande en una página de PDF sin perder la resolción, esto por medio de:

pdf("dendrogram.pdf", width=30, height = 15)

p <- fviz_dend(hc, k=4, cex = 1, k_colors = "jco")
print(p)
dev.off()

## png 
##   2

Manipular dendrogramas usando dendextend

El paquete dendextend provee funciones para cambiar la apariencia de los dendrogramas.

El código estándar para crear dendrogramas:

data <- scale(USArrests)
dist.res <- dist(data)
hc <- hclust(dist.res, method = "ward.D2")
dend <- as.dendrogram(hc)
plot(dend)

Código para crear un dendrograma con una operadora de encadenamiento:

library(dendextend)
dend <- USArrests [1:15, ] %>%
  scale %>%
  dist %>%
  hclust(method = "ward.D2") %>%
  as.dendrogram
plot(dend)

Función para personalizar los dendrogramas:

library(dendextend)

mycols <- c("#2E9FDF","#00AFBB", "#E7B800", "#FC4E07")
dend <- as.dendrogram(hc) %>%
  set("branches_lwd", 1) %>%
  set("branches_k_color", mycols, k = 4) %>%
  set("labels_colors", mycols, k =4) %>%
  set("labels_cex", 0.5)

fviz_dend(dend)

Capítulo 10: Mapa de calor: estatico e interactivo

Un mapa del calor es una manera de visualizar agrupaciones jerárquicas. Aquí los datos son transformados en colores.

Base de mapas de calor para R: heatmap()

df <- scale(mtcars)

heatmap(df, scale = "none")

También es posible especificar los colores de la paleta.

Usando colores personalizados:

col <- colorRampPalette(c("red", "white", "blue")) (256)

Usando la paleta de RColorBrewer

library("RColorBrewer")
col <- colorRampPalette(brewer.pal(10, "RdYlBu"))(256)

Ejemplificando:

library("RColorBrewer")
col <- colorRampPalette(brewer.pal(10, "RdYlBu"))(256)
heatmap(df, scale = "none", col = col, RowSideColors = rep(c("blue", "pink"), each = 16), ColSideColors = c(rep("purple", 5), rep("orange", 6)))

Mapas decalor mejorados: heatmap.2()

Para esta sección se utiliza la función heatmap.2()

library("gplots")

## 
## Attaching package: 'gplots'

## The following object is masked from 'package:stats':
## 
##     lowess

heatmap.2(df, scale = "none", col = bluered (100), trace = "none", density.info = "none")

Mapas de calor lindos : pheatmap()

Para realizar esto se debe instalar el paquete pheatmap

library("pheatmap")
pheatmap(df, cutree_rows = 4)

Mapas de calor interactivos: d3heatmap()

Para esta función se debe instalar el paquete d3heatmap.

Debido a que este paquete no existe en la versión R 4.3.2 (La que ha sido usada para este proyecto) no es posible aplicar el código, sin embargo este quedará escrito.

```{install.packages(“d3heatmap”)}

{library("d3heatmap")}
d3heatmap(scale(mtcars), colors = "RdYlBu", k_row = 4, k_col = 2)

Mejorando mapas de calor usando dendextend

El paquete dendextend es usado para mejorar funciones de otros paquetes.

Se comenzará definiendo el orden y la apariencia de las filas y columnas usando este paquete.

library(dendextend)

Rowv <- mtcars %>% scale %>% dist %>% hclust %>% as.dendrogram %>%
  set("branches_k_color", k = 3) %>%
  set("branches_lwd", 1.2)%>%
  ladderize()

Colv <- mtcars %>% scale %>% t %>% dist %>% hclust %>% as.dendrogram %>%
  set("branches_k_color", k = 2, value = c("orange", "blue")) %>%
  set("branches_lwd", 1.2) %>%
  ladderize()

El código posterior puede ser usado en las siguientes fuciones:

heatmap(scale(mtcars), Rowv = Rowv, Colv = Colv, scale = "none")

library(gplots)
heatmap.2(scale (mtcars), scale = "none", col = bluered(100), 
          Rowv = Rowv, Colv = Colv,
          trace = "none", density.info = "none")

{library("d3heatmap")} d3heatmap(scale(mtcars), colors = "RdBu", Rowv = Rowv, Colv = Colv)

Mapas de calor complejos

El paquete ComplexHeatmap provee una solución para organizar, y anotar múltiples mapas de calor, permite visualizar la asociación de datos entre diferentes fuentes u orígenes.

Este es instalado:

library(devtools)

## Loading required package: usethis

install_github("jokergoo/ComplexHeatmap")

## Skipping install of 'ComplexHeatmap' from a github remote, the SHA1 (7d95ca5c) has not changed since last install.
##   Use `force = TRUE` to force installation

Mapa de calor simple

El código para generar un mapa de calor simple es el siguiente:

library(ComplexHeatmap)

## Loading required package: grid

## ========================================
## ComplexHeatmap version 2.15.4
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
## 
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
##     genomic data. Bioinformatics 2016.
## 
## 
## The new InteractiveComplexHeatmap package can directly export static 
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================
## ! pheatmap() has been masked by ComplexHeatmap::pheatmap(). Most of the arguments
##    in the original pheatmap() are identically supported in the new function. You 
##    can still use the original function by explicitly calling pheatmap::pheatmap().

## 
## Attaching package: 'ComplexHeatmap'

## The following object is masked from 'package:pheatmap':
## 
##     pheatmap

Heatmap(df, name = "mtcars",
        column_title = "Variables", row_title = "Samples",
        row_names_gp = gpar(fontsize = 7))

Para especificar algunos colores personalizados se utiliza:

library(circlize)

## ========================================
## circlize version 0.4.15
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: https://jokergoo.github.io/circlize_book/book/
## 
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
##   in R. Bioinformatics 2014.
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(circlize))
## ========================================

## 
## Attaching package: 'circlize'

## The following object is masked from 'package:igraph':
## 
##     degree

mycols <- colorRamp2(breaks = c(-2, 0, 2),
                     colors = c("green", "white", "red"))
Heatmap(df, name = "mtcars", col = mycols)

Para usar la paleta de colores de RColorBrewer:

library("circlize")
library("RColorBrewer")
Heatmap(df, name = "mtcars", col = colorRamp2(c(-2, 0, 2), brewer.pal(n = 3, name = "RdBu")))

También se puede personalizar la apariencia de los dendogramas usando color_branches():

library(dendextend)
row_dend = hclust(dist(df))
col_dend = hclust(dist(t(df)))
Heatmap(df, name = "mtcars", 
        row_names_gp = gpar (fontsize = 6.5), cluster_rows = color_branches(row_dend, k = 4),
        cluster_columns = color_branches(col_dend, k = 2))

Dividiendo mapas de calor por filas

Es posible dividir el mapa de calor usando el algoritmo de k-means o variable de agrupación.

Para dividir el dendrograma con K-means:

set.seed(2)
Heatmap(df, name = "mtcars", k = 2)

Para dividir por medio de variable de agrupación:

Heatmap(df, name = "mtcars", split = mtcars$cyl,
        row_names_gp = gpar(fontsize = 7))

Para dividir combinando múltiples variables:

Heatmap(df, name = "mtcars",
        split = data.frame(cyl = mtcars$cyl, am = mtcars$am))

Anotación de mapas de calor

HeatmapAnnotation es usado para definir la anotación de filas o columnas.

HeatmapAnnotation(df, name, col, show_legend)

En esta estructura se evidencian;

df: data.frame con nombre de columnas
name: nombre de la anotación del mapa de calor
col: lista de colores

Anotación compleja

Las anotaciones complejas consisten en combinar mapas de calor con ptrps tipos de gráficos básicos para mostrar la distribución de los datos.

.hist = anno_histogram(df, gp = gpar(fill = "lightblue"))
.density = anno_density(df, type = "line", gp = gpar(col = "blue"))
ha_mix_top = HeatmapAnnotation(
  hist = .hist, density = .density,
  height = unit(3.8, "cm")
  )

.violin = anno_density(df, type = "violin", 
                       gp = gpar(fill = "lightblue"), which = "row")
.boxplot = anno_boxplot(df, which = "row")
ha_mix_right = HeatmapAnnotation(violin = .violin, bxplt = .boxplot,
                              which = "row", width = unit(4, "cm"))


Heatmap(df, name = "mtcars", 
        column_names_gp = gpar(fontsize = 8),
        top_annotation = ha_mix_top) + ha_mix_right

Combinando múltiples mapas de calor

ht1 = Heatmap(df, name = "ht1", km = 2,
              column_names_gp = gpar(fontsize = 9))

ht2 = Heatmap(df, name = "ht2",
              col = circlize::colorRamp2(c(-2, 0, 2), c("green", "white", "red")),
              column_names_gp = gpar(fontsize = 9))
ht1 + ht2

Aplicación para matriz de expresión de genes

En datos de expresión de genes, las filas spn genes y las columnas son muestras.

expr <- readRDS(paste0(system.file(package = "ComplexHeatmap"),
                       "/extdata/gene_expression.rds"))

mat <- as.matrix(expr[, grep("cell", colnames(expr))])

type <- gsub("s\\d+_", "", colnames(mat))
ha = HeatmapAnnotation(df = data.frame(type = type))


Heatmap(mat, name = "expression", km = 5, top_annotation = ha,
    show_row_names = FALSE, show_column_names = FALSE) +
Heatmap(expr$length, name = "length", width = unit(5, "mm"),
    col = circlize::colorRamp2(c(0, 100000), c("white", "orange"))) +
Heatmap(expr$type, name = "type", width = unit(5, "mm")) +
Heatmap(expr$chr, name = "chr", width = unit(5, "mm"),
    col = circlize::rand_color(length(unique(expr$chr))))

## There are 23 unique colors in the vector `col` and 23 unique values in
## `matrix`. `Heatmap()` will treat it as an exact discrete one-to-one
## mapping. If this is not what you want, slightly change the number of
## colors, e.g. by adding one more color or removing a color.

Actividad_4_R_cap_4_10

Valentina Londoño

2023-11-19

Capítulo 4: K-Means Clustering

Computar K-means clustering en R

Datos

Paquetes y funciones requeridos

Estimación del numéro optimo de grupos

Computando k-means clustering

Acceso a los resultados de la función kmeans ()

Visualización de grupos k-means

Capítulo 5: K-medoids

Computar PAM en R

Datos

Paquetes y funciones necesarios

Estimación del número óptimo de agrupaciones

Computar agrupaciones PAM

Acceso a los resultados de la función pam()

Visualización de las agrupaciones PAM

Capítulo 6: CLARA- Agrupando grandes aplicaciones

Computando CLARA en R

Datos y preparación

Paquetes y funciones requeridos

Computando CLARA

Visualizando agrupaciones CLARA

Capítulo 7: agrupación aglomerativa

Pasos para una agrupación aglomerativa herarquica

Estructura y preparación de los datos

Medidas de similaridad

Dendograma

Verificando el arbol de agrupaciones

Cortar el dendrograma en diferentes grupos

Agrupación en paquetes de R

Capítulo 8: Comparando dendrogramas

Preparación de datos

Comparando dendrogramas

Comparación visual de dos dendrogramas

Matriz de correlación entre una lista de dendrogramas

Capítulo 9: visualizando dendrogramas

Visualizando dendrogramas

Caso de dendrogramas con gran cantidad de datos

Zoom en un dendrograma

Trazar arboles a partir de dendrogramas

Guardar dendrogramas en una página de PDF

Manipular dendrogramas usando dendextend

Capítulo 10: Mapa de calor: estatico e interactivo

Base de mapas de calor para R: heatmap()

Mapas decalor mejorados: heatmap.2()

Mapas de calor lindos : pheatmap()

Mapas de calor interactivos: d3heatmap()

Mejorando mapas de calor usando dendextend

Mapas de calor complejos

Mapa de calor simple

Dividiendo mapas de calor por filas

Anotación de mapas de calor

Anotación compleja

Combinando múltiples mapas de calor

Aplicación para matriz de expresión de genes