Conjunto de datos

Consideremos el conjunto de datos \(\texttt{planets}\) del paquete \(\texttt{HSAUR2}\) que contiene una base de datos de 101 observaciones para distintas variables asociadas a datos de exoplanetas. El objetivo es encontrar un determinado número de categorías o clústers donde agrupar a los exoplanetas.

library(HSAUR2)
## Cargando paquete requerido: tools
## 
## Adjuntando el paquete: 'HSAUR2'
## The following object is masked from 'package:tidyr':
## 
##     household
data <- planets |> 
  glimpse()
## Rows: 101
## Columns: 3
## $ mass   <dbl> 0.120, 0.197, 0.210, 0.220, 0.230, 0.250, 0.340, 0.400, 0.420, …
## $ period <dbl> 4.95000, 3.97100, 44.28000, 75.80000, 6.40300, 3.02400, 2.98500…
## $ eccen  <dbl> 0.000, 0.000, 0.340, 0.280, 0.080, 0.020, 0.080, 0.498, 0.000, …

El primer paso que vamos a dar es visualizar los datos, en esta ocasión como disponemos de 3 variables, vamos a representar tanto los pares de variables, como las 3 a la vez con un gráfico situándo una de ellas, por ejemplo, con el color.

plot.mp <- data |> 
  ggplot() +
  aes(x = mass,
      y = period) +
  geom_point() + 
  theme_bw()

plot.me <- data |> 
  ggplot() +
  aes(x = mass,
      y = eccen) +
  geom_point() + 
  theme_bw()

plot.pe <- data |> 
  ggplot() +
  aes(x = period,
      y = eccen) +
  geom_point() + 
  theme_bw()

library(ggpubr)
ggarrange(plot.mp,
          plot.me,
          plot.pe,
          ncol=1,nrow=3)

Podemos observar que hay una gran cantidad de valores pequeños agrupados alrededor del 0, es por ello que sería ideal aplicar una determinada transformación para evitar este problema. Como se trata de variables positivas, podríamos aplicar la transformación logarítmica (debido también a que presentan largas colas).

data_log <- data |> 
  transmute(
    mass_log = log(mass+1), 
    period_log = log(period+1),
    eccen_log = log(eccen+1)) |> #Le sumaremos 1 para evitar la presencia de 0
  glimpse()
## Rows: 101
## Columns: 3
## $ mass_log   <dbl> 0.1133287, 0.1798184, 0.1906204, 0.1988509, 0.2070142, 0.22…
## $ period_log <dbl> 1.783391, 1.603621, 3.812865, 4.341205, 2.001885, 1.392276,…
## $ eccen_log  <dbl> 0.000000000, 0.000000000, 0.292669614, 0.246860078, 0.07696…
plot.mp <- data_log |> 
  ggplot() +
  aes(x = mass_log,
      y = period_log) +
  geom_point() + 
  labs(x = "log(mass)",
       y = "log(period)") +
  theme_bw()

plot.me <- data_log |> 
  ggplot() +
  aes(x = mass_log,
      y = eccen_log) +
  geom_point() + 
  labs(x = "log(mass)",
       y = "log(eccen)") +
  theme_bw()

plot.pe <- data_log |> 
  ggplot() +
  aes(x = period_log,
      y = eccen_log) +
  geom_point() + 
  labs(x = "log(period)",
       y = "log(eccen)") +
  theme_bw()

library(ggpubr)
ggarrange(plot.mp,
          plot.me,
          plot.pe,
          ncol=1,nrow=3)

Con esta transformación podemos observar que efectivamente existen relaciones entre las variables, por lo que son útiles para poder construir clústers.

data_log |> 
  ggplot() + 
  aes(x = mass_log,
      y = period_log,
      color = eccen_log) + 
  labs(x="log(mass)",y="log(period)") + 
  scale_color_continuous(name="log(eccen)") + 
  geom_point() + 
  theme_bw()

En este último gráfico sí que podríamos ver que se podrían crear diferentes clúster que nos permitan agrupar las observaciones. Representemos ahora un gráfico 3D, para ello vamos a emplear un nuevo paquete de gráficos, \(\texttt{plotly}\) que permite crear gráficos interactivos.

library(plotly)
## 
## Adjuntando el paquete: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
data_log |> 
  plot_ly(x = ~ mass_log,
          y = ~ period_log,
          z = ~ eccen_log,
          color = ~ eccen_log,
          type = "scatter3d",
          mode = "markers"
  ) |> 
  layout(scene = list(xaxis = list(title = "log(mass)"),
                      yaxis = list(title = "log(period)"),
                      zaxis = list(title = "log(eccen)"))
  )

Con este gráfico sí que se observa que hay al menos 3 clústers en los que agrupar las observaciones.

