library("purrr")
library("dplyr")
library("fastDummies")
library("ggplot2")
library("dendextend")
library("cluster")
library("tidyr")
library("factoextra")
library("gridExtra")
# Datasets
lineup <- readRDS("lineup.rds")
oes <- readRDS("oes.rds")
ws_customers <- readRDS("ws_customers.rds")
df <- USArrests

Conjunto de técnicasmediante la cual observaciones son divididas en subgrupos que comparte características en común.

Pasos de clustering

  1. Asegurarse que los datos estén listos para ser clusterizados (Sin datos datos faltantes , variables en escalas similares)
  2. Métrica de similitud apropiada para la clusterización
  3. Método de clusterización deseado para

Distancia entre observaciones

\[ distancias = 1 - similitud \]

La medida de distancia más común es la euclidiana

\[ d(p,q) = \sqrt{\sum_{i=1}^n{(q_{i}- p_{i})^2}} \] Distancia manhattan:

\[ d_{man} = \sum_{i=1}^n|(x_{i} - y_{i})| \]

Existen también medidas de distancias basadas en correlaciones. Las cuales son obtenidas restando la correlación a 1.

Distancia de correlación de perason

\[ d_{cor} = 1- \frac{\sum_{i=1}^n(x_{i} - \bar{x})(y_{i}-\bar{y})}{\sqrt{\sum_{i=1}^n(x_{i}-\bar{x})^2\sum_{i=1}^n(y_{i}-\bar{y})^2}} \]

Distancia de correlación de spearman

Calcula la correlación entre los rangos de \(x\) e \(y\) \[ d_{cor} = 1- \frac{\sum_{i=1}^n(x_{i}^{'} - \bar{x}^{'})(y_{i}^{'}-\bar{y}^{'})}{\sqrt{\sum_{i=1}^n(x_{i}^{'}-\bar{x}^{'})^2\sum_{i=1}^n(y_{i}^{'}-\bar{y}^{'})^2}} \]

Donde \(x_{i}^{'} = rank(x_{i})\) y \(y_{i}^{'} = rank(y_{i})\).

Distancia de correlación de kendall

Mide la correspondencia entre el ranking de las variabesl \(x\) e \(y\). El número de pares posibles de las variables es \(n(n-1)/2\), donde \(n\) es el tamaño de \(x\) e \(y\). Empezamos ordenando los pares por los valores de \(x\). Si \(x\) y \(y\) están correlacionados, entonces tendrían los mismos ordenes de ranks relativos. Ahora, para cada \(y_{i}\), contamos el número de \(y_{j}\) > \(y_{i}\) (pares concordantes (c)) y el número de \(y_{j}\) < \(y_{i}\) (pares discordantes (d))

La distancia de correlación de kendall está definida como:

\[ d_{kend}(x, y) = 1 - \frac{n_{c} - n_{d}}{\frac{1}{2}n(n-1)} \]

Es muy común trabajr con distancia euclidiana, pedo depende de los datos y el tipo de investigación qué medida es la más adecuada de usar

df <- na.omit(df)
df <- scale(df)
distance <- get_dist(df)
fviz_dist(distance, gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))

dist(lineup, method = 'euclidean')
##            1         2         3         4         5         6         7
## 2   4.123106                                                            
## 3  10.295630 13.453624                                                  
## 4  12.041595 10.295630 14.035669                                        
## 5  13.038405 14.866069 20.099751 24.839485                              
## 6  14.035669 13.341664 23.769729 23.409400  8.544004                    
## 7  16.278821 13.038405 26.400758 20.099751 18.027756 10.198039          
## 8  21.931712 25.495098 12.206556 25.298221 28.160256 34.000000 38.209946
## 9  22.022716 23.537205 13.601471 17.204651 33.541020 36.055513 36.055513
## 10 20.615528 18.439089 21.377558  8.602325 33.241540 30.886890 25.495098
## 11 24.000000 23.345235 33.376639 33.241540 14.764823 10.049876 16.278821
## 12 27.018512 28.160256 18.973666 20.615528 38.832976 41.000000 40.261644
##            8         9        10        11
## 2                                         
## 3                                         
## 4                                         
## 5                                         
## 6                                         
## 7                                         
## 8                                         
## 9  15.231546                              
## 10 31.144823 19.235384                    
## 11 42.720019 46.010868 40.311289          
## 12 19.416488  5.385165 20.518285 51.009803

