Ejemplo de uso de un Mapa Auto-Organizado (SOM) de Kohonen en R

Los Mapas Auto-Organizados (con las siglas SOM en inglés) son un tipo de red neuronal entrenada como aprendizaje no supervisado, de forma que se modifican repetidamente los pesos de dicha red en respuesta a patrones de activación hasta que una configuración final queda desarrollada.

El modelo, creado por Teuvo Kohonen en 1982, consiste en establecer una correspondencia entre los datos de entrada y un espacio bidimensional, creando mapas topológicos, de manera que datos similares entre si activen neuronas en zonas próximas. Con ello, se pueden clasificar datos para los que no se conoce a priori ningún tipo de ordenación.

Estructura de dos capas en SOM – Imagen vía http://www.cs.us.es/~fsancho/?e=76

Los algoritmos SOM tienen una arquitectura en 2 capas: por un lado, la capa de entrada, formada por N neuronas (una neurona por cada dato de entrada), que se encarga de recibir y transmitir a la capa de salida la información procedente del exterior. Por otro lado, la capa de salida, formada por M neuronas, que es la encargada de procesar la información, crear patrones e identificar las posibles categorías.

Durante la fase de entrenamiento los pesos de toda la vecindad son movidos en la misma dirección, por lo que elementos similares tienden a excitar neuronas adyacentes. Por tanto, se forma un mapa semántico donde muestras similares son mapeadas juntas, y las diferentes aparte.

Los Mapas de Kohonen nos servirán para reducción de dimensionalidad, visualización y clustering.

1. Datos de inicio

Disponemos de un dataset con información acerca de observaciones astronómicas realizadas sobre estrellas. Los datos son de acceso público y se pueden descargar del conocido repositorio de la UCI con el nombre de HTRU2 Data Set. HTRU2 es un conjunto de datos ya clasificados, que indican qué observaciones son estrellas del tipo púlsar y cuales no. Este tipo de estrellas son muy interesantes ya que emiten radiación muy intensa a intervalos cortos y regulares, y su correcta clasificación es importante para los astrónomos.

En nuestro conjunto de datos vamos a tener 17897 observaciones con 8 variables pertenecientes a cada estrella (features), además de una señal clasificadora binaria que descartaremos al tratarse de aprendizaje no supervisado.

Muestra de los primeros datos de nuestro dataframe estelar.

Las 8 variables son: media de la señal recibida, desviación estándar de la señal recibida, curtosis de la señal recibida, sesgo de la señal recibida, media de la curva DM-SNR, desviación estándar de la curva DM-SNR, curtosis de la curva DM-SNR y sesgo de la curva DM-SNR. La última columna, X0, es el clasificador.

2. Carga de librerías y datos

library(kohonen)
library(dplyr)
library(plot3D)
library(plot3Drgl)

La librería kohonen la usaremos para realizar los mapas auto-organizados. Con dplyr manipularemos el dataframe y con las dos librerías de plot3D realizaremos algunas interesantes representaciones gráficas en 3D.

df <- read.csv('HTRU_2.csv')
class <- df$X0 # Guardo el clasificador
n_originales <- colnames(df) # Guardamos nombres de columnas 
df <- select(df, -X0) # Eliminamos el clasificador
n_nuevos <- c('V1', 'V2', 'V3', 'V4', 'V5', 'V6', 'V7', 'V8') # Asignamos nuevos
colnames(df) <- n_nuevos # Cambiamos nombres
rm(n_nuevos)

set.seed(100)
df <- na.omit(df) 
Z <- scale(df,center=T,scale=T) # Estandarización

En primer lugar, cargamos los datos que están en formato csv. Después, extraigo la columna clasificadora que no nos interesa al ser aprendizaje no supervisado. Aprovechamos también para cambiar los nombres de todas las columnas simplemente por comodidad, ya que los nombres originales son largos y es difícil citarlos y/o representarlos. Los nombres originales los dejamos guardados en la variable n_originales para recuperarlos al final.

Establecemos una semilla con set.seed() para garantizar la repetibilidad de resultados. Con na.omit() nos aseguramos que no haya valores nulos y por último y muy importante, normalizamos nuestro conjunto de datos con el comando scale(), ya que usaremos distancias euclidianas para comparar datos.