Agrupamiento Jerárquico: Dendograma

En el Agrupamiento Jerárquico no es necesario predefinir el número de clústers, pero si es importante estandarizar previamente las variables.

data_log <- data_log |> 
  mutate(
    across(c(mass_log,
             period_log,
             eccen_log), 
           ~ as.numeric(scale(.))
    )
  ) |> 
  glimpse()
## Rows: 101
## Columns: 3
## $ mass_log   <dbl> -1.5124356, -1.4190589, -1.4038889, -1.3923301, -1.3808658,…
## $ period_log <dbl> -1.66550283, -1.75241650, -0.68431095, -0.42887428, -1.5598…
## $ eccen_log  <dbl> -1.45405728, -1.45405728, 0.35679448, 0.07335447, -0.977871…

Recordemos que este algoritmo se basa en encontrar jerarquías en los datos a partir de generar grupos basados en la cercanía de ellos, esto se puede visualizar a través del dendograma. Para ello, vamos a utilizar la función \(\texttt{hclust()}\) y la función \(\texttt{ggdendogram()}\) del paquete \(\texttt{ggdendro}\) para visualizar el dendograma.

Medidas de Disparidad entre Observaciones: Distancias

Pero para ello es necesario construir previamente la matriz de disparidad a partir de la distancia entre observaciones elegida, nosotros hemos visto la distancia euclídea y la distancia Manhattan. El cálculo de ambas se realiza mediante la función $.

dist_eucli <- dist(data_log,method="euclidean")

dist_manha <- dist(data_log,method="manhattan")

Medidas de Disparidad entre Clústers: Enlaces

A continuación, se debe decidir cuál es la distancia que mide la disparidad entre grupos, esto es, se debe elegir entre enlace único, enlace completo o enlace medio. Para ello se emplea el argumento \(\texttt{method}\) con las opciones:

  • \(\texttt{single}\): para enlace único (usa la distancia mínima entre los puntos del clúster)
  • \(\texttt{complete}\): para enlace completo (usa la distancia máxima entre los puntos del clúster)
  • \(\texttt{average}\): para enlace medio (emplea la distancia media entre clústers)
  • \(\texttt{ward.D}\): este enlace se construye minimizando la varianza dentro de los clústers.

Vamos a probar con cada uno de ellos, empleando la distancia euclídea.

library(ggdendro)
model_dendo.single <- dist_eucli |> 
  hclust(method = "single")

model_dendo.single |> 
  ggdendrogram(
    rotate = FALSE,
    labels = FALSE,
    theme_dendro = TRUE
  ) +
  labs(title = "Dendograma con enlace único")

model_dendo.complete <- dist_eucli |> 
  hclust(method = "complete")

model_dendo.complete |> 
  ggdendrogram(
    rotate = FALSE,
    labels = FALSE,
    theme_dendro = TRUE
  ) +
  labs(title = "Dendograma con enlace completo")

model_dendo.average <- dist_eucli |> 
  hclust(method = "average")

model_dendo.average |> 
  ggdendrogram(
    rotate = FALSE,
    labels = FALSE,
    theme_dendro = TRUE
  ) +
  labs(title = "Dendograma con enlace medio")

model_dendo.ward <- dist_eucli |> 
  hclust(method = "ward.D")

model_dendo.ward |> 
  ggdendrogram(
    rotate = FALSE,
    labels = FALSE,
    theme_dendro = TRUE
  ) +
  labs(title = "Dendograma minimizando la varianza")

Interpretación

  • En el eje horizontal se tiene cada uno de los datos que componen el conjunto de entrenamiento (el argumento \(\texttt{labels}\) permite identificar cada observación.)

  • En el eje vertical se represnta la distancia (euclídea o Manhattan) que existe entre cada clúster a medida que éstos se van jerarquizando. De modo que a baja altura indica que las observaciones son similares y una alta representa disparidad.

  • Cada línea vertical representa un agrupamiento que unifica los nodos que forman la rama. A la vista del gráfico, podría decirse que se pueden agrupar las observaciones en 3 clústers.

Representemos con colores los clústers para mejorar la visualización, para ello vamos a indicar a qué distancia debemos hacer la agrupación. Para ello vamos a emplear funciones del paquete \(\texttt{dendextend}\) que permite extender las funciones para construir dendogramas.

library(dendextend)
## 
## ---------------------
## Welcome to dendextend version 1.19.0
## 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))
## ---------------------
## 
## Adjuntando el paquete: 'dendextend'
## The following object is masked from 'package:ggdendro':
## 
##     theme_dendro
## The following object is masked from 'package:ggpubr':
## 
##     rotate
## The following object is masked from 'package:stats':
## 
##     cutree
model_dendo.ward |> 
  as.dendrogram() |> 
  color_branches(k = 3) |>  
  color_labels(k = 3) |>  #colorea según el clúster
  set("labels_cex",0.4) |> 
  plot() |>  
  abline(h = 15, lty=2)

