Qué es y cómo programar Kmeans desde 0 en R

Share on linkedin
Share on twitter
Share on email

El algoritmo K-means es uno de los algoritmos de machine learning no supervisados más utilizados. Sin embargo, ¿conoces exactamente cómo funciona y cómo y cuándo deberías usarlo? En este post vamos a programar un Kmeans desde 0 en R y vamos a descubrir todas sus particulares. ¡Vamos a ello!

La base del algoritmo k-means: optimizando la suma errores al cuadrado

La idea del algoritmo es agrupar cada observación (fila de una tabla) en k grupos, llamados cluserts, que decida el usuario. De ahí su nombre k-means. El número de grupos estará, por tanto, definido por el usuario, por lo que puede ser un poco subjetivo, aunque luego veremos qué metodos podemos usar para calcular el número de clusters óptimo. Pero, en cualquier caso, ¿cómo crea el algoritmo de Kmeans estos clusters?

Para ello, el algoritmo sigue los sigueentes pasos:

  1. Se inicializan k puntos aleatorios, llamados centroides.
  2. Para cada observación se calcula la suma de errores al cuadrado de esa observación respecto a cada uno de los k centroides.
  3. A cada observación se le asigna el centroide que menos error tenga.

Dicho así, hablar de errores suena un poco abstracto, así que veámos cómo funciona gráficamente (por cierto, la función circulo la usé también en el tutorial de programar una red neuronal desde 0 en R, por si quieres ver qué hace;).

Lo primero de todo, antes de meternos a programar nuestro Kmeans desde 0 en R, vamos a crear nuestros datos:

library(ggplot2)
set.seed(123)

# Creamos puntos dispersos sobre un punto
datos_1 <- circulo(x = 20,R = 10, centroX = 5, centroY = 30)
datos_2 <- circulo(x = 20,R = 10, centroX = 20, centroY = 10)
datos_3 <- circulo(x = 20,R = 10, centroX = 50, centroY = 50)


datos = data.frame(
  rbind(datos_1,datos_2, datos_3),
  numero_punto = 1:60
)

rm(datos_1,datos_2,datos_3, circulo)

Ahora vamos a crear una función que cree centroides y crearemos nuestros centroides. Además, dibujaremos de forma gráfica el error al cuadrado de cada de la primera observación sobre cada uno de los centroides.

crear_centroides <- function(n_centroides = 3,datos = datos, columns_skip){
  datos_limpios = datos[ , -which(names(datos) %in% columns_skip)]
  columnas <- c(letters[24:26],letters[1:23])

  x = matrix(ncol = ncol(datos_limpios) +1, nrow = n_centroides)
  x[1:n_centroides,1:ncol(datos_limpios)] = replicate(ncol(datos_limpios), runif(n_centroides,min(datos_limpios), max(datos_limpios))) 
  x[,ncol(datos_limpios) +1] = 1:n_centroides
  x = data.frame(x, stringsAsFactors = FALSE)
  colnames(x) <- c(columnas[1:ncol(x)-1], "n_centroide")
  return(x)
}

centroides = crear_centroides(n_centroides = 3, datos = datos, columns_skip = "numero_punto")

ggplot() +
  geom_text(aes(x,y, label = numero_punto), data = datos) + 
  geom_point(aes(x,y), col = "blue",alpha = 0.2, data = datos) + 
  geom_point(aes(x,y), data = centroides, size = 10, shape=21) +
  geom_text(aes(x,y, label = n_centroide), data = centroides, col = "black") +

  geom_segment(aes(x = datos[1,"x"], y =  datos[1,"y"],
                  xend =   centroides[1,"x"], yend =  centroides[1,"y"]), linetype = "dashed") +

  geom_segment(aes(x = datos[1,"x"], y =  datos[1,"y"],
                  xend =   centroides[2,"x"], yend =  centroides[2,"y"]), linetype = "dashed") +
  geom_segment(aes(x = datos[1,"x"], y =  datos[1,"y"],
                  xend =   centroides[3,"x"], yend =  centroides[3,"y"]), linetype = "dashed") +

  theme_minimal() +
  labs(title = "Error del punto 1 respecto a los 3 centroides")
Iniciación aleatoria de los centroides al programar el algoritmo kmeans desde 0 en R

