k-means: Base algébrica

La minimización de la traza de la Matriz de dispersión: k-Means

k-means es un algoritmo que crea clusters a través de la minimización de la traza de la Matriz de dispersión \(\mathbf{W}\) de los datos.

Vector columna y su transpuesta

Un vector columna es una matriz de dimensiones \(n\times 1\) que organiza \(n\) componentes numéricos en una sola columna vertical. \[v=\left[\begin{matrix}x_{1}\\ x_{2}\\ \vdots \\ x_{n}\end{matrix}\right]\]

Representa un conjunto ordenado utilizado en álgebra lineal, donde la transpuesta de un vector columna es un vector fila \(1\times n\).

Medición de la variación dentro de un cluster

Si representamos el vector de medias de un cluster i-ésimo, \[ \overline x_i\; \text o\; m_k= \frac{1}{n} \sum_{j=1}^{n_i} x_{ij} \]

Si representamos la distancia (error) entre un punto de un cluster \(k\) y su centroide (media) con un vector columna \((x-m_{k})\), al multiplicarlo por su transpuesta \((x-m_{k})^{T}\) (un vector fila), realizamos un producto externo donde, en este caso hipotizando un vector columna de 3 elementos:

Vector columna: \((x-m_{k})=\left[\begin{matrix}\Delta x_{1}\\ \Delta x_{2}\\ \Delta x_{3}\end{matrix}\right]\)
Vector fila (transpuesta): \((x-m_{k})^{T}=\left[\begin{matrix}\Delta x_{1}&\Delta x_{2}&\Delta x_{3}\end{matrix}\right]\)
Resultado: El producto \((x-m_{k})(x-m_{k})^{T}\) genera una matriz cuadrada donde los elementos de la diagonal son los cuadrados de las desviaciones y los elementos fuera de la diagonal son los productos cruzados (que indican la covarianza interna del cluster).

\[ \bf W=\left[\begin{matrix}\LARGE\sum \Delta x_{1}^{2}&\sum \Delta x_{1}\Delta x_{2}&\sum \Delta x_{1}\Delta x_{3}\\ \sum \Delta x_{2}\Delta x_{1} & \LARGE\sum \Delta x_{2}^{2}&\sum \Delta x_{2}\Delta x_{3}\\ \sum \Delta x_{3}\Delta x_{1} & \sum \Delta x_{3}\Delta x_{2} & \LARGE \sum \Delta x_{3}^{2}\end{matrix}\right] \]

La traza \(Tr(\bf W)\) es la suma de la diagonal principal: \(w_{11}+w_{22}+w_{33}\): \[ Tr(\bf W)=\Delta x_{1}^{2}+\sum \Delta x_{2}^{2}+\sum \Delta x_{3}^{2} \]

Matemáticamente la matriz \(\bf W\) se define como la suma de las matrices de dispersión de cada cluster \(k\), es la variación dentro de cada cluster:

\[ W=\sum _{k=1}^{K}\sum _{x\in C_{k}}(x-m_{k})(x-m_{k})^{T} \]

Donde \(x\) es un vector de datos y \(m_{k}\) es el centroide del cluster \(k\).

La variación total será

\[ \mathbf{T}=\sum _{k=1}^{K} \sum _{x\in C_{k}}^{n_{i}} ({x} - m)({x} - m)^T \]

Mientras, la variación entre clusters será:

\[ \mathbf B = \sum_{k=1}^{K} (m_k - m)(m_{k}- m)^T \]

Verificándose pues la igualdad \[ \mathbf{T=W+B} \]

Minimización de la traza de \(\mathbf{W}\)

Hemos visto que, para un hipotético cluster de 3 puntos, la matriz de dispersión \(\bf W\) sería

Y la traza \(Tr(\bf W)\) en esta matriz, que es la suma de la diagonal principal: \(w_{11}+w_{22}+w_{33}\), sería: \[ Tr(\bf W)=\Delta x_{1}^{2}+\sum \Delta x_{2}^{2}+\sum \Delta x_{3}^{2} \]