Búsqueda óptima del número de clústers \(k\)

En el dendograma podemos observar qué observaciones pertenecen a cada clúster, considerando la agrupación en 3 grupos. Sin embargo, podemos comprobar un valor óptimo para el número de clústers \(k\) a partir de la función \(\texttt{fviz_nbclust()}\) del paquete \(\texttt{factoextra}\).

library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
data_log |> 
  fviz_nbclust(FUNcluster = hcut, 
               method = "silhouette")

Podemos observar en el gráfico que el valor que mejor promedio presenta es considerando \(k = 2\) clústers. Equivalente a los dendogramas anteriores, este paquete también contiene la función \(\texttt{fviz_dend()}\) que representa el dendograma junto con las observaciones coloreadas según su clúster.

model_dendo.ward |> 
  fviz_dend(k = 2,
            rect = TRUE,
            horiz = FALSE,
            rect_border = "gray",
            rect_fill = FALSE,
            cex = 0.4,
            lwd = 0.2,
            k_colors = c("red","blue"),
            ggtheme = theme_bw()
  )
## 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.

## Dendogramas: otras representaciones

También es posible cambiar el tipo de dendograma, por ejemplo uno curioso es verlo como una red o circular con el argumento \(\texttt{type}\).

model_dendo.ward |> 
  fviz_dend(k = 2,
            rect = TRUE,
            horiz = FALSE,
            rect_border = "gray",
            rect_fill = FALSE,
            cex = 0.4,
            lwd = 0.2,
            k_colors = c("red","blue"),
            ggtheme = theme_bw(),
            type = "phylogenic"
  )

model_dendo.ward |> 
  fviz_dend(k = 2,
            rect = TRUE,
            horiz = FALSE,
            rect_border = "gray",
            rect_fill = FALSE,
            cex = 0.4,
            lwd = 0.2,
            k_colors = c("red","blue"),
            ggtheme = theme_bw(),
            type = "circular"
  )

## Poda del Dendograma y Construcción de Clústers

Finalmente, una vez decidido cuántos clúster crear, se puede “podar” el dendograma para situar las observaciones en cada uno de los clústers y añadirlos al conjunto de datos. Para ello se emplea la función \(\texttt{cutree()}\) indicando el número de clústers \(k\) en los que agrupar las observaciones.

data_log <- data_log |> 
  mutate(
    cluster_jerar = factor(model_dendo.ward |> 
                             cutree(k = 2))
  ) |> 
  glimpse()
## Rows: 101
## Columns: 4
## $ mass_log      <dbl> -1.5124356, -1.4190589, -1.4038889, -1.3923301, -1.38086…
## $ period_log    <dbl> -1.66550283, -1.75241650, -0.68431095, -0.42887428, -1.5…
## $ eccen_log     <dbl> -1.45405728, -1.45405728, 0.35679448, 0.07335447, -0.977…
## $ cluster_jerar <fct> 1, 1, 2, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 2, 1, 1, 2, 2, 2,…

Podemos volver a representar las observaciones pero ya clasificadas en los clústers.

library(plotly)

data_log |> 
  plot_ly(x = ~ mass_log,
          y = ~ period_log,
          z = ~ eccen_log,
          color = ~ cluster_jerar,
          type = "scatter3d",
          mode = "markers",
          colors = c("red","blue")
  ) |> 
  layout(scene = list(xaxis = list(title = "log(mass)"),
                      yaxis = list(title = "log(period)"),
                      zaxis = list(title = "log(eccen)"))
  )

Método K-medias

Búsqueda óptima del número de clústers \(K\)

Se trata de otra técnica de agrupación de observaciones en \(K\) grupos distintos (clústers). A diferencia del agrupamiento jerárquico, K-medias necesita conocer previamente el número de clúster a considerar. Este método construye los clúster minimizando la variación total dentro del clúster, por ello vamos a emplear la función \(\texttt{fviz_nbclust()}\) del paquete \(\texttt{factoextra}\) empleando los argumentos \(\texttt{FUNcluster}\) y \(\texttt{method}\) con los valores \(\texttt{kmeans}\) y \(\texttt{wss}\), respectivamente.

library(factoextra)

data_log |> 
  select(-cluster_jerar) |> #debemos eliminar esta variable
  fviz_nbclust(FUNcluster = kmeans, 
               method = "wss")