Como podemos ver, la suma de errores al cuadrado es igual a la distancia euclídea entre las observaciones. Esto se dará independientemente del número de variables que tengamos. Vamos, que lo que Kmeans indirectamente hace es encontrar el centroide que está linealmente más cerca de la observación.

Por si no te acuerdas, la distancia euclídea se basa en el teorema de Pitágoras y básicamente se traduce en la siguiente fórmula (en el caso de dos dimensiones):

Fórmula de la distancia euclidea

Kmeans optimiza el error al cuadrado, que justo coincide con la distancia Euclídea sin aplicar la raíz. Por eso mismo muchas veces se confunde y se piensa que Kmeans mide distancias… aunque no es del todo cierto.

En nuestro caso lo haremos todo desde cero, aunque podrías usar la función distance de la librería philentropy o incluso usar la función kmeans de R. Pero bueno, programarlo nosotros mismos nos permitirá ver:

  1. Qué problemáticas tiene este algoritmo y qué hay que hacer para aplicarlo correctamente.
  2. Entender a la perfección cómo funciona.

Ambos puntos, en mi opinión son muy interesantes. Así que, veamos cómo programar el algoritmo Kmeans desde cero en R.

Como hemos dicho, lo primero que hay que hacer para programar el algoritmo Kmeans desde 0 (en R o donde quieras) es calcular el error al cuadrado de las observaciones respecto a los centroides. Veámoslo cómo se hace con la primera observación.

# Distancia 1
(centroides[1,1] - datos[1,1])^2 + (centroides[1,2] - datos[1,2])^2 
(centroides[2,1] - datos[2,1])^2 + (centroides[2,2] - datos[2,2])^2  
(centroides[3,1] - datos[3,1])^2 + (centroides[3,2] - datos[3,2])^2 
## [1] 1036.13
## [2] 199.8957
## [3] 1251.229
prueba <- c(
  (centroides[1,1] - datos[1,1])^2 + (centroides[1,2] - datos[1,2])^2,
  (centroides[2,1] - datos[2,1])^2 + (centroides[2,2] - datos[2,2])^2,
  (centroides[3,1] - datos[3,1])^2 + (centroides[3,2] - datos[3,2])^2  
)

which(prueba == min(prueba))
## [1] 2

Como podemos ver, el clúster 2 es el linealmente más cercano a la observación 1. Ahora tendríamos que repetir ese proceso para cada uno de los puntos, como a continuación. Además del clúster, guardaremos también el error, que nos servirá para más adelante.

datos$error = NA
datos$cluster = NA

for(posicion in 1:length(datos$x)){
  x <- c(
    (centroides[1,1] - datos[posicion,1])^2 + (centroides[1,2] - datos[posicion,2])^2,  
    (centroides[2,1] - datos[posicion,1])^2 + (centroides[2,2] - datos[posicion,2])^2, 
    (centroides[3,1] - datos[posicion,1])^2 + (centroides[3,2] - datos[posicion,2])^2  
  )

  datos$error[posicion] = min(x)
  datos$cluster[posicion] = which(x == min(x))

}
error = sum(datos$error)

rm(x, posicion)

Vamos a ver en qué grupos ha clasificado nuestro Kmeans a cada individuo:

ggplot() + 
  geom_point(aes(x,y, col = as.factor(cluster), size = 7), data = datos) + 
    geom_point(aes(x,y), data = centroides, size = 10, shape=21) +
  geom_text(aes(x,y, label = n_centroide), data = centroides, col = "black") +
  theme_minimal() + theme(legend.position = "bottom") + guides(size = FALSE) 
Resultado del algoritmo kmeans en el primer entrenamiento.

Como vemos, el resultado dista bastante de los grupos que visualmente vemos. Pero tranquilo, que eso no es todo, todavía queda ir puliendo los clusters. ¡Veámos!

Determinando el clúster al que perteneces

Ahora que ya tenemos un primer acercamiento de a qué cluster pertenece cada observación, hay que hacer que el algoritmo de Kmeans aprenda y así ir mejorando los clusters. Para ello tenemos que recalcular los centroides, que se representan calcula como el promedio de para cada uno de los ejes.