Reescalamiento de variables

La estandarización es una forma común de escalar las variables. Con esto aseguramos que ahora las variables tienen media cero y desviación estandar 1.

\[ variables_{escalada} = \frac{x - \mu}{\sigma} \]

scale(lineup)
##                x          y
##  [1,] -0.1745532  0.1403885
##  [2,] -0.2380271 -0.3409435
##  [3,]  0.3967118  0.7420536
##  [4,]  0.3332379 -0.9426086
##  [5,] -0.8727660  0.9827196
##  [6,] -1.0631877  0.0200555
##  [7,] -0.9362399 -1.1832747
##  [8,]  0.8410291  1.9453838
##  [9,]  1.2218724  0.2607215
## [10,]  0.6506074 -1.7849398
## [11,] -1.6979266  0.1403885
## [12,]  1.5392419  0.0200555
## attr(,"scaled:center")
##          x          y 
##  1.7500000 -0.1666667 
## attr(,"scaled:scale")
##         x         y 
## 15.754509  8.310271

Pero esto funciona bien para variables continuas. Hay otra forma de hacer para variables categoricas:

Podemos conseguir medir la distancia con el índice de Jaccard:

\[ J(A, B) = \frac{A\cap B}{A\cup B} \] Entre más cercano a 1 más similud hay entre la observación A y B. Entonces la distancia sería:

\[ distancia(A, B) = 1 - J(A, B) \]

drinks <- tibble(wine    = c(TRUE, FALSE, TRUE),
                 beer    = c(TRUE, TRUE, FALSE),
                 whiskey = c(FALSE, TRUE, TRUE),
                 vodka   = c(FALSE, TRUE, FALSE))

dist(drinks, method = "binary")
##           1         2
## 2 0.7500000          
## 3 0.6666667 0.7500000

Si se desea realizar el procedimiento con variables con multiples categorías se deben crear variables dummy para cada categoría:

df <- tibble(color = c("red", "green", "blue", "blue"),
             sport = c("soccer", "hockey", "hockey", "soccer"))

df_dummified <- dummy_cols(df)
print(df_dummified)
## # A tibble: 4 x 7
##   color sport  color_blue color_green color_red sport_hockey sport_soccer
##   <chr> <chr>       <int>       <int>     <int>        <int>        <int>
## 1 red   soccer          0           0         1            0            1
## 2 green hockey          0           1         0            1            0
## 3 blue  hockey          1           0         0            1            0
## 4 blue  soccer          1           0         0            0            1
dist(df_dummified[, 3:ncol(df_dummified)], method = 'binary')
##           1         2         3
## 2 1.0000000                    
## 3 1.0000000 0.6666667          
## 4 0.6666667 1.0000000 0.6666667

Cluster Jerarquico

Este método de manera iterativa agrupa observaciones basado en sus distancias hasta que cada observación pertenece a un grupo más grande.

No requiere especificación previa de clusters, como si lo hace k-means

El cluster jerarquico puede ser dividido en dos tipos

Criterios en enlazamiento (Linkage criteria)

Una vez tengamos las distancias entre observaciones necesitamos algún método para juntar las oservaciones, de manera que sean lo más parecidas entre ellas y muy distintas a las observaciones de los demás grupos (clusters). En otras palabras, necesitamos un método que seleccione a una observación y decida que es la más cercana a un grupo en partícular.

  • Complete linkage: máxima distancia entre dos conjuntos

Ejemplo: Supongamos que tenemos 4 observaciones. La obs. 1 y 2 tienen la mínima distancia y armamos el primer cluster. Para seleccionar la observación que se anexará al cluster buscamos a la más cercana a nuestro primer cluster, nos queda entonces revisar entre la observación 3 y 4. comparamos las distancias entre la observación 3 y cada una de las observaciones del cluster y lo mismo para la observación 4. Las distancas encontradas entre la observación candidata y las observaciones dentro del cluster son revisadas, de estas seleccionamos la máxima y finalmente comparamos el mínimo de esta distancias “máxima” entre las observaciones candidatas para incluirla dentro del cluster.

  • Single linkage: mínima distancia entre dos conjuntos
  • Average linkage: Distancia promedio entre dos conjuntos