Por lo tanto, podemos decir que la traza es \(\sum _{j=1}^{k}Tr(W_{j})\). Esto es equivalente a minimizar: \[ Tr(\mathbf W)=\sum _{k=1}^{K}\sum _{x\in C_{k}}\|x-m_{k}\|^{2}=d_{i\in C_k}^2 \]

Donde \(\|x-m_{k}\|^{2}\) es la distancia euclidiana al cuadrado entre el indivíduo y la media del cluster al cual pertenece. La traza de la matriz \(\bf W\) es la métrica de cohesión de los clusters; cuanto menor sea su valor, más compactos son los grupos.

Esta solución es la que persigue k-means.

Otros métodos similares

La minimización del determinante de \(\mathbf{W}\)

Este método está basado en el cociente \(\frac{det(\bf T)}{det(\bf W)}\).

Diferencia con K-means (Traza de W): Mientras que el método K-means clásico minimiza la traza de \(W\) (suma de varianzas, buscando formas esféricas), minimizar el determinante de \(W\) permite que los clusters tengan formas elipsoidales alargadas.

Propiedades: Los algoritmos basados en la minimización del determinante de \(W\) son invariantes no solo a cambios de traslación, sino también a cambios de escala (invarianza afín), lo que los hace útiles cuando las variables tienen diferentes unidades de medida.

La maximización de la traza de \(\mathbf{BW^{-1}}\)

El criterio busca encontrar una partición de los datos que maximice el valor propio de la matriz, es decir:

Maximizar la distancia entre grupos \((\bf B)\) y, al mismo tiempo,
Minimizar la dispersión dentro de los grupos \((\bf W)\).

Al utilizar \(\bf BW^{-1}\), estamos buscando un ratio: maximizar la separación inter-grupos relativa a la dispersión intra-grupo. Esto es conceptualmente análogo al criterio de Fisher en el análisis discriminante lineal (LDA)¹, pero aplicado al aprendizaje no supervisado.

El algoritmo k-means con R

En R, el análisis de cluster con el algoritmo k-means se realiza con la función de stat \[ \texttt{kmeans(m,r)} \]

Donde \(\tt m\): Matriz de datos \(tt r\): Número de clusters que se quieren.

Ejemplos

Ejemplo 1

Volvemos al ejemplo utilizado en otras ocasiones, de EAAR, Cap. 5² sobre estatura en cm y sueldos en Pts de 5 individuos.

\[\begin{matrix} & \begin{matrix} \phantom{i} X_1\phantom{i} & X_2\phantom{i} \end{matrix} \\ \begin{matrix} 1\\ 2\\ 3\\ 4\\ 5\\ \end{matrix} & \begin{pmatrix} 180 & 175 \\ 170 & 180 \\ 175 & 165 \\ 189 & 100 \\ 185 & 110 \\ \end{pmatrix} \\ \end{matrix}\]

Estos dato se pueden representar como puntos en un gráfico de dispersión:

cl1 <- kmeans(x= df,centers = 2)

Composición de los clusters

cl1$cluster

## 1 2 3 4 5 
## 2 2 2 1 1

Centroides de los clusters

cl1$centers

##   X_1      X_2
## 1 187 105.0000
## 2 175 173.3333

Suma de cuadrados dentro de cada cluster \(\bf W\)

cl1$withinss

## [1]  58.0000 166.6667

Suma de los cuadrados \(T\)

cl1$totss

## [1] 6000.8

Suma de los cuadrados \(B\)

cl1$betweenss

## [1] 5776.133

\[\bf T=B+W\]

cl1$tot.withinss + cl1$betweenss

## [1] 6000.8

Representación gráfica

plot(df, col=cl1$cluster,pch=16)

ggplot(data = df, aes(x=X_1,y=X_2, colour = factor(cl1$cluster))) +
  geom_point() +
  geom_point(data = cl1$centers, aes(x=cl1$centers[,1],y= cl1$centers[,2]),
             size=4, color= c("brown","darkblue"))+
  guides(col = guide_legend(title = "Cluster")) +
  theme_classic() +
  labs(title = "k-means: Ejemplo 1",
       subtitle = "Clústers de puntos y centroides")+
  ylab(TeX("$X_2$")) +
  xlab(TeX("$X_1$"))