library(dplyr)
centroides <- datos %>%
  group_by(cluster) %>%
  summarize(
    x = mean(x), 
    y = mean(y)
    ) %>%
  mutate(n_centroide = cluster) %>%
  select(-cluster) %>%
  ungroup() %>%
  as.data.frame(.)


ggplot() + 
  geom_point(aes(x,y, col = as.factor(cluster), size = 7), data = datos) + 
    geom_point(aes(x,y), data = centroides, size = 10, shape=21) +
  geom_text(aes(x,y, label = n_centroide), data = centroides, col = "black") +
  theme_minimal() + theme(legend.position = "bottom") + guides(size = FALSE) 
Reposición de los clusters del algoritmo kmeans en el primer entrenamiento

Vemos cómo se han movido los centroides. Ahora tenemos que repetir el proceso, hasta que el error total no se reduzca más. Por eso mismo hemos calculado el error respecto a los centroides, para ir comprobando si este se reduce o no.

Este proceso lo podemos hacer de forma iterativa mediante un bucle while.

# Creamos un error previo para poder comparar en la primera iteración.
error = c(0,error)

i = 2
while(round(error[i],2) != round(error[i-1],2) ){

  # Calculamos la distancia
  for(posicion in 1:length(datos$x)){
  x <- unlist( c(
    (centroides[1,1] - datos[posicion,1])^2 + (centroides[1,2] - datos[posicion,2])^2,  
    (centroides[2,1] - datos[posicion,1])^2 + (centroides[2,2] - datos[posicion,2])^2, 
    (centroides[3,1] - datos[posicion,1])^2 + (centroides[3,2] - datos[posicion,2])^2  
  ))

  # Actualizamos campos
  datos$error[posicion] = min(x)
  datos$cluster[posicion] = which(x == min(x))
}

  # Actualizamos campos una vez terminado el cálculo de distancias
  error = c(error, sum(datos$error))

centroides <- datos %>%
  group_by(cluster) %>%
  summarize(
    x = mean(x), 
    y = mean(y)
    ) %>%
  mutate(n_centroide = cluster) %>%
  select(-cluster) %>%
  ungroup() %>%
  as.data.frame(.)


  i = i + 1

}
rm(i,x, posicion, error)

Con esto nuestro Kmeans ya habría clasificado a nuestras observaciones en distintos grupos. Veamos a ver cuáles son esos grupos:

ggplot() + 
  geom_point(aes(x,y, col = as.factor(cluster), size = 7), data = datos) + 
    geom_point(aes(x,y), data = centroides, size = 10, shape=21) +
  geom_text(aes(x,y, label = n_centroide), data = centroides, col = "black") +
  theme_minimal() + theme(legend.position = "bottom") + guides(size = FALSE) 
Resultado del algoritmo kmeans en el priemr aprendizaje

Vaya, parece que visualmente el algoritmo no ha conseguido clasificar los datos de forma correcta… Esto de programar nuestro Kmeans desde 0 en R no es tan bonito… ¿o sí? Y es que nos acabamos de encontrar con uno de los problemas que puede encontrarse el algoritmo Kmeans: no clasificar bien.

Por suerte, este problema tiene solución. Pero antes de ello, hay que entender bien las causas, que son dos:

  • Los datos no tienen la misma escala: imaginemos que queremos agrupar a personas en función de ingresos y años trabajados. Claro, los ingresos están en una escala mucho mayor (valores más altos), por lo que influirán más al calcular la distancia. Para evitar eso, siempre se recomienda estandarizar los datos antes de aplicar un clustering más info aquí.
  • Aleatoriedad de la posición inicial: los centroides se inicializan con valores aleatorios, por lo que podemos tener la mala suerte de que nuestros centroides se hayan inicializado con unos valores que perjudican el clustering, como en nuestro caso.

Como os decía, programar Kmeans (o cualquier otro algoritmo) desde 0 (en R) tiene sus ventajas: te encuentras los problemas del algoritmo y aprendes a utilizarlo mucho mejor.

En cualquier caso, ¿cómo resolvemos estos dos problemas? ¡Veamos!

Resolviendo posibles problemas de Kmeans

Estandarización de los datos