Podemos echar un último vistazo a los datos y comprobar que no haya anomalías con summary():

Resumen de los datos antes de normalizar.

3. Análisis básico SOM

Comenzamos directamente creando un mapa de Kohonen de 7×7, estructura hexagonal y 1000 iteraciones:

# SOM
carte <- som(Z, grid = somgrid(7,7,"hexagonal"), rlen = 1000)
plot(carte, shape = "straight")

Estamos creando el mapa en la variable carte. Con el comando plot() realizamos la primera representación básica que nos ofrece la librería kohonen:

Mapa de Kohonen 7×7 y estructura hexagonal.

Podemos observar como en la parte superior izquierda y superior central del mapa las variables V1 y V2 son las que mayor peso tienen. En la parte inferior derecha predominan las señales V3 y V4 y en la parte inferior izquierda mandan V7 y V8. V5 y V6 tienen mayor peso en la parte superior derecha.

Podemos comprobar si la solución converge con las 1000 iteraciones indicadas mediante una gráfica de training progress:

# Training progress for SOM
plot(carte, type="changes")

Procedemos ahora a representar la gráfica de densidad. El número de instancias en cada celda nos servirá para identificar áreas de alta densidad. La librería kohonen de R nos permitirá, además, especificar la escala de colores a utilizar:

# Mapa de densidad con tonos azules
degrade.bleu <- function(n){
   return(rgb(0,0.4,1,alpha=seq(0,1,1/n)))}

plot(carte, type="count", shape = "straight", palette.name = degrade.bleu)

La función degrade.bleu() toma el parámetro “n” como entrada, que es número de colores a producir. Se genera un conjunto de tonos azules más o menos opacos usando rgb(). Cuanto mayor es el número “n”, más oscuro es el color. Posteriormente, diujamos el mapa con el comando plot() especificando las opciones correctas. El resultado es el siguiente:

Mapa de densidad.

Podemos ver el número de elementos en cada nodo de la siguiente forma:

# Número de elementos en cada nodo
nb <- table(carte$unit.classif)
print(nb)
Número de elementos en cada nodo del mapa.

Nota: los nodos del mapa de Kohonen se enumeran de izquierda a derecha y de abajo a arriba.

Vemos que el primer nodo contiene 42 elementos, el segundo nodo 275, el tercer nodo 610 elementos… y el último nodo, el 49 (7×7 celdas son 49 nodos en total) contiene 126 elementos. A simple vista se observa que no hay ningún nodo sin elementos, aunque se podría comprobar directamente de la siguiente forma:

# Comprobación de nodos vacíos
print(length(nb))

Es importante aclarar que la elección del tamaño y configuración del mapa es a criterio del data scientist. Idealmente, la distribución debe ser lo más homogénea posible. El tamaño del mapa debe reducirse si hay muchas celdas vacías o, por el contrario, debemos aumentarlo si aparecen áreas de muy alta densidad.

Podemos comprobar también a qué nodo se ha asignado a cada observación de nuestro dataframe de la siguiente forma:

# Nodo asignado a cada registro 
print(carte$unit.classif) 
Asignación de nodos.

El primer registro de nuestro dataset ha sido asignado al nodo 36, el segundo registro al nodo 25, el tercer registro al nodo 44… y así sucesivamente. Pero OJO: los nodos no son clústers, sino que los clústers serán uno o más nodos como veremos en el punto 5.

Por último procedemos a representar la gráfica de distancias. Llamada “U-Matrix” (matriz de distancia unificada), es un mapa con la distancia euclidiana entre los vectores de cada neurona con su vecina representada con una degradación de colores.

# Plot distancia con los vecinos
plot(carte, type="dist.neighbours", shape = "straight")
Mapa de distancias o U-Matrix

Los nodos que forman el mismo grupo tienden a estar cerca. Las áreas del borde están delimitadas por nodos que están lejos el uno del otro.

4. Influencia de las variables