En esta solución con \(\tt ggplot\) hemos añadido la repesentación de los centroides (o medias), que puede ser interesante para ubicar mejor los clusters.

Ejemplo 2

Otro ejemplo de la librería de datos \(\tt cluster.datasets\) es el datasets nutrients.meat.fish.fowl.1959: Tabla con los niveles de nutrientes en carne, pescado y aves. Los niveles de nutrientes se midieron en una porción de 85 g de diversos alimentos. Esta es la Tabla 4.1 del Capítulo 4 de Hartigan (1975), página 86.

library(cluster.datasets)
data("nutrients.meat.fish.fowl.1959")

nut <- nutrients.meat.fish.fowl.1959

head(nut)

##              name energy protein fat calcium iron
## 1    Braised beef    340      20  28       9  2.6
## 2       Hamburger    245      21  17       9  2.7
## 3      Roast beef    420      15  39       7  2.0
## 4       Beefsteak    375      19  32       9  2.6
## 5     Canned beef    180      22  10      17  3.7
## 6 Broiled chicken    115      20   3       8  1.4

Transformamos en matriz (formato que puede ser procesado tanto por \(\tt keameans\) como por \(\tt hclus\)) manteniendo los nombres como \(\tt rownames\), y escalamos.

rownames(nut) <- nut$name

nut <- nut[,2:6]

nut <- as.matrix(nut)

nut <- scale(nut)

head(nut)

##                     energy    protein        fat    calcium       iron
## Braised beef     1.3101024  0.2352002  1.2897287 -0.4480464  0.1495365
## Hamburger        0.3714397  0.4704005  0.3125618 -0.4480464  0.2179685
## Roast beef       2.1005553 -0.9408009  2.2668955 -0.4736761 -0.2610553
## Beefsteak        1.6559256  0.0000000  1.6450621 -0.4480464  0.1495365
## Canned beef     -0.2708033  0.7056007 -0.3092717 -0.3455273  0.9022882
## Broiled chicken -0.9130462  0.2352002 -0.9311051 -0.4608612 -0.6716471

Creamos un cluster de 4 centroides

cl2 <- kmeans(nut, centers = 4)

cl2

## K-means clustering with 4 clusters of sizes 14, 8, 2, 3
## 
## Cluster means:
##       energy     protein        fat    calcium        iron
## 1 -0.4295996  0.40320040 -0.5123193 -0.3089133 -0.24150330
## 2  1.3286287 -0.05880006  1.3674579 -0.4512501  0.03833458
## 3 -1.4811842 -2.35200232 -1.1087718  0.4361807  2.27092763
## 4 -0.5507554 -0.15680015 -0.5165495  2.3541419 -0.48916187
## 
## Clustering vector:
##        Braised beef           Hamburger          Roast beef           Beefsteak 
##                   2                   1                   2                   2 
##         Canned beef     Broiled chicken      Canned chicken          Beef heart 
##                   1                   1                   1                   1 
##      Roast lamb leg Roast lamb shoulder          Smoked ham          Pork roast 
##                   2                   2                   2                   2 
##       Pork simmered         Beef tongue         Veal cutlet      Baked bluefish 
##                   2                   1                   1                   1 
##           Raw clams        Canned clams     Canned crabmeat       Fried haddock 
##                   3                   3                   1                   1 
##    Broiled mackerel     Canned mackerel         Fried perch       Canned salmon 
##                   1                   4                   1                   4 
##     Canned sardines         Canned tuna       Canned shrimp 
##                   4                   1                   1 
## 
## Within cluster sum of squares by cluster:
## [1] 28.9724028  4.3364978  0.5626097  6.9584784
##  (between_SS / total_SS =  68.6 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"

Representación con \(\tt plot\)

plot(nut[,1],nut[,2], col= cl2$cluster,pch=16,
     xlab = "Energy", ylab = "Protein")