Estos dos ultimos siguen la misma lógica, solo que con mínimo y promedio respectivamente, para seleccionar la observación que se anexará al cluster.

  • Centroid linkage: se computa la disimilitud entre el centroide para el cluster 1 (un vector de media de longitud \(p\) variables) y el centroide para el cluster 2.

  • Método de mínima varianza de ward: mínimiza la varianza total intracluster. En cada paso el par de clusters con mínima varianza intracluster es unido.

# Ejemplo de distancias entre jugadores en un campo de futbol

ggplot(lineup, aes(x = x, y = y)) + geom_point() + theme_classic()

# distancia euclidiana
dist_players <- dist(lineup, method = 'euclidean')

# Cluster jerarquico con el método complete
hc_player <- hclust(dist_players, method = 'complete')

# Determinamos a dónde pertenece cada observación

cluster_assigments <- cutree(hc_player, k = 2)

# asignamos los clusters

assigned_cluster <- lineup %>% mutate(cluster = as.factor(cluster_assigments))


ggplot(assigned_cluster, aes(x = x, y = y, color = cluster)) +
  geom_point() + theme_classic()

fviz_cluster(list(data = lineup, cluster = cluster_assigments))

Dendograma

  • height hace referencia a la distancia total de las observaciones en el cluster
plot(hc_player, cex = 0.6)
rect.hclust(hc_player, k = 2, border = 2:5)

Ahora, podemos revisar el impacto que tiene cada método de enlazamiento en el dendograma, y así con la intuición y criterio experto se puede tomar la decisión de cuál manejar.

hc_complete <- hclust(dist_players, method = 'complete')
hc_single <- hclust(dist_players, method = 'single')
hc_average <- hclust(dist_players, method = 'average')


par(mfrow = c(1,3))
plot(hc_complete, main = 'Complete Linkage')
plot(hc_single, main = 'Single Linkage')
plot(hc_average, main = 'Average Linkage')

Para agregar mayor claridad al dendograma podemos usar la librería dendextend para poder visualizar de forma más colorida los dendogramas, basados en los criterio de k clusters o de la distancia en la consideramos que deberían empezar a formarse los clusters.

# La interpretación acá sería:
# Cluster jerarquico con distancia maxima euclidiana donde la distancia de los individuos dentro no puede ser superior a 15
dend_players <- as.dendrogram(hc_player)
dend_colored <- color_branches(dend_players, h = 15)
plot(dend_colored)

dend_colored <- color_branches(dend_players, k = 3)
plot(dend_colored)

Si lo desearamos también podríamos cambiar el parámetro k de la función cutree() por h, de manera que los clusters los haga con un corte en la distancia (euclidiana en este caso) determinada por la altura del dendograma (height).

A medida que h es más bajo más clusters estamos permitiendo que se creen. Esto dependerá según la tarea que estemos desempeñando qué entendemos por grupos homogeneos.

  • Cuando se tienen multiples variables a clusterizar una buena idea es usa PCA para poder gráficarlas y de esa manera marcarles el cluster al que pertenecen y así revisar si los grupos están siendo bien clusterizados (Qué tan preciso pueda llegar a ser el gráfico depende también de la varianza recogida por los primeros dos componentes, que son los que se grafican)

  • También se podría usar un gráfico bidimensional por cada par de variables, pero esto deja de ser eficiente cuando las dimensiones crecen demasiado o no se tiene una intuición sobre las variables más representativas del dataset.

  • Análisis descriptivo de las variables también puede dar buena información sobre cómo las variables están asociadas entre si

Otras funciones alternativas para cluster jerarquico en R son agnes y diana (AGNES para agglomerative nesting y DIANA para divisive analysis).

la función agnes es similar a hclust, pero también arroja el coeficiente de aglomeración que mide el monto de estructura de cluster encontrado (valores cercanos a 1 sugieren una estructura de cluster fuerte)