En la primera gráfica del apartado 3 vimos como las diferentes áreas que componen el mapa topológico daban protagonismo para diferentes variables. Pero en lugar de hacer un solo gráfico para todas las variables, podemos hacer un gráfico para cada variable, tratando de resaltar los contrastes entre las áreas de alto y bajo valor. Esta descripción univariante es más fácil de entender que la anteriormente indicada:

# Mapa de densidad de colores de azul a rojo
coolBlueHotRed <- function(n, alpha = 1) {
  rainbow(n, end=4/6, alpha=alpha)[n:1]}

# Mapa de calor para cada variable
par(mfrow=c(3,3))
for (j in 1:ncol(df)){
  plot(carte,type="property",property=getCodes(carte,1)[,j], palette.name=coolBlueHotRed,main=colnames(df)[j],cex=0.5, shape = "straight")}
Mapa de calor de cada variable

Se observa la zona de alta influencia de cada variable por los colores rojo y naranja.Por ejemplo, V6 predomina en la parte derecha superior del mapa y V7 en la esquina inferior izquierda.

Podemos dar un paso más y calcular un gráfico de Variable Factor Map de forma similar a cómo se hace en el análisis de PCA. Calculamos la correlación ponderada para cada columna del vector de pesos de cada nodo:

# Correlación en función de coordenadas (x,y) del mapa
# v: variable (primera columna), w: weight, grille: Mapa de Kohonen
weighted.correlation <- function(v,w,grille){
  x <- grille$grid$pts[,"x"]
  y <- grille$grid$pts[,"y"]
  mx <- weighted.mean(x,w)
  my <- weighted.mean(y,w)
  mv <- weighted.mean(v,w)
  numx <- sum(w(x-mx)(v-mv))
  denomx <- sqrt(sum(w(x-mx)^2))sqrt(sum(w(v-mv)^2))   
  numy <- sum(w(y-my)(v-mv))   
  denomy <- sqrt(sum(w(y-my)^2))sqrt(sum(w(v-mv)^2))
  # Correlación para los dos ejes
  res <- c(numx/denomx,numy/denomy)
  return(res)
 }
 
# Correlación para cada columna del vector de nodos
CORMAP <- apply(getCodes(carte,1),2,weighted.correlation,w=nb,grille=carte)
print(CORMAP)

Obtenemos un par de valores para cada variable:

Correlación ponderada de cada columna.

Lo representamos en un diagrama de dispersión:

# Representación de cada variable en un mapa de factores
par(mfrow=c(1,1))
plot(CORMAP[1,],CORMAP[2,],xlim=c(-1,1),ylim=c(-1,1),type="n")
lines(c(-1,1),c(0,0))
lines(c(0,0),c(-1,1))
text(CORMAP[1,],CORMAP[2,],labels=colnames(Z),cex=0.75)
symbols(0,0,circles=1,inches=F,add=T)
Variables en un Variable Factor Map

El mapa de factores variables es una herramienta indicativa en este marco. Los cálculos no tienen significados matemáticos, especialmente porque se supone que el mapa topológico puede representar patrones no lineales.

Sin embargo, este tipo de gráfico sintético es más fácil de inspeccionar que las gráficas de pesos vectoriales o los mapas de calor cuando aumenta el número de variables y nodos.

En nuestro caso, vemos claramente como las variables se agrupan por pares en cada cuadrante.

Por último en esta sección, podemos calcular la importancia de cada variable a partir de la varianza ponderada según el tamaño de cada nodo, ya que las variables están estandarizadas y por tanto sus influencias son directamente comparables.

# Importancia de cada variable - varianza ponderada por el tamaño de la celda
sigma2 <- sqrt(apply(getCodes(carte,1),2,function(x,effectif)
   {m<-sum(effectif*(x-weighted.mean(x,effectif))^2)/(sum(effectif)-1)},
   effectif=nb))
 print(sort(sigma2,decreasing=T))
Relevancia de cada variable.

V4 y V7 (asimetría del perfil integrado y exceso de curtosis de la curva DM-SNR) tienen la mayor influencia mientras que V1 y V2 (media y desviación estándar del perfil integrado) son los menos relevantes del conjunto.

5. Clustering SOM