Representación con \(\tt ggplot\)

as.data.frame(nut) %>% 
  ggplot(aes(x=energy,y= protein, )) +
  geom_point(aes(colour= factor(cl2$cluster))) +
  geom_point(data = cl2$centers, aes(x=cl2$centers[,1],y= cl2$centers[,2],
                                    colour = rownames(cl2$centers)), size=5)+
  ggforce::geom_mark_ellipse(aes(group = factor(cl2$cluster), colour = factor(cl2$cluster)))+
  guides(col = guide_legend(title = "Cluster")) +
  theme_classic() +
  #ggrepel::geom_text_repel(aes(label = rownames(nut),colour = factor(cl2$cluster)))+
  labs(title = "k-means: Ejemplo 2",
       subtitle = "Clústers de puntos y centroides")+
  ylab("Protein") +
  xlab("Energy")

Determinación del número de clusters: El método del codo (elbow mwthod)

El método del codo analiza el porcentaje de varianza explicada en función del número de conglomerados: se debe seleccionar un número de conglomerados tal que añadir otro no mejore la modelización de los datos. Más precisamente, si se grafica el porcentaje de varianza explicada por los conglomerados en función del número de conglomerados, los primeros conglomerados aportarán mucha información (explicarán mucha varianza), pero en algún momento la ganancia marginal disminuirá, generando un ángulo en el gráfico. En este punto se selecciona el número de conglomerados, de ahí el “criterio del codo”. Este “codo” no siempre se puede identificar con precisión.

# Calcular la suma total de los cuadrados dentro del cluster 
# Total within-cluster sum of squares (WSS)


# FOR LOOP
n.max <- 10 # Número máximo de clusters a plantear

suma <- c()

for (i in 1:n.max) { # A partir de 2 clusters, claro... 
  suma[i] <- sum(kmeans(nut, centers = i, nstart = 50, iter.max = 15)$withinss)
}

suma

##  [1] 130.000000  85.939021  59.259501  40.829989  28.145711  21.780432
##  [7]  16.408219  12.250982  10.097726   8.127928

Con \(\tt plot\)

plot(1:n.max,suma, type = "b", frame = FALSE,
     ylab = "Suma de cuadrados dentro de los clusters", xlab = "Número de clusters",
     main = "k-means - Gráfico del 'elbow'")

Con \(\tt ggplot\)

data.frame(cbind(1:n.max,suma)) %>% 
  ggplot(aes(x=V1,y=suma)) +
  geom_point()+
  geom_line() +
  scale_x_continuous(breaks = c(1,2,3,4,5,6,7,8,9,10))+
  xlab("Número de clusters")+
  ylab("Suma de cuadrados dentro de los clusters")+
  labs(title = "k-means: Ejemplo 2",
       subtitle = "Gráfico del 'elbow'")

Podría interpretarse un ligero cambio de tendencia que podría apuntar a seleccionar 4 clusters. Sin embargo, yo no lo acabo de ver claro.

Determinación del número de clusters: agrupamiento jerárquico (método Ward)

d <- dist(nut)

nut.hcl <- hclust(d, method = "ward.D")

plot(nut.hcl)
rect.hclust(nut.hcl, k=4)

Directamente visualizando los 4 cluster identificado sobre los ejes de \(PC_1\) y \(PC_2\)

library(cluster)

clusplot(nut,cl2$cluster, 
         color = T, shade = T,
         main = "4 Clusters",labels = 2,col.p = cl2$cluster,cex.txt = 0.7)

Ejemplo 3: PCA previa, para determinar las variables que aportan más variabilidad

Se van a utilizar los mismos datos que en el Ejemplo 2.

Reducción del número de variables: Método Jollife

Autovalores

princomp(nut,cor = T)$sdev

##     Comp.1     Comp.2     Comp.3     Comp.4     Comp.5 
## 1.48199981 1.07046241 0.92071355 0.89915054 0.04002061

Peso autovalores sobre el total