Estandarizar los datos significa hacer que estos tengan promedio cero y desviación típica de uno. Para conseguirlo, a cada valor simplemente tenemos que:

  1. Restarle la media.
  2. Dividirlo entre su desviación típica.

En nuestro caso, este hecho seguramente no tenga un gran impacto, ya que los datos se han generado aleatoriamente de una misma distribución… Pero bueno, hacerlo cuesta poco y siempre es recomendable.

# Estandarización de X
datos$x <- datos$x - mean(datos$x) 
datos$x <- datos$x / sd(datos$x)

# Estandarización de y
datos$y <- datos$y - mean(datos$y) 
datos$y <- datos$y / sd(datos$y)

Ahora nuestros datos tienen la siguiente forma:

ggplot() + 
  geom_point(aes(x,y, col = as.factor(cluster), size = 7), data = datos) + 
    theme_minimal() + theme(legend.position = "bottom") + guides(size = FALSE) 
Reusltado del algortimo kmeans tras estandarizar. Programar el algoritmo kmeans desde 0 en R

Ahora que ya hemos estandarizado los datos (aunque no haya generado un gran cambio), vamos a ver cómo arreglar la aleatoriedad de los clústers.

Resolviendo la aleatoriedad de los centroides iniciales

Resolver el problema de la aleatoriedad inicial es una cuestión de estadística: si solo haces Kmeans una vez, puede que te pase, pero si lo haces muchas… es poco probable que pase siempre. Además, contamos con la suerte de que el algoritmo Kmeans es bastante rápido, por lo que no hay problema en correrlo varias veces.

Por tanto, la forma de resolver el problema de aleatoriedad de los centroides iniciales es realizar varios Kmeans.

Digamos que vamos a realizar 5 rondas. Seguramente en alguna de ellas los centroides se inicialicen de una forma más “normal” para nuestros datos. En ese caso, la suma total de los errores al cuadrado será menor que en el resto de Kmeans. Así, tendremos nuestro ganador.

Para hacerlo vamos a seguir los siguientes pasos:

  1. Creamos los centroides. Para ello usaremos la función crear_centroides que hemos creado anteriormente.
  2. Nos quedamos con los datos que nos interesan (restando las columnas columns_skip) y convertimos los datos en una matriz para poder trabajar matricialmente y que el proceso sea más eficiente.
  3. Para cada una de las filas, calculamos el error al cuadrado hacia cada uno de los centroides. Para ello, tendremos que convertir la observación en una matriz que tenga tantas filas como clusters hay.
  4. Encontramos la posición de los centroides que minimiza el error y asignamos el número de clúster y el error a esa observación. Repetimos este proceso para todas las observaciones.
  5. Al terminar de calcular errores, calculamos el error total de las observaciones sobre los centroides. Si es el error es igual al error anterior, detenemos el proceso. Sino, continuamos.
  6. Si el error no es el mismo, recalculamos la posición de los centroides.
  7. Cuando el error ya no se reduce más significa que los clusters están en su posición óptima. Por tanto, devolvemos la posición de los centroides, los errores y número de clúster de cada una de las observaciones.

Para poder realizar este proces, definiremos una función llamada elevar que devuelva el cuadrado de un número. Esto será necesario para elevar al cuadrado los resultados de las diferencias, ya que elevar los elementos de una matriz parece que no es algo se pueda hacer por defecto en R.