Los Mapas de Kohonen son un tipo de análisis de clústering, donde los nodos adyacentes tienen vectores de pesos similares. Podemos utilizar cada nodo del mapa como un pre-clúster para realizar un agrupamiento jerárquico (HAC):

# Matriz de distancia entre nodos
dc <- dist(getCodes(carte,1))
# HAC
cah <- hclust(dc,method="ward.D2",members=nb)
plot(cah,hang=-1,labels=F)

Nota: El parámetro “members” es crucial. De hecho, el número de instancias relacionadas con cada nodo no es el mismo, es decir, los nodos no tienen pesos idénticos. Debemos tenerlo en cuenta al calcular el método de Ward.

Dendrograma

Parece que una partición en 3, 4 ó 5 clústers es razonable. Probamos con 4:

# Visualización de los 4 clústeres
rect.hclust(cah,k=4)
Partición del dendrograma en 4 clústers.

Representamos a qué clúster pertenece cada nodo del Mapa de Kohonen:

# A qué clúster pertenece cada nodo del mapa de kohonen
groupes <- cutree(cah, k=4)
print(groupes)
Clúster al que pertenece cada nodo del mapa.

Se observa que el primer nodo queda asignado el clúster 1, el segundo nodo al clúster 2… y por último el nodo 49 al clúster 4.

Podemos ver las agrupaciones en el mapa:

# Visualización de los clústers en el mapa
plot(carte,type="mapping",bgcol=c("steelblue1","sienna1","yellowgreen","red")[groupes], 
      shape = "straight")
add.cluster.boundaries(carte,clustering=groupes)
Clusterización de los nodos.

Pasamos a asignar a cada registro del dataframe original su correspondiente clúster:

# Asignamos a cada registro su clúster
ind.groupe <- groupes[carte$unit.classif]
df$groupes <- ind.groupe

Con los datos ya clusterizados, podemos obtener datos estadísticos de cada grupo con las funciones group_by() y summarise(). La interpretación de los clústers queda a cargo de los astrónomos expertos en los datos.

6. Representación 3D de los clústeres

Dado que conocemos la pertenencia del grupo de individuos y las variables que explican la agrupación, es posible representarlos en un espacio tridimensional formado a partir de las principales variables que definen las 3 áreas identificadas previamente. Usaremos las 3 variables más significativas: V4, V7 y V8:

# Representación 3D de los puntos y su clúster asignado
points3D(df$V4, df$V7, df$V8,colvar=ind.groupe,
          col=c("steelblue1","sienna1","yellowgreen","red"),
          phi=35,theta=70)
Representación 3D de los clústers en función de las 3 variables más significativas.

Podemos obtener una representación 3D interactiva:

# Representación 3D animada
plotrgl(lighting = T)
Representación 3D interactiva.

7. Resultados finales

Por último, podemos recuperar la columna clasificadora en nuestro dataframe si fuera necesario y el nombre de las variables originales:

# Recuperamos clasificador y los nombres de las columnas
df$class <- class
col_last <- colnames(df)[9:10] # Dos últimos nombres
colnames(df) <- c(n_originales, col_last)
rm(col_last)
Dataframe final clasificado.

Para guardar el nuevo dataframe en un excel, ejecutamos lo siguiente:

library(openxlsx)
write.xlsx(df, file = "groupes.xlsx")
Para saber más:

Self-organizing map.
Representation of data using a Kohonen map, followed by a cluster analysis.
Self-Organising Maps for Customer Segmentation using R.
Los mapas auto-organizados de Kohonen (SOM).
Kohonen package | R Documentation.

8 comentarios en “Ejemplo de uso de un Mapa Auto-Organizado (SOM) de Kohonen en R”

  1. Muchas gracias por la publicación. Me gustaría hacer unas preguntas, si es posible. Puedo utilizar variables categóricas con este SOM, por ejemplo si utilizo una variable como el sexo (M/F) utilizaría variables dummy’s y, otra pregunta, podría utilizar variables de tipo factor, tengo una variable que incluye 5 categorías (A,B,C,D,E), lo cambio por (1,2,3,4,5) y lo transformo a factor, ¿Es esto posible?. Nuevamente, gracias, estoy iniciando en este tema.

Deja una respuesta