Como se observa, la varianza total disminuye a medida que aumenta el número de clúster. Esta curva parece mostrar un ``codo’’ en \(K = 3\). Una vez decidido el número de clúster se puede proceder a ejecutar la función \(\texttt{kmeans()}\) que implementa la técnica de K-medias. En esta función es necesario indicar el número de clústers con el argumento \(\texttt{center}\), así como el número de conjuntos iniciales aleatorias (\(\texttt{nstart}\)).

model_kmeans <- data_log |> 
  select(-cluster_jerar) |> 
  kmeans(centers = 3,
         nstart = 10)

model_kmeans
## K-means clustering with 3 clusters of sizes 42, 31, 28
## 
## Cluster means:
##     mass_log period_log  eccen_log
## 1 -0.2664146  0.5391742 -0.1177906
## 2  1.1094462  0.5019882  0.9471295
## 3 -0.8286936 -1.3645339 -0.8719217
## 
## Clustering vector:
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   1   3   1   1 
##  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
##   1   3   1   3   3   3   1   1   1   1   1   1   3   1   1   1   1   1   1   1 
##  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 
##   3   1   2   1   1   1   1   1   2   1   1   3   3   2   1   1   1   1   1   1 
##  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80 
##   1   1   1   1   2   1   1   2   2   3   2   2   3   1   2   3   2   2   2   2 
##  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 
##   2   1   1   1   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2 
## 101 
##   2 
## 
## Within cluster sum of squares by cluster:
## [1] 45.04624 44.45944 28.29321
##  (between_SS / total_SS =  60.7 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"

En esta salida podemos comprobar cuál es la media de cada grupo (para las 3 variables), así como las etiquetas de los clústers. Éstas pueden recuperarse a través \(\texttt{model.kmeans\$cluster}\).

table(model_kmeans$cluster)
## 
##  1  2  3 
## 42 31 28

Si mostramos las frecuencias absolutas de los clústers, podemos observar que hay 42 observaciones en el primer clúster, 31 en el segundo y 28 en el tercero.

Construcción de Clústers

A continuación, vamos a construir los clústers y visualizarlos.

data_log <- data_log |> 
  mutate(
    cluster_kmeans = factor(model_kmeans$cluster)
    ) |> 
  glimpse()
## Rows: 101
## Columns: 5
## $ mass_log       <dbl> -1.5124356, -1.4190589, -1.4038889, -1.3923301, -1.3808…
## $ period_log     <dbl> -1.66550283, -1.75241650, -0.68431095, -0.42887428, -1.…
## $ eccen_log      <dbl> -1.45405728, -1.45405728, 0.35679448, 0.07335447, -0.97…
## $ cluster_jerar  <fct> 1, 1, 2, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 2, 1, 1, 2, 2, 2…
## $ cluster_kmeans <fct> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 1, 3, 1…
library(plotly)

data_log |> 
  plot_ly(x = ~ mass_log,
          y = ~ period_log,
          z = ~ eccen_log,
          color = ~ cluster_kmeans,
          type = "scatter3d",
          mode = "markers",
          colors = c("red","blue","green")
  ) |> 
  layout(scene = list(xaxis = list(title = "log(mass)"),
                      yaxis = list(title = "log(period)"),
                      zaxis = list(title = "log(eccen)")))

Con los datos transformados es difícil tomar decisiones, por lo que vamos a calcular la media para cada uno de los clústers utilizando las variables originales.

data |> 
  mutate(
    cluster_kmeans = data_log |> select(cluster_kmeans)
    ) |>  
  group_by(cluster_kmeans) |>  
  summarise(
    means = across(everything(),mean)
    ) |> 
  pull(means)
## # A tibble: 3 × 3
##    mass period eccen
##   <dbl>  <dbl> <dbl>
## 1  1.97  881.  0.250
## 2  7.25  963.  0.483
## 3  1.02   17.4 0.105

Otra manera de visualizar los clústers creados es a través de la función \(\texttt{fviz_cluster()}\) del paquete \(\texttt{factoextra}\) que permit visualizar los resultados obtenidos junto con los datos originales en 2 dimensiones. En el caso de que existan más de dos variables (por ejemplo, nuestro caso), emplea la técnica de Componentes Principales para reducir la dimensión a 2.

library(factoextra)

model_kmeans |> 
  fviz_cluster(
    data = data,
    geom = "text",
             palette = c("red","blue","green"),
             ellipse = TRUE,
             ellipse.type = "euclid",
             star.plot = TRUE,
             repel = TRUE,
             ggtheme = theme_bw()
  )

En el gráfico podemos observar que cuando se ha aplicado la técnica de Componentes Principales para reducir la dimensión del conjunto de datos, la primera dimensión recoge el 49.5% de la información, mientras que la segunda dimensión solo recoge el 28.6% de la información.

En cuanto al gráfico, se pueden ver los 3 clústers creados, junto con las observaciones que pertenecen a cada uno de ellos.