agnes_hc <- agnes(lineup, method = "complete")

agnes_hc$ac
## [1] 0.822312

Esto nos puede ayudar a encontrar métodos de clustering que identifican estructuras mas fuertes.

# methods to assess
m <- c( "average", "single", "complete", "ward")
names(m) <- c( "average", "single", "complete", "ward")

# function to compute coefficient
ac <- function(x) {
  agnes(lineup, method = x)$ac
}

map_dbl(m, ac)
##   average    single  complete      ward 
## 0.6863075 0.4114650 0.8223120 0.8553583

La función diana es similar.

# compute divisive hierarchical clustering
hc4 <- diana(lineup)

# Divise coefficient; amount of clustering structure found
hc4$dc
## [1] 0.822312
## [1] 0.8514345

# plot dendrogram
pltree(hc4, cex = 0.6, hang = -1, main = "Dendrogram of diana")

Podemos comparar también dos dendogramas (Por ejemplo para ver cómo cambian debido a métodos de enlazamiento)

También nos provee de una medida de entrelazamiento: valores bajos (cercanos a cero) de entrelazamiento indican buena alineación

# Compute distance matrix
res.dist <- dist(lineup, method = "euclidean")

# Compute 2 hierarchical clusterings
hc1 <- hclust(res.dist, method = "complete")
hc2 <- hclust(res.dist, method = "ward.D2")

# Create two dendrograms
dend1 <- as.dendrogram (hc1)
dend2 <- as.dendrogram (hc2)

tanglegram(dend1, dend2)

dend_list <- dendlist(dend1, dend2)
entanglement(dend_list)
## [1] 0.04123658

Podemos obtener el número óptimo de clusters con el método de elbow.

fviz_nbclust(lineup, FUN = hcut, method = "wss")

Caso: gasto de clientes en el super mercado

Revisemos si hay patrones que distingan a los clientes en una especie de “perfiles” basados en sus patrones de consumo

customers_spend <- ws_customers

# Revisemos las descriptivas
# Las escalas de medida son en dinero, así que no requieren preprocesamiento
summary(customers_spend)
##       Milk          Grocery          Frozen     
##  Min.   :  333   Min.   : 1330   Min.   :  264  
##  1st Qu.: 1375   1st Qu.: 2743   1st Qu.:  824  
##  Median : 2335   Median : 5332   Median : 1455  
##  Mean   : 4831   Mean   : 7830   Mean   : 2870  
##  3rd Qu.: 5302   3rd Qu.:10790   3rd Qu.: 3046  
##  Max.   :25071   Max.   :26839   Max.   :15601
# distancia euclidiana
dist_customers <- dist(customers_spend, method = "euclidean")
# Método de enlazamiento cmpleto
hc_customers <- hclust(dist_customers, method = 'complete')

# Cortamos los cluster en grupos cuya distancia no sea mayor a 15000
clust_customers <- cutree(hc_customers, h = 15000)

segment_customers <- mutate(customers_spend, cluster = clust_customers)

# Contamos el número de cliente que hay por cluster
count(segment_customers,cluster)
## # A tibble: 4 x 2
##   cluster     n
##     <int> <int>
## 1       1     5
## 2       2    29
## 3       3     5
## 4       4     6
# PArametros del dentrograma
dend_customers <- as.dendrogram(hc_customers)
dend_colored <- color_branches(dend_customers, h =15000)

# Dendograma 

plot(dend_colored)

# Revisamos los gastos promedio de cada grupo por producto
segment_customers %>% 
  group_by(cluster) %>% 
  summarise_all(list(mean))
## # A tibble: 4 x 4
##   cluster   Milk Grocery Frozen
##     <int>  <dbl>   <dbl>  <dbl>
## 1       1 16950   12891.   991.
## 2       2  2513.   5229.  1796.
## 3       3 10452.  22551.  1355.
## 4       4  1250.   3917. 10889.

Cluster de K-means

Consiste en definir clusters tal que la variación total intracluster es minimizada.Hay multiples algoritmos par alcanzar este objetivo, el estandar para k-means es la función objetivo que define la suma de las distancias euclidianas al cuadrado entre los elementos del centroide correspondiente:

\[ W(C_{k}) = \sum_{x_{i}\in C_{k}}{(x_{i} - \mu_{k})^2} \] Donde \(x_{i}\) es un punto de datos que pertenece al cluster \(C_{k}\) y \(\mu_{k}\) es el valor medio de los datos asignados al cluster \(C_{k}\)

Definimos la suma total instracluster como:

\[ SS = \sum_{k=1}^K W(C_{k}) = \sum_{k=1}^K \sum_{x_{i}\in C_{k}}{(x_{i} - \mu_{k})^2} \]

Pasos k-means

Establecemos un número k de clusters Creamos k centroides aleatorios (promedios de cada variable ubicado en un espacio p-dimensional, donde p es el número de variables) * Se establece distancias entre cada observación y cada centroida y se le asigna un cluster basado en la mínima distancia al centroido * Se calculan de nuevo los centroides con la nueva iteración de observaciones * Nos detenemos hasta que cada cluster se estabilice

kmeans_cluster  <- kmeans(lineup, centers = 2)

kmeans_cluster
## K-means clustering with 2 clusters of sizes 6, 6
## 
## Cluster means:
##           x          y
## 1 -11.33333 -0.5000000
## 2  14.83333  0.1666667
## 
## Clustering vector:
##  [1] 1 1 2 2 1 1 1 2 2 2 1 2
## 
## Within cluster sum of squares by cluster:
## [1] 570.8333 863.6667
##  (between_SS / total_SS =  58.9 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"
kmeans_cluster$cluster
##  [1] 1 1 2 2 1 1 1 2 2 2 1 2
kmeans_lineup <- lineup %>% mutate(cluster = kmeans_cluster$cluster)

ggplot(kmeans_lineup, aes(x = x, y = y, color = factor(cluster))) + geom_point()

# Otra opción

fviz_cluster(kmeans_cluster, data = lineup)

# Con etiquetas

lineup %>%
  as_tibble() %>%
  mutate(cluster = kmeans_cluster$cluster,
         state = row.names(lineup)) %>%
  ggplot(aes(x, y, color = factor(cluster), label =  row.names(lineup))) +
  geom_text()

k2 <- kmeans(lineup, centers = 2, nstart = 25)
k3 <- kmeans(lineup, centers = 3, nstart = 25)
k4 <- kmeans(lineup, centers = 4, nstart = 25)
k5 <- kmeans(lineup, centers = 5, nstart = 25)

# plots to compare
p1 <- fviz_cluster(k2, geom = "point",  data = lineup) + ggtitle("k = 2")
p2 <- fviz_cluster(k3, geom = "point",  data = lineup) + ggtitle("k = 3")
p3 <- fviz_cluster(k4, geom = "point",  data = lineup) + ggtitle("k = 4")
p4 <- fviz_cluster(k5, geom = "point",  data = lineup) + ggtitle("k = 5")


grid.arrange(p1, p2, p3, p4, nrow = 2)

Evaluación del número adecuado de K

  • La primera opción siempre es el contexto y criterio experto ante la pregunta que se desea responder.

  • Otra opción más técnica es el “elbow method” o método de codo: La idea es calcular la suma cuadrada intracluster para cada cluster, que es la suma de las distancias euclidianas entre cada observación y el centroide correspondiente al cluster al cual dicha obervación está asignada. Hacemos este cálculo para cada k valores de cluster y gráficamos, lo cual tiene una forma de codo. Haremos la elección de un K que disminuya significativamente el suma total intracluster

# Podemos hacer la prueba para k = 1 hasta n, donde n es el # de observaciones

k_clusters <- nrow(lineup) - 1

within_ss <- vector(mode = "numeric", k_clusters)


for (i in seq_along(within_ss)) {
  
  kmeans_iteration <- kmeans(lineup, centers = i)
  within_ss[i]     <- kmeans_iteration$tot.withinss
  
}


qplot(1:11, within_ss, geom = "point") +
  geom_line() + theme_classic() +
  xlab("K")

# Podemos hacer lo mismo con programación funcional

within_ss <- map_dbl(1:11, function(k){
  
  kmeans_iteration <- kmeans(lineup, centers = k)
  kmeans_iteration$tot.withinss
  
})