cumsum(princomp(nut,cor = T)$sdev) / sum(princomp(nut)$sdev)

##    Comp.1    Comp.2    Comp.3    Comp.4    Comp.5 
## 0.3422738 0.5895015 0.8021440 1.0098064 1.0190493

Loadings de las PC

princomp(nut, cor = T)$loadings[,]

##             Comp.1       Comp.2     Comp.3    Comp.4       Comp.5
## energy   0.6542825  0.085330063  0.1506898 0.1970344  0.709297704
## protein  0.1507319 -0.689332416 -0.4626844 0.5264918 -0.104068653
## fat      0.6399107  0.199787072  0.2174827 0.1317095 -0.697103521
## calcium -0.3549785 -0.003144364  0.6513825 0.6705754  0.003160074
## iron    -0.1170424  0.691096838 -0.5400166 0.4657952  0.010157612

Para seleccionar las variables con este método Método Jollife

Seleccionar un subconjunto entre las variables p de esta manera:

Elegir las q Componentes Principales cuyos autovalores sean mayores que 0.7: \(PC_1+PC_2+PC_3=0.80\)
Dentro de dichas \(CP\) elegidas, seleccionar la variable con carga mayor que no haya sido aún seleccionada: \[ \begin{array} {lll} PC_1 & Energy & 0.65 \\ PC_2 & Iron & 0.69 \\ PC3 & Calcium & 0.65 \\ \end{array} \]

Representación de k-means utilizando las variables identificadas con el \(PCA\)

nutPC <- nut[,c(1,4:5)]

n.max <- 10 # Número máximo de clusters a plantear

suma <- c()

for (i in 1:n.max) { # A partir de 2 clusters, claro... 
  suma[i] <- sum(kmeans(nutPC, centers = i, nstart = 50, iter.max = 15)$withinss)
}

plot(1:n.max,suma, type = "b", frame = FALSE,
     ylab = "Suma de cuadrados dentro de los clusters", xlab = "Número de clusters",
     main = "k-means - Gráfico del 'elbow' - Ejemplo 3")

Seleccionando las variables con el Método PC / Jollylife, el codo queda mucho más evidente en \(k=4\).

Realizamos la agregación con \(\tt kmeans\), utilizando las 3 variables seleccionadas: \(\tt energy,iron,calcium\)L

cl2PC <- kmeans(nutPC,centers = 4, nstart = 50, iter.max = 15)
cl2PC

## K-means clustering with 4 clusters of sizes 1, 13, 4, 9
## 
## Cluster means:
##       energy     calcium        iron
## 1 -0.2708033  4.13968254  0.08110456
## 2 -0.5406213 -0.01726906 -0.65059113
## 3 -0.9253971  0.03571556  1.96298376
## 4  1.2222743 -0.45089411  0.05829390
## 
## Clustering vector:
##        Braised beef           Hamburger          Roast beef           Beefsteak 
##                   4                   4                   4                   4 
##         Canned beef     Broiled chicken      Canned chicken          Beef heart 
##                   3                   2                   2                   3 
##      Roast lamb leg Roast lamb shoulder          Smoked ham          Pork roast 
##                   4                   4                   4                   4 
##       Pork simmered         Beef tongue         Veal cutlet      Baked bluefish 
##                   4                   2                   2                   2 
##           Raw clams        Canned clams     Canned crabmeat       Fried haddock 
##                   3                   3                   2                   2 
##    Broiled mackerel     Canned mackerel         Fried perch       Canned salmon 
##                   2                   2                   2                   2 
##     Canned sardines         Canned tuna       Canned shrimp 
##                   1                   2                   2 
## 
## Within cluster sum of squares by cluster:
## [1]  0.000000 11.344966  3.529996  2.452567
##  (between_SS / total_SS =  77.8 %)
## 
## Available components:
## 
## [1] "cluster"      "centers"      "totss"        "withinss"     "tot.withinss"
## [6] "betweenss"    "size"         "iter"         "ifault"