kmeans_personalizada <- function(numero_iteraciones, datos, n_centroides, columns_skip){
  # Inicializamos una lista donde guardaremos los errores y clusters 
resultados = list()

for(iteracion in 1:numero_iteraciones){

  set.seed(iteracion)
  # 1. Creamos lo centroides
  centroides = crear_centroides(n_centroides = n_centroides, datos = datos, 
                                columns_skip = columns_skip)

  centroides <- as.matrix(centroides)
  # centroides <- as.matrix(centroides[,1:ncol(centroides)-1])

  # 2. Converitmos los datos en matriz
  datos_limpios = datos[ , -which(names(datos) %in% columns_skip)]
  datos_limpios = as.matrix(datos_limpios)

  error = c(-10,0)

  i = 2
  centroides_que_no_estan = c(1)
  while(round(error[i],2) != round(error[i-1],2) | length(centroides_que_no_estan)> 0 ){

    for(posicion in 1:length(datos$x)){

    #3. calculamos el error al cuadrado hacia cada uno de los centroides
    repeticiones = rep(datos_limpios[posicion,], times = n_centroides)
    centroide_posicion = matrix(repeticiones, nrow = n_centroides, byrow = TRUE)

    errores = centroides[,1:ncol(centroides)-1]  - centroide_posicion
    errores = sapply(errores, elevar)
    errores = matrix(errores, nrow = n_centroides)
    errores = rowSums(errores)

    #4. Encontramos el centroide que minimiza el error
    datos$error[posicion] = min(errores)
    datos$cluster[posicion] = which(errores == min(errores))

  }

    #5. Calculamos el error total de las observaciones sobre los centroides
  error = c(error, sum(datos$error))

  #6. Recalculamos la posición de los centroides
  centroides <- datos %>%
    group_by(cluster) %>%
    summarize(
      x = mean(x), 
      y = mean(y)
      ) %>%
    mutate(n_centroide = cluster) %>%
    select(-cluster) %>%
    ungroup() %>%
    as.matrix(.)

  centroides_que_no_estan = setdiff(seq(1:n_centroides), unique(datos$cluster))

  # Comprobamos que no se ha perdido ningún cluster. Sino, lo incluímos
  if(length(centroides_que_no_estan)> 0){
    centroides_que_no_estan = setdiff(seq(1:n_centroides), unique(datos$cluster))

    centroides_nuevos = crear_centroides(n_centroides = length(centroides_que_no_estan), datos = datos, 
                                columns_skip = columns_skip)
    centroides_nuevos <- as.matrix(centroides_nuevos[,1:ncol(centroides_nuevos)-1])

    centroides_nuevos <- cbind(centroides_nuevos, centroides_que_no_estan)
    colnames(centroides_nuevos)[colnames(centroides_nuevos) == "centroides_que_no_estan"] <- "n_centroide"
    centroides <- rbind(centroides, centroides_nuevos)

  }

  i = i + 1

  }

  #7.Devolvemos la posición de los centroides, los errores y número de cluster de cada una de las observaciones
  lista <- list(datos[,c("error","cluster")], centroides)
  resultados[[iteracion]] = lista 
  iteracion = iteracion + 1

}

  return(resultados)
}

resultados <- kmeans_personalizada(numero_iteraciones = 5,
                     n_centroides = 3, datos = datos, 
                     columns_skip = c("numero_punto","error", "cluster"))

Ahora que ya tenemos los cinco Kmeans creados, vamos a encontrar aquel Kmeans que tenga la menor suma de errores al cuadrado.

iteraciones = c()

for(i in 1:length(resultados)){
  x = sum(resultados[[i]][[1]]$error)
  iteraciones = c(iteraciones,x)
}

iteraciones
## [1] 8.251086 8.251086 8.251086 8.251086 8.251086

Como vemos en este caso todas los nuevos Kmeans que hemos hecho han dado el mismo error, es decir, han llegado al mismo resultado. Veamos a ver cuál es ese resultado:

centroides <- as.data.frame(resultados[[5]][[2]])
datos[,c("error","cluster")] = resultados[[5]][[1]]

ggplot() + 
  geom_point(aes(x,y, col = as.factor(cluster), size = 7), data = datos) + 
    geom_point(aes(x,y), data = centroides, size = 10, shape=21) +
  geom_text(aes(x,y, label = n_centroide), data = centroides, col = "black") +
  theme_minimal() + theme(legend.position = "bottom") + guides(size = FALSE) 
Resultado del aprendizaje del algortimo kmeans. Programar el algoritmo kmeans desde 0 en R

¡Perfecto! Nuestro algoritmo de Kmeans ha aprendido y ha conseguido clusterizar correctamente a todos las observaciones en su clúster correspondiente.

Por último solo quedaría definir una cosa: conocer el número de clusters a elegir.

Kmeans: cómo elegir el número de clusters ideal

En este ejemplo el número de clusters a elegir no ha sido al azar, yo mismo he creado tres grupos diferenciados. Sin embargo, en la vida real no sabemos cuántos clusters hay, de ahí la relevancia del aprendizaje no supervisado.