tibble(k = 1:11, Within_SS = within_ss)
## # A tibble: 11 x 2
##        k Within_SS
##    <int>     <dbl>
##  1     1    3490. 
##  2     2    1434. 
##  3     3     881. 
##  4     4     622. 
##  5     5     482. 
##  6     6     265. 
##  7     7     171  
##  8     8     134  
##  9     9      59.5
## 10    10      83  
## 11    11      14.5

Análisis de silueta

Es usado para deteminar qué tan bien las observaciones se ajustan al cluster correspondiente y puede ser asumido como un método adicional para estimar el valor de K.

Ancho de l silueta es la medidad usada para cada observacion:

  • Distancia intracluster \(C_{i}\): Distancia euclidiana promedio de la observación a cada otra observación en el mismo cluster
  • Distancia al vecino más cercano \(N_{i}\): distancia promedio de la observación a los puntos del cluster vecino más cercano. La distancia promedio de los otros cluster más pequeña a la observación es usada para determinar el vecino más cercano.

\[ S_{i} = 1 - \frac{C_i}{N_{i}}, \ si \ \ C_{i} < N_{i} \\ S_{i} = 0, \ \ \ \ \ \ \ \ \ \ \ \ si \ \ C_{i} = N{i} \\ S_{i} = \frac{N_{i}}{C_{i}}, \ \ \ \ \ \ \ \ \ si \ \ C_{i} > N_{i} \]

Esto quiere decir, que valores cercanos 1 indican que la observaci´pn encaja muy bien dentro del cluster y valores cercanos a cero indican que una obserbación podría estar en el borde de dos clusters y posiblemente pertenezca a uno de los dos.

Valores de -1 o cercanos a él indican que una observación pertenecería al cluster vecino más cercano

s_i <- pam(lineup, k = 3)

s_i$silinfo$widths
##    cluster neighbor    sil_width
## 4        1        2  0.465320054
## 2        1        3  0.321729341
## 10       1        2  0.311385893
## 1        1        3  0.271890169
## 9        2        1  0.443606497
## 8        2        1  0.398547473
## 12       2        1  0.393982685
## 3        2        1 -0.009151755
## 11       3        1  0.546797052
## 6        3        1  0.529967901
## 5        3        1  0.359014657
## 7        3        1  0.207878188
# Podemos visualizarlas también

plot(s_i)

#

s_i$silinfo$avg.width
## [1] 0.353414
silhouette_width <- map_dbl(2:11, function(k){
  
  silhouette_k <- pam(lineup, k = k)
  silhouette_k$silinfo$avg.width
  
})

tibble(k = 2:11, silhouette_width = silhouette_width)
## # A tibble: 10 x 2
##        k silhouette_width
##    <int>            <dbl>
##  1     2           0.416 
##  2     3           0.353 
##  3     4           0.354 
##  4     5           0.372 
##  5     6           0.344 
##  6     7           0.324 
##  7     8           0.328 
##  8     9           0.255 
##  9    10           0.210 
## 10    11           0.0999
qplot(2:11, silhouette_width, geom = "point") +
  geom_line() + theme_classic() +
  xlab("K")

Segmentación de mercado: aplicación cluster k-emans

# Usamos los datos de consumo por productos
sil_width <- map_dbl(2:10,  function(k){
  model <- pam(x = customers_spend, k = k)
  model$silinfo$avg.width
})


sil_df <- data.frame(
  k = 2:10,
  sil_width = sil_width
)


ggplot(sil_df, aes(x = k, y = sil_width)) +
  geom_line() +
  scale_x_continuous(breaks = 2:10) + geom_point() + theme_classic()

Revisamos los perfiles

set.seed(42)


model_customers <- kmeans(customers_spend, centers = 2)


clust_customers <- model_customers$cluster


segment_customers <- mutate(customers_spend, cluster = clust_customers)


count(segment_customers, cluster)
## # A tibble: 2 x 2
##   cluster     n
##     <int> <int>
## 1       1    35
## 2       2    10
segment_customers %>% 
  group_by(cluster) %>% 
  summarise_all(list(mean))