as.data.frame(nutPC) %>% 
  ggplot(aes(x=energy, y= iron)) +
  geom_point(aes(colour= factor(cl2PC$cluster))) +
  geom_point(data = cl2PC$centers, aes(x=cl2PC$centers[,1],y= cl2PC$centers[,3],
                                    colour = rownames(cl2PC$centers)), size=4)+
  ggforce::geom_mark_ellipse(aes(group = factor(cl2PC$cluster), colour = factor(cl2PC$cluster)))+
  guides(col = guide_legend(title = "Cluster")) +
  theme_classic() +
  #ggrepel::geom_text_repel(aes(label = rownames(nut),colour = factor(cl2PC$cluster)))+
  labs(title = "k-means: Ejemplo 3",
       subtitle = "Clústers de puntos y centroides: variables seleccionadas con PCA (Método Jollylife)")+
  ylab("Iron") +
  xlab("Energy")

La selección prévia permite una visualización mejor en un gáfico de disperión, poniendo en los ejes las dos primeras variables identificadas con el método Lollylife.

Agrupamiento jerárquico (método completo) usando las variables seleccionadas con la \(PCA\)

d <- dist(nutPC)

nut.hcl <- hclust(d, method = "complete")

plot(nut.hcl)
rect.hclust(nut.hcl, k=4)

Esta metodología es la que tiene el resultado más ajustado respecto al del algoritmo k-means: la única diferencia evidente es la del ítem Canned beaf que cambia de cluster. Sin embargo, está decentrado respecto a su agrupación también en k-means.

Ejemplo 3 con \(BIC\) (Bayesian Information Criterion) y la libería \(\tt mclust\)

library(mclust)

mPCmcl <- Mclust(nutPC)

summary(mPCmcl, parameters = T)

## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust VVI (diagonal, varying volume and shape) model with 3 components: 
## 
##  log-likelihood  n df       BIC       ICL
##       -29.44168 27 20 -124.8001 -125.4608
## 
## Clustering table:
##  1  2  3 
## 11  8  8 
## 
## Mixing probabilities:
##         1         2         3 
## 0.4072285 0.3073134 0.2854580 
## 
## Means:
##                [,1]       [,2]       [,3]
## energy   0.97822401 -0.7983098 -0.5360837
## calcium -0.45270162  0.9466365 -0.3732986
## iron     0.07485667  0.7505871 -0.9148430
## 
## Variances:
## [,,1]
##           energy      calcium       iron
## energy  0.478114 0.000000e+00 0.00000000
## calcium 0.000000 9.764191e-05 0.00000000
## iron    0.000000 0.000000e+00 0.01741992
## [,,2]
##            energy  calcium     iron
## energy  0.2142837 0.000000 0.000000
## calcium 0.0000000 1.821872 0.000000
## iron    0.0000000 0.000000 1.708232
## [,,3]
##            energy   calcium       iron
## energy  0.1220356 0.0000000 0.00000000
## calcium 0.0000000 0.0154543 0.00000000
## iron    0.0000000 0.0000000 0.05807942

plot(mPCmcl, what = "classification")

library(factoextra)

fviz_mclust(mPCmcl, what = "classification", labelsize=7)

Gracias a haber realizado una PCA previas el algoritmo puede identificar 3 clusters (en pruebas que hemos hecho, sin PCA, serían 8 clusters pequeños y muy confusos) que se distribuyen de manera diáfana tanto en los planos de las variables como en el de las dos primeras Componentes Principales.

Conclusión sobre ejemplos 2 y 3

Realizar una \(PCA\) previa al proceso de clustering, puede mejorarlo bastante, aumentando su claridad y disminuyendo el ruido provocado por componentes secundarias. Con el método Jollylife. Podemos seleccionar de manera eficaz las variables que más aportan a la variabilidad del conjunto.

LDA: Ver capítulo sobre PCA, Análisis de Proyección Óptima, donde se utiliza el LDA.↩︎
A. García Perez, Estadística Aplicada Avanzada con R, Madrid, 2022↩︎

Cluster Analysis 2 - Algoritmo k-means

Stefano Sanfilippo

2026-02-12