Para saber el número de clusters vamos a realizar un cluster con distintos valores de k. En este caso, haremos 30 clusters. De cada uno de estos clusters guardaremos el error del mismo y lo dibujaremos en una gráfica.

En este gráfica el error irá disminuyendo: cuantos más clusters haya, menor será el error dentro del cluster y, por tanto, menor va a ser el error total. De hecho, cuando haya tantos clusters como individuos el error será 0, ya que el error de una persona sobre si misma es cero.

Sin embargo, habrá un punto en el que el error disminuya, pero no con tanta pendiente: ese será el número de clusters que debemos elegir. Esta forma de elegir el número de clusters se denomina el método del codo, ya que buscamos el “codo” de la gráfica.

errores = data.frame(
  numero_clusters = seq(1:(nrow(datos)/2)),
  error = NA
)

for(i in 1:length(errores$numero_clusters)){
  cluster <- kmeans_personalizada(numero_iteraciones = 1,
                     n_centroides = i, datos = datos, 
                     columns_skip = c("numero_punto","error", "cluster"))

  errores$error[i] = sum(cluster[[1]][[1]]$error)

}

ggplot(errores, aes(numero_clusters, error)) + geom_line() + geom_point() + 
  theme_minimal() + labs(title = "Error Total según el número de clusters") + scale_x_continuous(breaks = seq(1,30))
Cómo definir el número de clusters mediante la regla del codo con el algoritmo k-means desde 0 en R

Como vemos el error total de los clusters disminuye a medida que aumentamos el número de clusters. Sin embargo, así como pasar de 1 a 2 clusters disminuye mucho el error y pasar de 2 a 3 clústers también, vemos como pasar de 3 a 4 clusters no disminuye mucho el error. De esta manera, podemos elegir el número que debemos fijar en k para que esos clústers representen a la realidad.

¡Ya has aprendido todo lo que hay que saber sobre el algoritmo Kmeans y has aprendido a programar tu propio algoritmo Kmeans desde 0 en R! Hemos visto desde el inicio cómo funciona, hemos visto cómo elegir el número de clusters, los problemas que puede haber… Sin embargo hay más… Y es que, saber cómo funciona el algoritmo está bien, pero… ¿cuándo es recomendable usarlo? ¿Qué otras alternativas hay? ¡Te lo explico!

Ventajas, limitaciones y usos del algoritmo Kmeans

Ventajas del algoritmo Kmeans

Entre las principales ventajas del algoritmo Kmeans son que es sencillo de entender, es útil para diversos tipos de clústers, asegura la convergencia… pero en mi opinión, el aspecto clave del algoritmo Kmeans es que es muy fácil de implementar en datasets grandes, con muchas observaciones y variables.

Y es que la alternativa principal (ojo, no la única) a los métodos de clústering por partición, como es Kmeans, es el clústering jerárquico. Sin embargo, los algoritmos de clustering jerárquicos no son viables en datasets grandes, ya que requieren muchísimos cálculos.

Por tanto, si quieres aplicar un cluster sobre un dataset grandes, te recomiendo usar el algoritmo Kmeans. O, al menos, usar Kmeans eligiendo un gran número de clusters y a partir de ahí combinarlo con el clustering jerárquico.

Desventajas del algoritmo Kmeans

Una de las principales desventajas del algoritmo Kmeans es la aleatoriedad. Ya hemos visto que debido a que los centroides se inicializan de forma aleatoria, podemos terminar teniendo unos clusters que no representan correctamente a la realidad. Es cierto que podemos “arreglar” este problema haciendo varias pruebas… pero esto no siempre es viable, ya que cuanto mayor sea el número de clusters, más pruebas habrá que hacer.

Además existen otros aspectos que suponen ciertas pegas para el algoritmo Kmeans como son la subjetividad a la hora de elegir el número de clusters o el efecto de los outliers en los clusters.

En definitiva, en este post hemos aprendido a cómo programar el algoritmo Kmeans desde 0, las dificultades que este algoritmo presentar y cuándo es mejor que lo apliquemos. Y todo, derivado de un lector que me recomendó escribir al respecto.

Así que si os interesa que escriba sobre algo en concreto os animo a responder la encuesta que os sale en la página o escribirme por linkedin. ¡Nos vemos en la siguiente!