## # A tibble: 2 x 4
##   cluster   Milk Grocery Frozen
##     <int>  <dbl>   <dbl>  <dbl>
## 1       1  2296.    5004  3354.
## 2       2 13701.   17721  1173

Comparación entre los dos métodos de cluster

Criterio Cluster Jerarquico K-means
Distancia usada: Cualquiera Solo euclidiana
Resultados estables: SI NO
Evaluación del número de clusters: Dendograma , Silueta, Codo Silueta, Codo, (PCA tal vez)
Complejidad en la computación: relativamente alta Relativamente baja

Otros métodos populares

  • K-medoides
  • DBScan
  • Optics

Caso de estudio: Salarios anuales por tipo de ocupación

En este dataset reviamos si hay patrones entre los salarios medios anuales por ocupación y probaremos tanto cluster jerarquico omo kmeans

# Calculate Euclidean distance between the occupations
dist_oes <- dist(oes, method = 'euclidean')

# Generate an average linkage analysis 
hc_oes <- hclust(dist_oes, method = 'average')

# Create a dendrogram object from the hclust variable
dend_oes <- as.dendrogram(hc_oes)

# Plot the dendrogram
plot(dend_oes)

# Color branches by cluster formed from the cut at a height of 100000
dend_colored <- color_branches(dend_oes, h = 100000)

# Plot the colored dendrogram
plot(dend_colored)

dist_oes <- dist(oes, method = 'euclidean')
hc_oes <- hclust(dist_oes, method = 'average')

# Use rownames_to_column to move the rownames into a column of the data frame
df_oes <- tibble::rownames_to_column(as.data.frame(oes), var = 'occupation')

# Create a cluster assignment vector at h = 100,000
cut_oes <- cutree(hc_oes, h = 100000)

# Generate the segmented the oes data frame
clust_oes <- mutate(df_oes, cluster = cut_oes)

# Create a tidy data frame by gathering the year and values into two columns
gathered_oes <- gather(data = clust_oes, 
                       key = year, 
                       value = mean_salary, 
                       -occupation, -cluster)
# View the clustering assignments by sorting the cluster assignment vector
sort(cut_oes)
##                 Management                      Legal 
##                          1                          1 
##        Business Operations           Computer Science 
##                          2                          2 
##   Architecture/Engineering  Life/Physical/Social Sci. 
##                          2                          2 
##   Healthcare Practitioners         Community Services 
##                          2                          3 
## Education/Training/Library  Arts/Design/Entertainment 
##                          3                          3 
##         Healthcare Support         Protective Service 
##                          3                          3 
##           Food Preparation  Grounds Cleaning & Maint. 
##                          3                          3 
##              Personal Care                      Sales 
##                          3                          3 
##      Office Administrative   Farming/Fishing/Forestry 
##                          3                          3 
##               Construction Installation/Repair/Maint. 
##                          3                          3 
##                 Production      Transportation/Moving 
##                          3                          3
# Plot the relationship between mean_salary and year and color the lines by the assigned cluster
ggplot(gathered_oes, aes(x = year, y = mean_salary, color = factor(cluster))) + 
    geom_line(aes(group = occupation))

K-means clustering

# Use map_dbl to run many models with varying value of k (centers)
tot_withinss <- map_dbl(1:10,  function(k){
  model <- kmeans(x = oes, centers = k)
  model$tot.withinss
})

# Generate a data frame containing both k and tot_withinss
elbow_df <- data.frame(
  k = 1:10,
  tot_withinss = tot_withinss
)

# Plot the elbow plot
ggplot(elbow_df, aes(x = k, y = tot_withinss)) +
  geom_line() +
  scale_x_continuous(breaks = 1:10)

# Use map_dbl to run many models with varying value of k
sil_width <- map_dbl(2:10,  function(k){
  model <- pam(oes, k = k)
  model$silinfo$avg.width
})

# Generate a data frame containing both k and sil_width
sil_df <- data.frame(
  k = 2:10,
  sil_width = sil_width
)

# Plot the relationship between k and sil_width
ggplot(sil_df, aes(x = k, y = sil_width)) +
  geom_line() +
  scale_x_continuous(breaks = 2:10)