Nos disponemos a generar un juego de datos aleatorio con el objetivo de experimentar y así entender mejor el funcionamiento del algoritmo k-means.
En este ejemplo vamos a generar un conjunto de muestras aleatorias para posteriormente usar el algoritmo k-means para agruparlas. Se crearán las muestras alrededor de dos puntos concretos. Por lo tanto, lo lógico será agrupar en dos clústers.
Puesto que inicialmente, en un problema real, no se conoce cual es el número más idóneo de clústers k, vamos a probar primero con dos (el valor óptimo) y posteriormente con 4 y 8 clústers. Para evaluar la calidad de cada proceso de agrupación vamos a usar la silueta media. La silueta de cada muestra evalúa como de bien o mal está clasificada la muestra en el clúster al que ha sido asignada. Para ello se usa una fórmula que tiene en cuenta la distancia a las muestras de su clúster y la distancia a las muestras del clúster vecino más cercano.
A la hora de probar el código que se muestra, es importante tener en cuenta que las muestras se generan de forma aleatoria y también que el algoritmo k-means tiene una inicialización aleatoria. Por lo tanto, en cada ejecución se obtendrá unos resultados ligeramente diferentes.
Lo primero que hacemos es cargar la librería cluster que contiene las funciones que se necesitan
if (!require('cluster')) install.packages('cluster')
library(cluster)Generamos las muestras de forma aleatoria tomando como centro los puntos [0,0] y [5,5].
n <- 150 # número de muestras
p <- 2 # dimensión
sigma <- 1 # varianza de la distribución
mean1 <- 0 # centro del primer grupo
mean2 <- 5 # centro del segundo grupo
n1 <- round(n/2) # número de muestras del primer grupo
n2 <- round(n/2) # número de muestras del segundo grupo
x1 <- matrix(rnorm(n1*p,mean=mean1,sd=sigma),n1,p)
x2 <- matrix(rnorm(n2*p,mean=mean2,sd=sigma),n2,p)Juntamos todas las muestras generadas y las mostramos en una gráfica
x <- rbind(x1,x2)
plot (x, xlab="Grupo 1", ylab="Grupo 2")
Como se puede comprobar las muestras están claramente separadas en dos grupos, además vemos puntos extremos o outliers fuera de las 2 zonas más densamente pobladas de puntos.
Si se quiere complicar el problema se puede modificar los puntos centrales (mean1 y mean2) haciendo que estén más próximos y/o ampliar la varianza (sigma) para que las muestras estén más dispersas.
A continuación vamos a aplicar el algoritmo k-means con 2, 4 y 8 clústers
# Función para ajustar k-means y mostrar los resultados
ajustar_kmeans <- function(k, vdata) {
fit <- kmeans(vdata, centers = k)
return(fit$cluster)
}
# Ajustos amb k = 2, 4, 8
clusters_2 <- ajustar_kmeans(2, x)
clusters_4 <- ajustar_kmeans(4, x)
clusters_8 <- ajustar_kmeans(8, x)Las variables y_cluster2, y_cluster4 e y_cluster8 contienen para cada muestra el identificador del clúster a las que han sido asignadas. Por ejemplo, en el caso de los k=2 las muestras se han asignado al clúster 1 ó al 2
clusters_2 ## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [75] 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [112] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [149] 2 2Para visualizar los clústers podemos usar la función clusplot. Vemos la agrupación con 2 clústers y observemos como prácticamente no hay valores extremos y realmente los dos clústers generados son homogéneos.
clusplot(x, clusters_2, color=TRUE, shade=TRUE, labels=2, lines=0)
con 4 observamos como uno de los dos clúster lo ha dividido en 3.
clusplot(x, clusters_4, color=TRUE, shade=TRUE, labels=2, lines=0)
y con 8. El algoritmo obedece y nos genera 8 clústers aunque como se aprecia visualmente no tenga demasiada consistencia.
clusplot(x, clusters_8, color=TRUE, shade=TRUE, labels=2, lines=0)
También podemos visualizar el resultado del proceso de agrupamiento con el siguiente código para el caso de 2 clústers. El uso de colores facilita la identificación visual de clústers.
# Función para graficar los resultados de los clusters
graficar_clusters <- function(vdata, y_cluster, k) {
colors <- rainbow(k)
plot(vdata, col = colors[y_cluster], xlab = "Dimensión 1", ylab = "Dimensión 2",
main = paste("Clusters con k =", k))
legend("topright", legend = 1:k, fill = colors, title = "Clusters")
}
# Gráfico para k = 2
graficar_clusters(x, clusters_2, 2) 
para 4
# Gráfico para k = 4
graficar_clusters(x, clusters_4, 4) 
y para 8
# Gráfico para k = 8
graficar_clusters(x, clusters_8, 8) 
Ahora vamos a evaluar la calidad del proceso de agregación. Para ello usaremos la función silhouette que calcula la silueta de cada muestra
# Función para calcular y mostrar la silueta
calcular_silueta <- function(y_cluster, vdata) {
distances <- daisy(vdata)
silueta <- silhouette(y_cluster, distances)
mean_sil <- mean(silueta[, 3])
return(mean_sil)
}
sil_2 <- calcular_silueta(clusters_2, x)
sil_4 <- calcular_silueta(clusters_4, x)
sil_8 <- calcular_silueta(clusters_8, x)La función silhouette devuelve para cada muestra, el clúster dónde ha sido asignado, el clúster vecino y el valor de la silueta. Por lo tanto, calculando la media de la tercera columna podemos obtener una estimación de la calidad del agrupamiento
# Mostrar valores de siluetas medianas
cat("Silueta mediana para k=2:", sil_2, "\n")## Silueta mediana para k=2: 0.7319339cat("Silueta mediana para k=4:", sil_4, "\n")## Silueta mediana para k=4: 0.5505738cat("Silueta mediana para k=8:", sil_8, "\n")## Silueta mediana para k=8: 0.3340561Como se puede comprobar, agrupar con dos clúster es mejor que en 4 o en 8, lo cual es lógico teniendo en cuenta como se han generado los datos.
Una buena práctica para entender mejor el juego de datos, consiste en poner nombre a cada uno de los clústers identificados. Lo veremos más claramente en el siguiente ejemplo que parte de datos reales.
A continuación vamos a ver otro ejemplo de cómo se usan los modelos de agregación. Para ello usaremos el data set penguins contenido en el paquete R palmerpenguins. Esta base de datos se encuentra descrita en https://cran.r-project.org/web/packages/palmerpenguins/index.html y contiene mediciones de tamaño, observaciones de puestas y proporciones de isótopos sanguíneos de tres especies de pingüinos observadas en tres islas del archipiélago Palmer, en la Antártida, durante un período de estudio de tres años.
El propósito de aplicar un modelo de agregación sobre este juego de datos es el de identificar patrones, es decir, quisiéramos verificar si a partir de la caracterización de los pingüinos (mediciones diversas) es realmente posible establecer la especie a la que pertenecen.
Este dataset está previamente trabajado para que los datos estén limpios y sin errores. De no ser así antes de nada deberíamos buscar errores, valores nulos u outliers. Deberíamos tratar de discretizar o eliminar columnas. Incluso realizar este último paso varias veces para comprobar los diferentes resultados y elegir el que mejor rendimiento nos dé. De todos modos contiene algún valor nulo que procederemos a ignorar.
Vamos a visualizar la estructura y resumen de los datos
if (!require('palmerpenguins')) install.packages('palmerpenguins')
library(palmerpenguins)
palmerpenguins::penguins## # A tibble: 344 × 8
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Adelie Torgersen 39.1 18.7 181 3750
## 2 Adelie Torgersen 39.5 17.4 186 3800
## 3 Adelie Torgersen 40.3 18 195 3250
## 4 Adelie Torgersen NA NA NA NA
## 5 Adelie Torgersen 36.7 19.3 193 3450
## 6 Adelie Torgersen 39.3 20.6 190 3650
## 7 Adelie Torgersen 38.9 17.8 181 3625
## 8 Adelie Torgersen 39.2 19.6 195 4675
## 9 Adelie Torgersen 34.1 18.1 193 3475
## 10 Adelie Torgersen 42 20.2 190 4250
## # ℹ 334 more rows
## # ℹ 2 more variables: sex <fct>, year <int>summary(penguins)## species island bill_length_mm bill_depth_mm
## Adelie :152 Biscoe :168 Min. :32.10 Min. :13.10
## Chinstrap: 68 Dream :124 1st Qu.:39.23 1st Qu.:15.60
## Gentoo :124 Torgersen: 52 Median :44.45 Median :17.30
## Mean :43.92 Mean :17.15
## 3rd Qu.:48.50 3rd Qu.:18.70
## Max. :59.60 Max. :21.50
## NA's :2 NA's :2
## flipper_length_mm body_mass_g sex year
## Min. :172.0 Min. :2700 female:165 Min. :2007
## 1st Qu.:190.0 1st Qu.:3550 male :168 1st Qu.:2007
## Median :197.0 Median :4050 NA's : 11 Median :2008
## Mean :200.9 Mean :4202 Mean :2008
## 3rd Qu.:213.0 3rd Qu.:4750 3rd Qu.:2009
## Max. :231.0 Max. :6300 Max. :2009
## NA's :2 NA's :2Como se puede comprobar, esta base de datos está pensada para
problemas de clasificación supervisada que pretende clasificar cada tipo
de pingüino en una de las tres clases o especies existentes (Adelie,
Gentoo o Chinstrap). Como en este ejemplo vamos a usar un método no
supervisado, transformaremos el problema supervisado original en uno
no supervisado. Para conseguirlo no usaremos la columna
species, que es la variable que se quiere predecir. Por lo
tanto, intentaremos encontrar agrupaciones usando únicamente los cuatro
atributos numéricos que caracterizan a cada especie de pingüino.
x <- na.omit(penguins[,3:6]) Cargamos los datos y nos quedamos
únicamente con las cuatro columnas que definen a cada especie.
x_clean <- na.omit(penguins[,c(1,3,4,5,6)])
x_species <- x_clean[,1]
x <- x_clean[,2:5]Planteamos ahora un ejemplo más realista en el que inicialmente no conocemos el número óptimo de clústers. Empecemos probamos con varios valores.
# Silhouette method
max_k <- 10
silhouettes <- numeric(max_k)
for (k in 2:max_k) {
y_cluster <- kmeans(x, centers = k)$cluster
silhouettes[k] <- calcular_silueta(y_cluster, x)
}Mostramos en un gráfica los valores de las siluetas media de cada prueba para comprobar que número de clústers es el mejor.
plot(2:max_k, silhouettes[2:max_k], type = "b", pch = 19, frame = FALSE,
xlab = "Número de clusters (k)", ylab = "Silhouette Score",
main = "Silhouette para determinar el número óptimo de k")
Los valores de la silueta pueden fluctuar en el rango [-1,1], siendo valores cercanos a 1 indicativos de homogeneidad en los grupos y por el contrario valores de la silueta cercanos a -1 son indicativos de poca homogeneidad en los grupos, de modo que quisiéramos encontrarnos en un rango razonablemente cerca de 1.
En el caso de nuestro juego de datos, a pesar de que uno esperaría obtener un valor óptimo para k = 3, parece que del gráfico se desprende que es mejor k = 2.
Sin embargo, merece la pena observar que a partir de k = 3 la pérdida de homogeneidad es relativamente pequeña ya que se mantiene estable en el rango [0.50, 0.56]. Este hecho podría ser un argumento para seleccionar k = 3.
Otra forma de evaluar cual es el mejor número de clústers es considerar el mejor modelo, aquel que ofrece la menor suma de los cuadrados de las distancias de los puntos de cada grupo con respecto a su centro (withinss), con la mayor separación entre centros de grupos (betweenss). Como se puede comprobar es una idea conceptualmente similar a la silueta. Una manera común de hacer la selección del número de clústers consiste en aplicar el método elbow (codo), que no es más que la selección del número de clústers en base a la inspección de la gráfica que se obtiene al iterar con el mismo conjunto de datos para distintos valores del número de clústers. Se seleccionará el valor que se encuentra en el “codo” de la curva.
# Función para calcular la inercia intracluster (Within-cluster sum of squares)
inercia_intracluster <- function(vdata, max_k) {
wss <- numeric(max_k)
for (k in 1:max_k) {
wss[k] <- sum(kmeans(vdata, centers = k)$withinss)
}
return(wss)
}
# Elbow method
wss <- inercia_intracluster(x, max_k)
plot(1:max_k, wss, type = "b", pch = 19, frame = FALSE,
xlab = "Número de clusters (k)", ylab = "Within-cluster sum of squares",
main = "Elbow Method para determinar el número óptimo de k")
En este caso el número óptimo de clústers son 4 que es cuando la curva comienza a estabilizarse.
También se puede usar la función kmeansruns del paquete fpc que ejecuta el algoritmo kmeans con un conjunto de valores, para después seleccionar el valor del número de clústers que mejor funcione de acuerdo a dos criterios: la silueta media (“asw”) y Calinski-Harabasz (“ch”).
if (!require('fpc')) install.packages('fpc')
library(fpc)
fit_ch <- kmeansruns(x, krange = 1:10, criterion = "ch")
fit_asw <- kmeansruns(x, krange = 1:10, criterion = "asw") Podemos comprobar el valor con el que se ha obtenido el mejor resultado y también mostrar el resultado obtenido para todos los valores de k usando ambos criterios
fit_ch$bestk## [1] 10fit_asw$bestk## [1] 2plot(1:10,fit_ch$crit,type="o",col="blue",pch=0,xlab="Número de clústers",ylab="Criterio Calinski-Harabasz")
plot(1:10,fit_asw$crit,type="o",col="blue",pch=0,xlab="Número de clústers",ylab="Criterio silueta media")
Los resultados son muy parecidos a los que hemos obtenido anteriormente. Con el criterio de la silueta media se obtienen dos clústers y con el Calinski-Harabasz se obtienen 3.
Como se ha comprobado, conocer el número óptimo de clústers no es un problema fácil. Tampoco lo es la evaluación de los modelos de agregación.
Como en el caso que estudiamos sabemos que los datos pueden ser agrupados en 3 clases o especies, vamos a ver cómo se ha comportado kmeans en el caso de pedirle 3 clústers. Para eso comparamos visualmente los campos dos a dos, con el valor real que sabemos está almacenado en el campo “species” del dataset original.
(Aclaramos que obviamente no acostumbra a pasar que conozcamos de forma previa el número de clústers óptimo. Este ejemplo lo planteamos con finalidades didácticas y con voluntad de experimentar)
penguins3clusters <- kmeans(x, 3)
# bill_length y bill_depth
plot(x[c(1,2)],
col = penguins3clusters$cluster, # Color by k-means cluster
main = "Clasificación k-means", # Plot title
xlab = "bill_length",
ylab = "bill_depth")
# Add a legend for the clusters
legend("topright",
legend = unique(penguins3clusters$cluster), # Cluster labels
col = unique(penguins3clusters$cluster), # Colors corresponding to clusters
pch = 19, # Point symbol in the legend
title = "K-means Clusters")
# Plot with real species classification
plot(x[c(1,2)],
col = as.factor(x_species$species), # Color by real species
main = "Clasificación real", # Plot title
xlab = "bill_length",
ylab = "bill_depth")
# Add a legend for species
legend("topright",
legend = unique(x_species$species), # Species labels
col = unique(as.factor(x_species$species)), # Colors corresponding to species
pch = 19, # Point symbol in the legend
title = "Species")
Podemos observar en el siguiente gráfico, que flipper_length y body_mass no son buenos indicadores para diferenciar a las tres subespecies, dado que dos de las subespecies están demasiado mezcladas para poder diferenciar nada.
# flipper_length y body_mass
plot(x[c(3,4)],
col = penguins3clusters$cluster, # Color by k-means cluster
main = "Clasificación k-means", # Plot title
xlab = "flipper_length",
ylab = "body_mass")
# Add a legend for the clusters
legend("topright",
legend = unique(penguins3clusters$cluster), # Cluster labels
col = unique(penguins3clusters$cluster), # Colors corresponding to clusters
pch = 19, # Point symbol in the legend
title = "K-means Clusters")
# Plot with real species classification
plot(x[c(3,4)],
col = as.factor(x_species$species), # Color by real species
main = "Clasificación real", # Plot title
xlab = "flipper_length",
ylab = "body_mass")
# Add a legend for species
legend("topright",
legend = unique(x_species$species), # Species labels
col = unique(as.factor(x_species$species)), # Colors corresponding to species
pch = 19, # Point symbol in the legend
title = "Species")
# bill_length y flipper_length
plot(x[c(1,3)],
col = penguins3clusters$cluster, # Color by k-means cluster
main = "Clasificación k-means", # Plot title
xlab = "bill_length",
ylab = "flipper_length")
# Add a legend for the clusters
legend("topright",
legend = unique(penguins3clusters$cluster), # Cluster labels
col = unique(penguins3clusters$cluster), # Colors corresponding to clusters
pch = 19, # Point symbol in the legend
title = "K-means Clusters")
# Plot with real species classification
plot(x[c(1,3)],
col = as.factor(x_species$species), # Color by real species
main = "Clasificación real", # Plot title
xlab = "bill_length",
ylab = "flipper_length")
# Add a legend for species
legend("topright",
legend = unique(x_species$species), # Species labels
col = unique(as.factor(x_species$species)), # Colors corresponding to species
pch = 19, # Point symbol in the legend
title = "Species")
Las dos medidas de bill parecen lograr mejores resultados al dividir las tres especies de pingüinos. El grupo formado por los puntos negros que ha encontrado el algoritmo coincide con los de la especie Adelie. Los otros dos grupos sin embargo se entremezclan algo más, y hay ciertos puntos que se clasifican como Gentoo (verde) cuando en realidad son Chinstrap (rojo).
Una buena técnica que ayuda a entender los grupos que se han formado, es mirar de darles un nombre. Cómo por ejemplo:
Esto nos ayuda a entender cómo están formados los grupos y a referirnos a ellos en análisis posteriores.
Todo esto nos indica que el número de grupos o clúsers en un juego de datos no es un aspecto que podamos asegurar que siempre vamos a encontrar de forma precisa y objetiva, bien al contrario es un ámbito que requiere de análisis en sí mismo.
Os compartimos en el siguiente enlace un material didáctico
complementario que os puede ayudar a profundizar en el tema de la
selección del número de clústers más adecuado para un juego de
datos:
datascience.recursos.uoc.edu
Como continuación del estudio podríamos seguir experimentando combinando en gráficos similares a los anteriores. En definitiva se trataría en este punto de profundizar más en el conocimiento de las propiedades de las diferentes características o columnas del juego de datos.
En este ejemplo vamos ha trabajar los algoritmos DBSCAN y OPTICS como métodos de clustering que permiten la generación de grupos no radiales a diferencia de k-means. Veremos que su parámetro de entrada más relevante es minPts que define la mínima densidad aceptada alrededor de un centroide.
Incrementar este parámetro nos permitirá reducir el ruido (observaciones no asignadas a ningún cluster), en cualquier caso empezaremos por construir nuestro propio juego de datos en el que dibujaremos 4 zonas de puntos diferenciadas.
if (!require('dbscan')) install.packages('dbscan')
library(dbscan)
set.seed(2)
n <- 400
x <- cbind(
x = runif(4, 0, 1) + rnorm(n, sd=0.1),
y = runif(4, 0, 1) + rnorm(n, sd=0.1)
)
plot(x, col=rep(1:4, time = 100))
Una de las primeras actividades que realiza el algoritmo es ordenar las observaciones de forma que los puntos más cercanos se conviertan en vecinos en el ordenamiento. Se podría pensar como una representación numérica del dendograma de una agrupación jerárquica.
### Lanzamos el algoritmo OPTICS dejando el parámetro eps con su valor por defecto y fijando el criterio de vecindad en 10
res <- optics(x, minPts = 10)
res## OPTICS ordering/clustering for 400 objects.
## Parameters: minPts = 10, eps = 0.193786846197958, eps_cl = NA, xi = NA
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
## eps_cl, xi### Obtenemos la ordenación de las observaciones o puntos
res$order## [1] 1 363 209 349 337 301 357 333 321 285 281 253 241 177 153 57 257 29
## [19] 77 169 105 293 229 145 181 385 393 377 317 381 185 117 101 9 73 237
## [37] 397 369 365 273 305 245 249 309 157 345 213 205 97 49 33 41 193 149
## [55] 17 83 389 25 121 329 5 161 341 217 189 141 85 53 225 313 289 261
## [73] 221 173 69 61 297 125 81 133 129 197 109 137 59 93 165 89 21 13
## [91] 277 191 203 379 399 375 351 311 235 231 227 71 11 299 271 291 147 55
## [109] 23 323 219 275 47 263 3 367 331 175 87 339 319 251 247 171 111 223
## [127] 51 63 343 303 207 151 391 359 287 283 215 143 131 115 99 31 183 43
## [145] 243 199 79 27 295 67 347 255 239 195 187 139 107 39 119 179 395 371
## [163] 201 123 159 91 211 355 103 327 95 7 167 35 267 155 387 383 335 315
## [181] 259 135 15 113 279 373 4 353 265 127 45 37 19 276 224 361 260 288
## [199] 336 368 348 292 268 252 120 108 96 88 32 16 340 156 388 372 356 332
## [217] 304 220 188 168 136 124 56 236 28 244 392 184 76 380 232 100 116 112
## [235] 256 72 8 280 64 52 208 172 152 148 360 352 192 160 144 284 216 48
## [253] 84 92 36 20 212 272 264 200 128 80 180 364 196 12 132 40 324 308
## [271] 176 164 68 316 312 384 300 344 328 248 204 140 296 24 320 228 60 44
## [289] 233 65 400 376 240 163 104 396 307 75 14 325 269 262 234 382 294 206
## [307] 198 374 310 362 318 386 358 330 278 210 298 282 122 98 34 26 174 142
## [325] 46 6 62 118 190 202 114 322 286 38 242 394 342 266 162 130 30 182
## [343] 2 74 314 290 246 194 170 126 158 378 350 254 226 214 70 18 10 366
## [361] 354 186 150 86 306 102 338 346 134 250 138 94 78 390 274 58 42 258
## [379] 66 90 146 370 222 218 326 82 110 270 334 178 166 398 22 50 238 106
## [397] 154 302 230 54Otro paso muy interesante del algoritmo es la generación de un diagrama de alcanzabilidad o reachability plot, en el que se aprecia de una forma visual la distancia de alcanzabilidad de cada punto.
Conceptos clave para interpretar el gráfico
Como interpretar un Reachability Plot
- Eje X (orden de procesamiento): Refleja el orden que
los puntos han sido visitados según el algoritmo OPTICS.
- Eje Y (reachability distance): Indica como de
accesible es cada punto desde el vecino denos más cercano.
Clústers densos = valles
Cuando hay una zona donde los valores de reachability son bajos
durante un intervalo largo, indica un clúster denso, los puntos son
cercanos entre ellos y fácilmente accesibles. Cuanto más ancho y
profundo sea el valle, más grande y denso es el clúster.
Separación entre clústers = picos
- Cuando el valor de reachability sube repentinamente (pico),
indica que el algoritmo está entrando o saliendo de un clúster. Este
pico es como un límite natural entre grupos.
- Si el reachability distance sube mucho, puede indicar ruido o
puntos extremos.
Ruido o outliers
Puntos con reachability distance muy altas (al final de un
pico, o sin valle claro) pueden ser considerados ruido, no forman parte
de ningún clúster denso.
### Gráfica de alcanzabilidad
plot(res)
En este gràfico se aprecian claramente cuatro valles, que nos indican las cuatro zonas con más densidad de puntos (clusters).
Veamos otra representación del diagrama de alcanzabilidad, donde podemos observar las trazas de las distancias entre puntos cercanos del mismo cluster y entre clusters distintos.
### Dibujo de las trazas que relacionan puntos
plot(x, col = "grey")
polygon(x[res$order,])
Otro ejercicio interesante a realizar es extraer una agrupación de la ordenación realizada por OPTICS similar a lo que DBSCAN hubiera generado estableciendo el parámetro eps en eps_cl = 0.065. En este sentido animamos al estudiante a experimentar con diferentes valores de este parámetro.
### Extracción de un clustering DBSCAN cortando la alcanzabilidad en el valor eps_cl
res <- extractDBSCAN(res, eps_cl = .065)
res## OPTICS ordering/clustering for 400 objects.
## Parameters: minPts = 10, eps = 0.193786846197958, eps_cl = 0.065, xi = NA
## The clustering contains 4 cluster(s) and 92 noise points.
##
## 0 1 2 3 4
## 92 81 84 72 71
##
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
## eps_cl, xi, clusterplot(res) ## negro indica ruido
Observamos en el gráfico anterior como se han coloreado los 4 clusters y en negro se mantienen los valores outliers o extremos.
Seguimos adelante con una representación gráfica que nos muestra los clusters mediante formas convexas.
hullplot(x, res)
Repetimos el experimento anterior incrementando el parámetro epc_c, veamos como el efecto que produce es la concentración de clusters ya que flexibilizamos la condición de densidad.
### Incrementamos el parámetro eps
res <- extractDBSCAN(res, eps_cl = .1)
res## OPTICS ordering/clustering for 400 objects.
## Parameters: minPts = 10, eps = 0.193786846197958, eps_cl = 0.1, xi = NA
## The clustering contains 2 cluster(s) and 9 noise points.
##
## 0 1 2
## 9 295 96
##
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
## eps_cl, xi, clusterplot(res)
hullplot(x, res)
Veamos ahora una variante de la extracción DBSCN anterior. En ella el parámetro xi nos va a servir para clasificar los clusters en función del cambio en la densidad relativa de los mismos.
### Extracción del clustering jerárquico en función de la variación de la densidad por el método xi
res <- extractXi(res, xi = 0.05)
res## OPTICS ordering/clustering for 400 objects.
## Parameters: minPts = 10, eps = 0.193786846197958, eps_cl = NA, xi = 0.05
## The clustering contains 7 cluster(s) and 1 noise points.
##
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
## eps_cl, xi, cluster, clusters_xiplot(res)
hullplot(x, res)
#Elimibna mos las variables de entorno
rm(list = ls())Los ejercicios se realizarán en base al juego de datos Hawks presente en el paquete R Stat2Data.
Los estudiantes y el profesorado del Cornell College en Mount Vernon, Iowa, recogieron datos durante muchos años en el mirador de halcones del lago MacBride, cerca de Iowa City, en el estado de Iowa. El conjunto de datos que analizamos aquí es un subconjunto del conjunto de datos original, utilizando sólo aquellas especies para las que había más de 10 observaciones. Los datos se recogieron en muestras aleatorias de tres especies diferentes de halcones: Colirrojo, Gavilán y Halcón de Cooper.
Hemos seleccionado este juego de datos por su parecido con el juego de datos penguins y por su potencial a la hora de aplicarle algoritmos de minería de datos no supervisados.
if (!require('Stat2Data')) install.packages('Stat2Data')
library(Stat2Data)
data("Hawks")
summary(Hawks)## Month Day Year CaptureTime ReleaseTime
## Min. : 8.000 Min. : 1.00 Min. :1992 11:35 : 14 :842
## 1st Qu.: 9.000 1st Qu.: 9.00 1st Qu.:1995 13:30 : 14 11:00 : 2
## Median :10.000 Median :16.00 Median :1999 11:45 : 13 11:35 : 2
## Mean : 9.843 Mean :15.74 Mean :1998 12:10 : 13 12:05 : 2
## 3rd Qu.:10.000 3rd Qu.:23.00 3rd Qu.:2001 14:00 : 13 12:50 : 2
## Max. :11.000 Max. :31.00 Max. :2003 13:05 : 12 13:32 : 2
## (Other):829 (Other): 56
## BandNumber Species Age Sex Wing Weight
## : 2 CH: 70 A:224 :576 Min. : 37.2 Min. : 56.0
## 1142-09240: 1 RT:577 I:684 F:174 1st Qu.:202.0 1st Qu.: 185.0
## 1142-09241: 1 SS:261 M:158 Median :370.0 Median : 970.0
## 1142-09242: 1 Mean :315.6 Mean : 772.1
## 1142-18229: 1 3rd Qu.:390.0 3rd Qu.:1120.0
## 1142-19209: 1 Max. :480.0 Max. :2030.0
## (Other) :901 NA's :1 NA's :10
## Culmen Hallux Tail StandardTail
## Min. : 8.6 Min. : 9.50 Min. :119.0 Min. :115.0
## 1st Qu.:12.8 1st Qu.: 15.10 1st Qu.:160.0 1st Qu.:162.0
## Median :25.5 Median : 29.40 Median :214.0 Median :215.0
## Mean :21.8 Mean : 26.41 Mean :198.8 Mean :199.2
## 3rd Qu.:27.3 3rd Qu.: 31.40 3rd Qu.:225.0 3rd Qu.:226.0
## Max. :39.2 Max. :341.40 Max. :288.0 Max. :335.0
## NA's :7 NA's :6 NA's :337
## Tarsus WingPitFat KeelFat Crop
## Min. :24.70 Min. :0.0000 Min. :0.000 Min. :0.0000
## 1st Qu.:55.60 1st Qu.:0.0000 1st Qu.:2.000 1st Qu.:0.0000
## Median :79.30 Median :1.0000 Median :2.000 Median :0.0000
## Mean :71.95 Mean :0.7922 Mean :2.184 Mean :0.2345
## 3rd Qu.:87.00 3rd Qu.:1.0000 3rd Qu.:3.000 3rd Qu.:0.2500
## Max. :94.00 Max. :3.0000 Max. :4.000 Max. :5.0000
## NA's :833 NA's :831 NA's :341 NA's :343if (!require('DT')) install.packages('DT')
library(DT)
datatable(head(Hawks, 10),
caption = "Primeras 10 filas del dataset Hawks",
options = list(pageLength = 5))Tenemos la opción de descargar el archivo Hawks en local:
# Descargamos y guardamos una copia en local
download.file("https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/Stat2Data/Hawks.csv",
destfile = "Hawks.csv")
# Ahora puedes cargarlo desde el archivo guardado
hawks <- read.csv("Hawks.csv")Presenta el juego de datos, nombre y significado de cada columna, así como las distribuciones de sus valores.
El conjunto de datos Hawks se recopiló con el propósito de estudiar las características físicas y biológicas de diferentes especies de halcones en Iowa. Analizar estas variables puede ayudar a los investigadores a entender mejor las diferencias entre especies, las variaciones según la edad y el sexo, y otros factores relevantes en la biología de estas aves rapaces.
Fuente:
# Fuente:
browseURL("https://search.r-project.org/CRAN/refmans/Stat2Data/html/Hawks.html?")El conjunto de datos Hawks forma parte del paquete de R Stat2Data y contiene información recopilada sobre tres especies de halcones: el halcón de Cooper (CH), el halcón de cola roja (RT) y el gavilán (SS).
Descripción de las columnas y sus significados:
Month: Mes de captura, codificado del 8 (agosto) al 12 (diciembre).
Day: Día del mes en que se realizó la captura.
Year: Año de la captura, entre 1992 y 2003.
CaptureTime: Hora de captura en formato HH:MM.
ReleaseTime: Hora de liberación en formato HH:MM.
BandNumber: Código de identificación del anillo colocado en el ave.
Species: Especie del halcón; CH para halcón de Cooper, RT para halcón de cola roja y SS para gavilán.
Age: Edad del halcón; A para adulto y I para inmaduro.
Sex: Sexo del halcón; F para hembra y M para macho.
Wing: Longitud del ala primaria desde la punta hasta la muñeca donde se une, medida en milímetros.
Weight: Peso corporal en gramos.
Culmen: Longitud de la parte superior del pico desde la punta hasta donde se une con la parte carnosa del ave, medida en milímetros.
Hallux^: Longitud del espolón (garra trasera) utilizada para cazar, medida en milímetros.
Tail: Medida relacionada con la longitud de la cola, en milímetros.
StandardTail: Medida estándar de la longitud de la cola, en milímetros.
Tarsus: Longitud del tarso (hueso básico del pie), en milímetros.
WingPitFat: Cantidad de grasa en la axila del ala.
KeelFat: Cantidad de grasa en la quilla (esternón), evaluada al tacto.
Crop: Cantidad de material en el buche, codificada de 1 (lleno) a 0 (vacío).
# Tamaño del dataset Hawks
dim(Hawks)## [1] 908 19colnames(Hawks)## [1] "Month" "Day" "Year" "CaptureTime" "ReleaseTime"
## [6] "BandNumber" "Species" "Age" "Sex" "Wing"
## [11] "Weight" "Culmen" "Hallux" "Tail" "StandardTail"
## [16] "Tarsus" "WingPitFat" "KeelFat" "Crop"El dataset Hawks contiedne 908 observaciones / registros y 19 variables / columnas.
# Observamos la estructura de los datos una vez limpiados
structure_Hawks = str(Hawks)## 'data.frame': 908 obs. of 19 variables:
## $ Month : int 9 9 9 9 9 9 9 9 9 9 ...
## $ Day : int 19 22 23 23 27 28 28 29 29 30 ...
## $ Year : int 1992 1992 1992 1992 1992 1992 1992 1992 1992 1992 ...
## $ CaptureTime : Factor w/ 308 levels " ","1:15","1:31",..: 181 25 138 42 62 71 181 88 261 192 ...
## $ ReleaseTime : Factor w/ 60 levels ""," ","10:20",..: 1 2 2 2 2 2 2 2 2 2 ...
## $ BandNumber : Factor w/ 907 levels " ","1142-09240",..: 856 857 858 809 437 280 859 860 861 281 ...
## $ Species : Factor w/ 3 levels "CH","RT","SS": 2 2 2 1 3 2 2 2 2 2 ...
## $ Age : Factor w/ 2 levels "A","I": 2 2 2 2 2 2 2 1 1 2 ...
## $ Sex : Factor w/ 3 levels "","F","M": 1 1 1 2 2 1 1 1 1 1 ...
## $ Wing : num 385 376 381 265 205 412 370 375 412 405 ...
## $ Weight : int 920 930 990 470 170 1090 960 855 1210 1120 ...
## $ Culmen : num 25.7 NA 26.7 18.7 12.5 28.5 25.3 27.2 29.3 26 ...
## $ Hallux : num 30.1 NA 31.3 23.5 14.3 32.2 30.1 30 31.3 30.2 ...
## $ Tail : int 219 221 235 220 157 230 212 243 210 238 ...
## $ StandardTail: int NA NA NA NA NA NA NA NA NA NA ...
## $ Tarsus : num NA NA NA NA NA NA NA NA NA NA ...
## $ WingPitFat : int NA NA NA NA NA NA NA NA NA NA ...
## $ KeelFat : num NA NA NA NA NA NA NA NA NA NA ...
## $ Crop : num NA NA NA NA NA NA NA NA NA NA ...Número total de valores NA
# Verificamos la existencia de valors NA por columnas
colSums(is.na(Hawks))## Month Day Year CaptureTime ReleaseTime BandNumber
## 0 0 0 0 0 0
## Species Age Sex Wing Weight Culmen
## 0 0 0 1 10 7
## Hallux Tail StandardTail Tarsus WingPitFat KeelFat
## 6 0 337 833 831 341
## Crop
## 343Número total de valores vacíos
# Verificamos la existencia de valors nulos o vacíos por columnas
colSums(Hawks=="")## Month Day Year CaptureTime ReleaseTime BandNumber
## 0 0 0 0 1 0
## Species Age Sex Wing Weight Culmen
## 0 0 576 NA NA NA
## Hallux Tail StandardTail Tarsus WingPitFat KeelFat
## NA 0 NA NA NA NA
## Crop
## NAObservamos que existen más del 35% de los valores de las columnas StandardTail, KeelFat y Crop faltantes o nulos, además de más del 90% en las columnas Tarsus y WingPitFat. Por este motivo, nos vemos obligados a descartar del estudio las columnas mencionadas para no crear sesgos en los datos. Trataremos los valores faltantes de las demás columnas aplicando la media o el mode de la variable tratada, según el tipo de dato.
# Detección de registros duplicados
duplicados <- duplicated(Hawks)
sum(duplicados)## [1] 0Además creamos un nuevo dataset con las columnas que vamos a trabajar sin las variables con valores nulos o vacios.
# Creamos un nuevo dataset con solo las columnas clave a estudiar
Hawks2 <- Hawks[, c("Month", "Day", "Year", "CaptureTime",
"BandNumber", "Species", "Age",
"Wing", "Weight", "Culmen", "Hallux", "Tail")]# Reemplazamos los valores NA o vacios por la media o mode de la columna que corresponda
cols_na <- c("Wing", "Weight", "Culmen", "Hallux")
for (col in cols_na) {
media_valor <- mean(Hawks[[col]], na.rm = TRUE)
Hawks2[[col]][is.na(Hawks2[[col]])] <- media_valor
}# Verificamos la existencia de valors NA por columnas
colSums(is.na(Hawks2))## Month Day Year CaptureTime BandNumber Species
## 0 0 0 0 0 0
## Age Wing Weight Culmen Hallux Tail
## 0 0 0 0 0 0# Observamos la estructura de los datos una vez limpiados
structure_Hawks2 = str(Hawks2)## 'data.frame': 908 obs. of 12 variables:
## $ Month : int 9 9 9 9 9 9 9 9 9 9 ...
## $ Day : int 19 22 23 23 27 28 28 29 29 30 ...
## $ Year : int 1992 1992 1992 1992 1992 1992 1992 1992 1992 1992 ...
## $ CaptureTime: Factor w/ 308 levels " ","1:15","1:31",..: 181 25 138 42 62 71 181 88 261 192 ...
## $ BandNumber : Factor w/ 907 levels " ","1142-09240",..: 856 857 858 809 437 280 859 860 861 281 ...
## $ Species : Factor w/ 3 levels "CH","RT","SS": 2 2 2 1 3 2 2 2 2 2 ...
## $ Age : Factor w/ 2 levels "A","I": 2 2 2 2 2 2 2 1 1 2 ...
## $ Wing : num 385 376 381 265 205 412 370 375 412 405 ...
## $ Weight : num 920 930 990 470 170 1090 960 855 1210 1120 ...
## $ Culmen : num 25.7 21.8 26.7 18.7 12.5 ...
## $ Hallux : num 30.1 26.4 31.3 23.5 14.3 ...
## $ Tail : int 219 221 235 220 157 230 212 243 210 238 ...Graficamos las variables numéricas para observar sus distribuciones en base a la Species y a la varia ble Age (variable categórica) y
# Distribución de variables numéricas individualmente
library(ggplot2)
vars_numericas <- c("Wing", "Weight", "Culmen", "Hallux", "Tail")
colors <- c("lightgreen", "tomato", "skyblue", "purple", "pink")
for (i in seq_along(vars_numericas)) {
var <- vars_numericas[i]
print(
ggplot(Hawks2, aes_string(x = var)) +
geom_histogram(bins = 30, fill = colors[i], color = "black") +
ggtitle(paste("Distribución de", var)) +
theme_minimal()
)
}
# Distribución de variables numéricas s/ Species
vars_numericas <- c("Wing", "Weight", "Culmen", "Hallux", "Tail")
colors <- c("lightgreen", "tomato", "skyblue", "purple", "pink")
# 1. Distribución de numéricas por ESPECIE
for (i in seq_along(vars_numericas)) {
var <- vars_numericas[i]
print(
ggplot(Hawks2, aes_string(x = "Species", y = var, fill = "Species")) +
geom_boxplot() +
scale_fill_manual(values = colors[1:3]) +
ggtitle(paste("Distribución de", var, "por especie")) +
theme_minimal()
)
}
# Distribución de variables numéricas s/ Age
for (i in seq_along(vars_numericas)) {
var <- vars_numericas[i]
print(
ggplot(Hawks2, aes_string(x = "Age", y = var, fill = "Age")) +
geom_boxplot() +
scale_fill_manual(values = colors[4:5]) +
ggtitle(paste("Distribución de", var, "por edad (A = Adulto, I = Inmaduro)")) +
theme_minimal()
)
}
Wing La distribución de la longitud del ala (Wing) no es unimodal, sino bimodal, con dos picos bien definidos: - Un primer grupo concentrado entre 180 y 240 mm. - Un segundo grupo más denso entre 360 y 420 mm. Esto sugiere la presencia de dos poblaciones diferenciadas, probablemente especies con distinta morfología, o diferenciada por la variable sexo (la cual se ha dejado fuera del alcance del estudio por la falta de consistencia de sus datos). La bimodalidad es un fuerte indicador de heterogeneidad estructural en la muestra.
Se confirma que según el histograma, RT (Red-tailed Hawk) muestra las alas más largas, con una mediana cerca de 390 mm y valores entre 350–450 mm. CH (Cooper’s Hawk) tiene alas más cortas, con mediana cercana a 260 mm. Y, SS (Sharp-shinned Hawk) tiene las alas más pequeñas, centradas en torno a 200 mm. Esta segmentación es consistente: RT es la más grande y SS la más pequeña. Las diferencias son claras y significativas, con poco solapamiento.
El análisis por edad muestra una diferencia menos marcada que por especies. Aunque, los adultos (A) tienen alas más largas que la mediana comparada con los inmaduros (I), existe un solapamiento considerable. Ambas categorías presentan rangos amplios, lo cual sugiere que la variable Wing está más influida por la especie que por la edad. Aunque hay desarrollo de las alas con la edad, la diferencia no es tan pronunciada como entre especies.
Weight
La distribución del peso muestra dos grupos diferenciados al igual que en Wing: El primer pico se observa entre 100–300 g, que indica una agrupación de halcones más pequeños. El segundo grupo, más amplio, se concentra entre 900 y 1300 g, el cual corresponde a individuos significativamente más grandes. Esta bimodalidad indica nuevamente la presencia de al menos dos especies con diferente peso. También se aprecian valores extremos por encima de los 2000 g y por debajo de los 100 g, los cuales podrían representar outliers o errores de medición.
Según la especie, el histograma muestra que RT (Red-tailed Hawk) presenta los valores más altos, con la mediana cerca de 1100 g y una distribución entre 800–1300 g. Algunos valores superan los 2000 g. La especie CH (Cooper’s Hawk) tiene valores intermedios, con una mediana alrededor de 400 g. Y, por último, SS (Sharp-shinned Hawk) es la más ligera, con pesos entre 100 y 300 g, y una mediana inferior a 250 g. Esta clara separación entre especies evidencia que el peso es un excelente discriminador entre halcones.
La distribución por edad muestra que los Los adultos (A) y los inmaduros (I) presentan una mediana similar. Ambas categorías tienen una distribución extensa de valores, lo que indica que, a diferencia del caso de Wing, el peso no varía significativamente entre edades dentro de la misma especie. La variación parece deberse más a la especie que a la madurez.
Culmen
La distribución de Culmen presenta una clara bimodalidad, como en los casos anteriores. Un primer grupo con valores entre 8 y 15 mm, lo que sugiere una especie con pico corto (posiblemente SS) y un segundo grupo con alta densidad entre 25 y 30 mm, representando aves con picos más largos. Esta bimodalidad refleja probablemente diferencias por especie. Existen algunos valores dispersos entre 15 y 25 mm y algunos outliers por encima de 35 mm.
RT (Red-tailed Hawk) tiene los valores más altos de Culmen, con una mediana cercana a 27 mm. Muestra una distribución casi homogénea y valores atípicos por encima de 35 mm. CH (Cooper’s Hawk) presenta valores intermedios, con una mediana alrededor de 17 mm. Y, SS (Sharp-shinned Hawk) con los picos más cortos, centrados en torno a 12 mm. Se confirma que Culmen es una variable diferenciadora por especies. Existe muy poco solapamiento entre categorías.
Tanto los adultos (A) como los inmaduros (I) presentan distribuciones muy similares.Las medianas no difieren excesivamente, con una mínima tendencia a valores mayores para los adultos. Gran dispersión en ambos grupos.La longitud del pico es sensible a la especie y poco dependiente de la edad.
Hallux
La distribución está altamente sesgada hacia la izquierda. La mayoría de las observaciones están concentradas entre 10 y 40 mm. Se observan outliers extremos, incluso mayores a 300 mm, lo cual podría deberse a errores de medición o registros atípicos.
RT (Red-tailed Hawk) presenta los valores más elevados de Hallux, con una mediana alrededor de 35 mm. Se observan outliers extemos por encima de 300 mm. CH (Cooper’s Hawk) tiene valores más moderados, con una mediana en torno a 30 mm. Y, por último observamos que SS (Sharp-shinned Hawk) tiene el Hallux más corto, con valores por debajo de 20 mm en la mayoría de los casos. Se confirma que la variable Hallux es útil para la diferenciación entre especies, aunque la presencia de outliers puede distorsionar los resultados.
No se observan diferencias significativas entre adultos (A) e inmaduros (I). Las medianas son muy similares, y ambos grupos comparten una amplia dispersión. La presencia de valores atípicos también se distribuye en ambos grupos.
Tail La variable Tail presenta una distribución bimodal, con grupos en 130-170 mm y 210-240 mm. Por especie, RT tiene colas más largas, CH y SS presentan valores medios y bajos, respectivamente. Y, según edad, los Inmaduros tienden a tener una cola ligeramente más larga en la mediana, aunque la diferencia no es significativa respecto a los adultos.
# Estrudio de correslación de las variables numéricas
if (!require('corrplot')) install.packages('corrplot')
library(corrplot)
# Si tienes solo los nombres de columnas numéricas
vars_numericas <- c("Wing", "Weight", "Culmen", "Hallux", "Tail")
hawks_cor <- Hawks2[, vars_numericas]
corr_matrix <- cor(hawks_cor, use = "complete.obs")
# Visualizamos la matriz de correlación
corrplot(corr_matrix,
method = "color",
type = "upper",
tl.col = "black",
tl.srt = 30,
order = "AOE",
number.cex = 0.75,
sig.level = 0.01,
addCoef.col = "white")
###Justificación de los resultados
La matriz de correlación permite entender las relaciones lineales entre las variables numéricas del conjunto de datos de Hawks en base en el coeficiente de correlación de Pearson, valores cercanos a 1 indican alta correlación positiva, valores cercanos a 0 indican ausencia de relación lineal.
Culmen, Weight, Wing, Tail presentan entre ellas correlaciones muy altas ≥ 0.90. Lo que indica que están relacionadas con el tamaño de los halcones. Un halcón con alas más largas también suele tener mayor peso, longitud del pico y de la cola.
library(dplyr)
library(tidyr)
# Analisis multidimensional entre especies y año que fueron capturados
species_year <- Hawks2 %>%
group_by(Species, Year) %>%
summarise(total_observations = n(), .groups = "drop")
# Convertimos los años en columnas (análisis multidimensional)
species_wide <- species_year %>%
pivot_wider(names_from = Year, values_from = total_observations,
values_fill = list(total_observations = 0))
# Tabla a mostrar
print(species_wide)## # A tibble: 3 × 13
## Species `1992` `1993` `1994` `1995` `1996` `1997` `1998` `1999` `2000` `2001`
## <fct> <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
## 1 CH 1 4 4 5 3 8 2 6 8 7
## 2 RT 35 29 90 54 14 43 26 56 80 36
## 3 SS 3 6 35 15 1 29 6 24 29 41
## # ℹ 2 more variables: `2002` <int>, `2003` <int>El análisis de la distribución de especies a lo largo del tiempo nos muestra que el Red-tailed Hawk (RT) es la especie más observada en 1994 (90 registros) y en los años siguientes. Le sigue el Sharp-shinned Hawk (SS), con una tendencia más variable sobretodo en 1994, 2001 y 2003. El Cooper’s Hawk (CH) es la menos frecuente, aunque muestra un crecimiento leve desde 1992 hasta 2003. El año 1994 destaca como el de mayor actividad en general, mientras que 1996 presenta el número más bajo de registros. Estos resultados sugieren tanto patrones migratorios como posibles diferencias en el esfuerzo de muestreo según el año y la especie.
Realiza un estudio aplicando el método K-means, similar al de los ejemplos 1 y 2
Usaremos k = 3, ya que hay tres especies en el dataset (RT, CH, SS).
# Método k-means
hawks_numeric <- Hawks2[, vars_numericas]
k3_result <- kmeans(hawks_numeric, centers = 3)
Hawks2$cluster <- as.factor(k3_result$cluster)La media de cada variable numérica por clúster, lo que te permitirá caracterizar los grupos.
# Estadísticos descriptivos por clúster
aggregate(hawks_numeric, by = list(Cluster = Hawks2$cluster), FUN = mean)## Cluster Wing Weight Culmen Hallux Tail
## 1 1 373.2752 989.1030 26.04471 29.88227 217.5123
## 2 2 198.8460 196.6116 12.84358 16.83398 158.3945
## 3 3 395.2533 1279.2336 28.21262 35.09141 228.5841# Tabla cruzada entre clusters y especies
table(Hawks2$cluster, Hawks2$Species)##
## CH RT SS
## 1 4 357 6
## 2 66 6 255
## 3 0 214 0ggplot(Hawks2, aes(x = Wing, y = Weight, color = cluster)) +
geom_point(alpha = 0.6) +
labs(title = "Clustering K-means (k=3) en Hawks2", x = "Wing (mm)", y = "Weight (g)")
ggplot(Hawks2, aes(x = Wing, y = Weight, color = Species)) +
geom_point(alpha = 0.6) +
labs(title = "Distribución por especie", x = "Wing (mm)", y = "Weight (g)")
K-means con k = 3, teniendo en cuenta las diferentes especies de halcones: - CH = Cooper’s Hawk - RT = Red-Tailed Hawk - SS = Sharp-Shinned Hawk Asignamos cada halcón al clúster más cercano en función de las distancias a los centroides. E iteramos hasta estabilizar las asignaciones de los clústeres.
En el primer gráfico se muestra la separación por cluster wingd vs weight. En elk segundo gráfico, se muestra la distribución por especies con el objetivo de contrastar. Y, observamos que el Cluster 1 incluye halcones ligeros y pequeños, correspondiendo en su mayoría a la especie SS. El Clúster 2 agrega halcones grandes y pesados, dominado por la especie RT. Y, por último el Clúster 3 contiene valores intermedios, mezclando algunos CH y RT. Podemos señalar que aunque los clusteres no son puros tiene una alta correspondencia entre las variables Wing y Weight contra las diferentes especies.
# K-means con k = 2, 3, 5
hawks_numeric <- Hawks2[, vars_numericas]
# Función para aplicar K-means
kmeans_ajust <- function(k, vdata) {
fit <- kmeans(vdata, centers = k)
return(fit$cluster)
}
# Aplicamos K-means con distintos valores de k
clusters_2 <- kmeans_ajust(2, hawks_numeric)
clusters_3 <- kmeans_ajust(3, hawks_numeric)
clusters_5 <- kmeans_ajust(5, hawks_numeric)# Graficamos los clusteres
clusplot(hawks_numeric, clusters_2, color = TRUE, shade = TRUE, labels = 2, lines = 0, main = "K-means con k = 2")
clusplot(hawks_numeric, clusters_3, color = TRUE, shade = TRUE, labels = 2, lines = 0, main = "K-means con k = 3")
clusplot(hawks_numeric, clusters_5, color = TRUE, shade = TRUE, labels = 2, lines = 0, main = "K-means con k = 5")
Aplicamos K-means al dataset Hawks2 con k = 2, 3 y 5 para identificar agrupaciones. El mejor resultado se obtuvo con k = 3, ya que coincide con las tres especies reales (CH, RT, SS) y refleja de forma coherente su distinció. Con k = 2, los grupos son muy generales, y con k = 5, el modelo sobreajusta creando divisiones poco consistentes. Así, que k = 3 nos ofrece la segmentación más adecuada en este caso.
# Método silhouette
# Calculamos la matriz de distancias
dist_hawks <- daisy(hawks_numeric)
# Función para calcular la silueta media
calcular_silueta <- function(clusters, dist_matrix) {
sil <- silhouette(clusters, dist_matrix)
return(mean(sil[, 3]))
}
# Siluetas para k = 2, 3 y 5
sil_2 <- calcular_silueta(clusters_2, dist_hawks)
sil_3 <- calcular_silueta(clusters_3, dist_hawks)
sil_5 <- calcular_silueta(clusters_5, dist_hawks)
# Mostrar resultados
cat("Silueta media para k = 2:", round(sil_2, 3), "\n")## Silueta media para k = 2: 0.797cat("Silueta media para k = 3:", round(sil_3, 3), "\n")## Silueta media para k = 3: 0.65cat("Silueta media para k = 5:", round(sil_5, 3), "\n")## Silueta media para k = 5: 0.573K = 2 (0.797): Es el valor más alto, indica buena separación y cohesión de los grupos. Los clústers son compactos y bien diferenciados. K = 3 (0.65): Aunque aceptable, es menos óptimo. Existe solapamiento entre grupos o estructuras menos claras. K = 5 (0.573): Es el más bajo, lo que sugiere una partición menos definida, posible sobreajuste. El número de clústers más adecuado según la silueta media para este conjunto de datos numéricos es k=2, ya que proporciona la estructura más clara sin solapamientos. La silueta media más cercana a 1 indica que el número de clústers elegido reflejará una buena cohesión y externamente bien diferenciados.
# Evaluamos y verificamos la silueta para distintos valores de k y comaparamos con los resultados anteriores
max_k <- 10
silhouettes <- numeric(max_k)
for (k in 2:max_k) {
y_cluster <- kmeans(hawks_numeric, centers = k)$cluster
silhouettes[k] <- calcular_silueta(y_cluster, dist_hawks)
}
# Graficamos siluetas medias
plot(2:max_k, silhouettes[2:max_k], type = "b", pch = 19, frame = FALSE,
xlab = "Número de clusters (k)", ylab = "Silhouette media",
main = "Evaluación del número óptimo de clústers con método de Silhouette",
col = "tomato")
El gráfico muestra los valores de silueta media para diferentes números de clústers (k). El método de la silueta respalda que k = 2 es el número óptimo de clústers, con la mejor. Esto concuerda con las visualizaciones previas y la validez del modelo de K-means aplicado.
# Método Elbow
# Función para calcular la suma de cuadrados intra-cluster (WSS)
inercia_intracluster <- function(vdata, max_k) {
wss <- numeric(max_k)
for (k in 1:max_k) {
wss[k] <- sum(kmeans(vdata, centers = k)$withinss)
}
return(wss)
}
# Aplicamos la función a hawks_numeric
max_k <- 10
wss <- inercia_intracluster(hawks_numeric, max_k)
# Graficamos el método del codo
plot(1:max_k, wss, type = "b", pch = 19, col = "darkorange",
xlab = "Número de clusters (k)", ylab = "Suma de cuadrados intra-cluster (WSS)",
main = "Elbow Method para determinar el número óptimo de clústers")
El método del codo (Elbow Method) nos determinará el número óptimo de clústers en un análisis de agrupamiento como el K-means. Consiste en calcular la suma de cuadrados intra-clúster (WSS) para distintos valores de k, la variación total dentro de cada clúster. El punto de inflexión o codo es el número óptimo de clústers. El codo más claro parece estar en k=2 o k=3, lo que sugiere que esos valores son los más adecuados para segmentar los datos. Estos resultados coinciden con los obtenidos con los métodos anteriores.
Con el juego de datos proporcionado realiza un estudio aplicando DBSCAN y OPTICS, similar al del ejemplo 3
A diferencia de K-means, que tiende a formar grupos esféricos, DBSCAN y OPTICS identifican clústers de cualquier forma. Son adecuados cuando hay ruido, es decir presencia de outliers y formas complejas de agrupamiento, como es en nuestro caso. - minPts: número mínimo de puntos para formar una región dentro del radio eps como un punto central de clúster. Si minPts es bajo, se generarán más clústers y menos ruido, pero podrían ser poco consistentes. Si minPts es alto, solo zonas con mucha densidad formarán clústers, el resto será considerado como ruido. - eps (DBSCAN): radio de vecindad para formar un clúster. Si eps es un valor grande, el algoritmo agrupará todo en un único clúster. Si es pequeño, se crearán muchos grupos y habrían puntos quedarán fuera.
# Instalamos la biblioteca dbscan
if (!require('dbscan')) install.packages('dbscan')
library(dbscan)# Aplicamos DBSCAN con diferentes valores de eps
db_05 <- dbscan::dbscan(hawks_numeric, eps = 0.5, minPts = 10)
db_50 <- dbscan::dbscan(hawks_numeric, eps = 50, minPts = 10)
db_75 <- dbscan::dbscan(hawks_numeric, eps = 75, minPts = 10)
db_100 <- dbscan::dbscan(hawks_numeric, eps = 100, minPts = 10)
db_150 <- dbscan::dbscan(hawks_numeric, eps = 150, minPts = 10)
db_175 <- dbscan::dbscan(hawks_numeric, eps = 175, minPts = 10)
db_200 <- dbscan::dbscan(hawks_numeric, eps = 200, minPts = 10)
# Ver distribución de clústers
table(db_05$cluster)##
## 0
## 908table(db_50$cluster)##
## 0 1 2 3
## 41 552 63 252table(db_75$cluster)##
## 0 1 2 3
## 28 563 65 252table(db_100$cluster)##
## 0 1 2
## 18 567 323table(db_150$cluster)##
## 0 1 2
## 6 575 327table(db_175$cluster)##
## 0 1
## 4 904table(db_200$cluster)##
## 0 1
## 3 905eps = 0.5 Todos los puntos se consideran ruido. No hay agrupaciones válidas.
eps = 50–75 Comienza a detectar estructura real en los datos, identificando 3 clústers principales y reduciendo el ruido a 41.
eps = 100–150: Los clústers se fusionan, reduciendo el número total a 2 grandes clústers. Ruido casi nulo.
eps > 175: Prácticamente todos los puntos se agrupan en un único clúster, lo cual indica que eps es demasiado alto y se pierde la diferenciación.
# Aplicamos OPTICS con eps = 150 y minPts = 10
opt_res <- optics(hawks_numeric, eps = 150, minPts = 10)
# Visualizamos la distancia de alcanzabilidad
plot(opt_res, main = "Reachability plot - OPTICS eps = 150")
# Extraemos los clústers usando una altura de corte eps_cl
opt_clusters <- extractDBSCAN(opt_res, eps_cl = 150)
# Mostramos los clústers con líneas convexas
hullplot(hawks_numeric, opt_clusters, main = "Clustering OPTICS con eps_cl = 150")
# Aplicamos OPTICS con eps = 100 y minPts = 10
opt_res <- optics(hawks_numeric, eps = 100, minPts = 10)
# Visualizamos la distancia de alcanzabilidad
plot(opt_res, main = "Reachability plot - OPTICS eps = 100")
# Extraemos los clústers usando una altura de corte eps_cl
opt_clusters <- extractDBSCAN(opt_res, eps_cl = 100)
# Mostramos los clústers con líneas convexas
hullplot(hawks_numeric, opt_clusters, main = "Clustering OPTICS con eps_cl = 100")
En este análisis se aplicó el algoritmo OPTICS al dataset
hawks_numeric para descubrir agrupaciones basadas en densidad. El
parámetro eps = 150 define la distancia máxima considerada para formar
una vecindad, y minPts = 10 indica el número mínimo de puntos que deben
encontrarse dentro de esa distancia para considerar que hay una densidad
suficiente para formar un clúster.
Una vez aplicado OPTICS aplicado con eps = 150 generó el Reachability
Plot, donde se observan varios valles profundos, indicativos de regiones
densas (clústers). Los Clústers extraídos con eps_cl = 150 usando
extractDBSCAN() fueron dos clústers principales más una serie de puntos
dispersos considerados ruido, representados como puntos negros. En el
Reachability Plot, representan transiciones entre grupos o puntos poco
accesibles, mientras que los valles amplios y profundos representan
agrupaciones consistentes. En el gráfico de hullplot, se visualizan
claramente dos grandes clústers convexos bien separados, con unos pocos
outliers fuera de ellos. Comparando con eps_cl = 100, que genera también
dos clústers pero con mayor cantidad de puntos considerados como ruido,
el valor de 150 logra una mejor agrupación, se reducen los outliers
manteniendo una separación definida entre grupos. Por lo tanto, eps_cl =
150 parece ser el valor más adecuado para esta estructura de datos.
# OPTICS sin definir eps_cl, solo minPts
opt_res <- optics(hawks_numeric, minPts = 10)
# Extraemos los clusters con el método jerárquico usando xi
opt_xi <- extractXi(opt_res, xi = 0.05)
# Resultado
print(opt_xi)## OPTICS ordering/clustering for 908 objects.
## Parameters: minPts = 10, eps = 532.082042546072, eps_cl = NA, xi = 0.05
## The clustering contains 13 cluster(s) and 1 noise points.
##
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
## eps_cl, xi, clusters_xi, cluster# Visualización del reachability plot con los clústers generados por xi
plot(opt_xi, main = "Reachability plot con xi = 0.05")
# Graficamos los clústers detectados (polígonos convexos)
hullplot(hawks_numeric, opt_xi, main = "Clústers detectados con extractXi (xi = 0.05)")
Explorando la agrupación jerárquica a través del parámetro xi = 0.05. De esta forma permite detectar automáticamente los clústers observando los cambios relativos a la densidad de datos, en lugar de fijarlo manualmente un umbral como eps_cl. Reachability Plot con xi = 0.05, muestra múltiples valles de baja distancia de alcance, lo que indica la presencia clusters densos. Las áreas que se elevan las barras corresponden a zonas menos densas o ruido. Se identificaron 13 clusters y 1 punto de ruido. Hullplot representa los clústers con polígonos convexos en el espacio de componentes principales (PC1 y PC2).Cada color representa un clúster distinto.
# opt_res <- optics(hawks_numeric, minPts = 10)
# Lista de valores de xi a probar
xi_values <- c(0.01, 0.03, 0.05, 0.1, 0.15)
# Aplicamos extractXi para cada xi y visualizamos los resultados
for (xi in xi_values) {
cat(">>> xi =", xi, "\n")
# Extraer clústers con el valor actual de xi
opt_xi <- extractXi(opt_res, xi = xi)
# Imprimir resumen del resultado
print(opt_xi)
# Graficar reachability plot con clústers detectados
plot(opt_xi, main = paste("Reachability plot con xi =", xi))
# Graficar los clústers con polígonos convexos
hullplot(hawks_numeric, opt_xi, main = paste("Clústers detectados con extractXi (xi =", xi, ")"))
}## >>> xi = 0.01
## OPTICS ordering/clustering for 908 objects.
## Parameters: minPts = 10, eps = 532.082042546072, eps_cl = NA, xi = 0.01
## The clustering contains 33 cluster(s) and 1 noise points.
##
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
## eps_cl, xi, clusters_xi, cluster
## >>> xi = 0.03
## OPTICS ordering/clustering for 908 objects.
## Parameters: minPts = 10, eps = 532.082042546072, eps_cl = NA, xi = 0.03
## The clustering contains 18 cluster(s) and 1 noise points.
##
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
## eps_cl, xi, clusters_xi, cluster
## >>> xi = 0.05
## OPTICS ordering/clustering for 908 objects.
## Parameters: minPts = 10, eps = 532.082042546072, eps_cl = NA, xi = 0.05
## The clustering contains 13 cluster(s) and 1 noise points.
##
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
## eps_cl, xi, clusters_xi, cluster
## >>> xi = 0.1
## OPTICS ordering/clustering for 908 objects.
## Parameters: minPts = 10, eps = 532.082042546072, eps_cl = NA, xi = 0.1
## The clustering contains 8 cluster(s) and 1 noise points.
##
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
## eps_cl, xi, clusters_xi, cluster
## >>> xi = 0.15
## OPTICS ordering/clustering for 908 objects.
## Parameters: minPts = 10, eps = 532.082042546072, eps_cl = NA, xi = 0.15
## The clustering contains 5 cluster(s) and 1 noise points.
##
## Available fields: order, reachdist, coredist, predecessor, minPts, eps,
## eps_cl, xi, clusters_xi, cluster
xi = 0.01: Alta sensibilidad. El algoritmo detecta muchos clústers, incluso pequeños grupos locales. Es útil para capturar estructuras muy finas pero puede producir sobreajuste (demasiados clústers no significativos). xi = 0.05 (valor común). Balance entre sensibilidad y generalización. Aquí detectas 13 clústers, lo que probablemente representa subgrupos reales con menor ruido. Es un valor recomendado para exploraciones iniciales. xi = 0.10 o 0.15: Baja sensibilidad. Se detectan solo clústers grandes y compactos. Se pierde detalle, pero puede ser útil para una vista más general o si se desea evitar fragmentación.
Realiza una comparativa de los métodos k-means y DBSCAN
Para realizar el análisis del conjunto de datos Hawks, se aplicó los métodos de clustering: K-means y DBSCAN, cada uno con supuestos y objetivos distintos.
K-means se aplicado sobre cinco variables numéricas (Wing, Weight, Culmen, Hallux, Tail) con datos reales. El método exige especificar el número de clústers, por lo que utilizamos el método de la silueta y el método del codo o Elbow para determinar el valor de k. Se observó que k=3 ofrecía una buena separación entre grupos, con una silueta media de 0.765. Los clústers generados se correspondían parcialmente con las especies de halcones, diferenciando claramente al halcón de cola roja (Red-Tailed) debido a sus mayores dimensiones y peso. Sin embargo, todos los puntos fueron clasificados, no se trataron los valores atípicos, lo que puede llevar a una agrupación forzada en presencia de outliers.
En contraste, DBSCAN, basado en densidades, no requiere conocer k previamente. Se aplicó con distintos valores de eps y minPts, y se observó que con eps=150 y minPts=10 el algoritmo identificaba dos grupos principales y algunos puntos como ruido. Esta capacidad para detectar outliers es una de las fortalezas de DBSCAN. Además, los clústers obtenidos reflejaban estructuras de densidad más naturales y no dependían de formas geométricas estrictas como en K-means. El algoritmo capturó las agrupaciones principales basadas en medidas físicas sin imponer una estructura, lo cual se visualizó claramente en los gráficos hullplot y en el reachability plot del análisis OPTICS.
Conclusión Mientras que K-means proporciona una clasificación eficiente cuando se conoce la estructura o número de grupos, DBSCAN resulta más robusto ante ruido, y más flexible ante clústers de forma irregular, tampoco requiere especificar k, adaptándose mejor a datos como Hawks, donde existen especies con tamaños bien diferenciados y posibles valores extremos.