Introducción

El agrupamiento (clustering/partitioning) se refiere a la segmentación de un conjunto de elementos en subconjuntos “naturales”.

Una partición \(\mathcal{C}\) de un conjunto finito \(S\) se refiere a una descomposición de \(S\) en \(K\) subconjuntos \(C_1,\ldots,C_K\) de \(S\) tales que: \[ C_k\neq\Phi\,,\qquad C_k\cap C_\ell = \Phi\,,\qquad \cup_{k=1}^K C_k = S\,, \] para \(k,\ell\in\{1,\ldots,K\}\) con \(k\neq \ell\).

La partición de redes (network partitioning) también conocida como detección de comunidades (commuty detection), constituye una metodología no supervisada para encontrar subconjuntos de vértices “homogéneos” respeto a sus patrones relacionales.

Los algoritmos de agrupamiento de grafos establecen una partición \(\mathcal{C}=\{C_1,\ldots,C_K\}\) del conjunto de vértices \(V\) de un grafo \(G=(V,E)\) tal que el conjunto de aristas conectando vértices de \(C_k\) con vértices de \(C_{\ell}\) sea relativamente “pequeño” comparado con el conjunto de aristas conectando vértices dentro de \(C_k\).

Los métodos que se discuten a continuación son algorítmicos y no constituyen modelos estadísticos, y por lo tanto, no permiten cuantificar la incertidumbre asociada con la detección de las comunidades como la formación de enlaces.

Agrupación jerárquica

Estos métodos adoptan un enfoque computacional intensivo para explorar el espacio de todas las posibles particiones, modificando iterativamente las particiones candidatas.

Los métodos difieren principalmente en la métrica para evaluar la calidad de las agrupaciones y los algoritmos para optimizar tal métrica.

  • Agrupamiento aglomerativo (agglomerative clustering): ampliación de las comunidades a través de fusión (merging).
  • Agrupamiento divisivo (divisive clustering): refinamiento de las comunidades a través de división (spliting).

En cada etapa del algoritmo, la partición actual se modifica de manera que se minimice/maximice una función de perdida/utilidad específica mediante la acción menos costosa de fusión/división.

La modulariad de una red respecto a una partición mide qué tan buena es la división o qué tan separados están los diferentes tipos de vértices entre sí: \[ \textsf{mod}(\mathcal{C}) = \frac{1}{2m} \sum_{i,j:i\neq j} \left(y_{i,j} - \tfrac{1}{2m}d_id_j\right)\delta_{c_i,c_j}\,, \] donde:

  • \(\mathbf{Y}=[y_{i,j}]\) : matriz de adyacencia.
  • \(m\) : tamaño del grafo (número de aristas).
  • \(d_i\) : grado del vértice \(i\).
  • \(c_i\) : grupo de la partición al que pertenece el vértice \(i\) (membership o cluster assignment).
  • \(\delta_{x,y} = 1\) si \(x=y\) y \(\delta_{x,y} = 0\) en otro caso.

\(\textsf{mod}(\mathcal{C}) = 0\) si la formación de aristas ocurre sin tener en cuenta las comunidades de la partición de referencia.

Valores más grandes de \(\textsf{mod}(\mathcal{C})\) indican una estructura de comunidades más fuerte.

cluster_fast_greedy en igraph.

Clauset, A., Newman, M. E., & Moore, C. (2004). Finding community structure in very large networks. Physical review E, 70(6), 066111.

Particionamiento espectral

Este método consiste en aplicar un método de agrupación estándar (e.g., \(k\)-medias), a una representación de los vértices construida a partir de los vectores propios más relevantes de la matriz laplaciana del grafo.

Si \(\mathbf{Y} = [y_{i,j}]\) es la matriz de adyacencia de un grafo no dirigido, la matriz laplaciana se define como \[ \mathbf{L} = \mathbf{D} - \mathbf{Y}, \] donde \(\mathbf{D} = \textsf{diag}(d_1,\ldots,d_n)\) es la matriz diagonal de grados, con entradas \[ d_i = \sum_j y_{i,j}. \]

La idea central de este enfoque es que la estructura espectral de \(\mathbf{L}\) contiene información sobre la conectividad del grafo, lo que permite identificar grupos de vértices fuertemente relacionados entre sí y débilmente conectados con el resto de la red.

Por otra parte, también es común utilizar la matriz de modularidad. Si \(\mathbf{Y}\) es la matriz de adyacencia del grafo, la matriz de modularidad se define como \(\mathbf{B} = [b_{i,j}]\), donde \[ b_{i,j} = y_{i,j} - \frac{1}{2m}\,d_i d_j, \] y \(m\) denota el número total de aristas del grafo.

La matriz laplaciana y la matriz de modularidad son objetos diferentes. La primera se utiliza principalmente en partición espectral y agrupación espectral, mientras que la segunda se emplea sobre todo en métodos de detección de comunidades basados en maximización de modularidad.

Cúales vectores propios elegir?

Con la matriz laplaciana \(\mathbf{L}\), los vectores propios más relevantes son los asociados a los autovalores más pequeños (la estructura de comunidades o particiones “suaves” del grafo aparece en la parte baja del espectro de \(\mathbf{L}\)). En general, se descarta el vector trivial asociado a \(\lambda_1=0\) y, para identificar \(K\) grupos, se usan los vectores correspondientes a \(\lambda_2,\ldots,\lambda_K\).

Con la matriz de modularidad \(\mathbf{B}\), los vectores propios más relevantes suelen ser los asociados a los autovalores más grandes, especialmente los positivos (indican direcciones en las que la estructura de comunidades se aparta más de lo esperado bajo el modelo nulo).

En ambos casos, una guía útil para elegir cuántos vectores usar es observar el eigengap, es decir, un salto grande entre autovalores consecutivos.

cluster_leading_eigen en igraph.

Resumen

Algoritmo Función en igraph Idea
Fast-greedy cluster_fast_greedy Optimizar una métrica de modularidad
Edge-betweenness cluster_edge_betweenness Optimizar una métrica de aristas basada en caminos más cortos
Leading eigenvector cluster_leading_eigen Calcular el vector propio principal no negativo de la matriz de modularidad
Louvain cluster_louvain Optimizar una métrica de modularidad múltinivel
Walktrap cluster_walktrap Caminatas aleatorias cortas tienden a permanecer en la misma comunidad
Label propagation cluster_label_prop Etiquetar vértices con etiquetas únicas y actualizarlas por votación mayoritaria en la vecindad del vértice
InfoMAP cluster_infomap Optimizar la longitud esperada de una trayectoria de una caminata aleatoria
Spinglass cluster_spinglass Modelo spin-glass y simulated annealing
Optimal cluster_optimal Optimizar una métrica de modularidad

Validación

Cuando se tiene conocimiento de alguna noción de pertenencia a una clase definida externamente, resulta interesante interesante comparar y contrastar las asignaciones resultantes con las que se derivan de la partición.

compare en igraph.

Comparación de dos particiones:

  • rand : the Rand index (Rand 1971).
  • adjusted.rand : adjusted Rand index (Hubert and Arabie 1985).
  • vi : variation of information (VI) metric (Meila 2003).
  • nmi : normalized mutual information measure (Danon et al. 2005).
  • split.join : split-join distance (can Dongen 2000).

Por ejemplo, el índice Rand tiene un valor entre 0 y 1, donde 0 indica que las dos agrupaciones de datos no coinciden en ningún par de puntos y 1 indica que las agrupaciones de datos son exactamente iguales: \[ \textsf{RI}(X,Y) = \frac{a+b}{a+b+c+d} \] donde:

  • \(X\) y \(Y\): particiones del conjunto de números enteros \(S=\{1,\ldots,n\}\).
  • \(a\): el número de pares de elementos en \(S\) que están en el mismo subconjunto en \(X\) y en el mismo subconjunto en \(Y\).
  • \(b\): el número de pares de elementos en \(S\) que están en diferentes subconjuntos en \(X\) y en diferentes subconjuntos en \(Y\).
  • \(c\): el número de pares de elementos en \(S\) que están en el mismo subconjunto en \(X\) y en diferentes subconjuntos en \(Y\).
  • \(d\): el número de pares de elementos en \(S\) que están en diferentes subconjuntos en \(X\) y en el mismo subconjunto en \(Y\).

Ejemplo: Interacciones sociales

Red de interacciones sociales entre los miembros de un club de karate.

Estos datos fueron recolectados para estudiar la fragmentación que sufrió el club en dos clubes diferentes debido a una disputa entre el director y el administrador.

\(y_{i,j} = 1\) si los miembros \(i\) y \(j\) tuvieron una interacción social en el club y \(y_{i,j} = 0\) en otro caso.

Una descripción completa de los datos se puede encontrar aquí.

Disponible en el paquete igraphdata de R.

Zachary, W. W. (1977). An information flow model for conflict and fission in small groups. Journal of anthropological research, 33(4), 452-473.

# librerías
suppressMessages(suppressWarnings(library(igraph)))
suppressMessages(suppressWarnings(library(igraphdata)))
suppressMessages(suppressWarnings(library(sand)))
# datos
data(karate)
karate <- upgrade_graph(karate)
# orden
vcount(karate)
## [1] 34
# tamaño
ecount(karate)
## [1] 78
# dirigida?
is_directed(karate)
## [1] FALSE
# ponderada?
is_weighted(karate)
## [1] TRUE
# visualización
par(mfrow = c(1,2))

set.seed(123)
l = layout_with_dh(karate)
plot(karate, 
     layout = l, 
     vertex.size = 14, 
     vertex.frame.color = "black", 
     vertex.label.color = "black", 
     main = "")

corrplot::corrplot(corr = as.matrix(as_adjacency_matrix(graph = karate, names = F)), 
                   col.lim = c(0,1), 
                   method = "color", 
                   tl.col = "black", 
                   addgrid.col = "gray", 
                   cl.pos = "n")

# aplicación del algoritmo
kc <- cluster_fast_greedy(karate)
# estructura
str(kc)
## Class 'communities'  hidden list of 6
##  $ merges    : num [1:33, 1:2] 26 7 17 14 3 4 8 1 30 25 ...
##  $ modularity: num [1:34] -0.0511 -0.02356 -0.00362 0.02084 0.03786 ...
##  $ membership: num [1:34] 2 2 2 2 3 3 3 2 1 1 ...
##  $ names     : chr [1:34] "Mr Hi" "Actor 2" "Actor 3" "Actor 4" ...
##  $ algorithm : chr "fast greedy"
##  $ vcount    : num 34
# algoritmo
kc$algorithm
## [1] "fast greedy"
# algoritmo jerárquico?
is_hierarchical(kc)
## [1] TRUE
# asignaciones
kc$membership
##  [1] 2 2 2 2 3 3 3 2 1 1 3 2 2 2 1 1 3 2 1 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1
# tamaños
sizes(kc)
## Community sizes
##  1  2  3 
## 18 11  5
table(kc$membership)
## 
##  1  2  3 
## 18 11  5
# tamaño de la partición
length(kc)
## [1] 3
length(table(kc$membership))
## [1] 3
# visualización
par(mfrow = c(1,2))
set.seed(1)
plot(x = kc, y = karate, vertex.size = 12)
set.seed(1)
plot(x = kc, y = karate, mark.groups = NULL, edge.color = "darkgray", vertex.size = 12)

# visualización
par(mfrow = c(1,2))
plot(karate, 
     layout = l,
     vertex.size = 12, 
     vertex.frame.color = "black", 
     vertex.label.color = "black", 
     main = "Facciones")

plot(karate, 
     layout = l, 
     vertex.size = 12, 
     vertex.frame.color = "black", 
     vertex.label.color = "black", 
     vertex.color = kc$membership, 
     main = "Fast greedy")

# funciones
get_adjacency_ordered <- function(xi, A) {
  xi2 <- xi[order(xi)]
  indices <- order(xi)
  d <- NULL
  for (i in 1:(length(xi)-1)) if (xi2[i] != xi2[i+1]) d <- c(d, i)
  list(A = A[indices,indices], d = d)
}

heat.plot0 <- function (mat, show.grid = FALSE, cex.axis, tick, labs, col.axis, ...) { 
  JJ <- dim(mat)[1]
  colorscale <- c("white", rev(heat.colors(100)))
  if(missing(labs))     labs <- 1:JJ
  if(missing(col.axis)) col.axis <- rep("black", JJ)
  if(missing(cex.axis)) cex.axis <- 1
  if(missing(tick))     tick <- TRUE
  # adjacency matrix
  image(seq(1, JJ), seq(1, JJ), mat, axes = FALSE, xlab = "", ylab = "", col = colorscale[seq(floor(100*min(mat)), floor(100*max(mat)))], ...)
  for (j in 1:JJ) {
    axis(1, at = j, labels = labs[j], las = 2, cex.axis = cex.axis, tick, col.axis = col.axis[j], col.ticks = col.axis[j])
    axis(2, at = j, labels = labs[j], las = 2, cex.axis = cex.axis, tick, col.axis = col.axis[j], col.ticks = col.axis[j])
  }
  box()
  if (show.grid) grid(nx = JJ, ny = JJ)
}

# asignaciones 
xi <- kc$membership
# asignaciones ordenadas 
xi2 <- xi[order(xi)]
# matriz de adyacencia original
Y <- as_adjacency_matrix(graph = karate, sparse = F)
# matriz de adyacencia ordenada y lineas divisorias de acuerdo con las comunidades
tmp <- get_adjacency_ordered(xi = xi, A = Y)
A <- tmp$A
d <- tmp$d
# visualización
par(mfrow = c(1,2))
heat.plot0(mat = Y)
heat.plot0(mat = A, col.axis = c("darkgoldenrod3","deepskyblue3","forestgreen")[xi2], labs = (1:vcount(karate))[order(xi)])
abline(v = d+.5, h = d+.5)

# visualización
plot_dendrogram(x = kc, mode = "phylo")

Ejemplo: Interacciones sociales (cont.)

# algoritmos
kc_fast_greedy <- cluster_fast_greedy(karate)
kc_leading_eigen <- cluster_leading_eigen(karate)
kc_walktrap <- cluster_walktrap(karate)
kc_louvain <- cluster_louvain(karate)
kc_label_prop <- cluster_label_prop(karate)
kc_spinglass <- cluster_spinglass(karate)
kc_optimal <- cluster_optimal(karate)
kc_infomap <- cluster_infomap(karate)
# gráficos
igraph_options(vertex.size = 10, vertex.frame.color = "black")
par(mfrow = c(3,3))
plot(karate, vertex.label = NA, layout = l, vertex.color = kc_fast_greedy$membership, main = paste0("fast greedy: ", "Mod = ", round(modularity(kc_fast_greedy), 4)))
plot(karate, vertex.label = NA, layout = l, vertex.color = kc_leading_eigen$membership, main = paste0("leading eigen: ", "Mod = ", round(modularity(kc_leading_eigen), 4)))
plot(karate, vertex.label = NA, layout = l, vertex.color = kc_walktrap$membership, main = paste0("walktrap: ", "Mod = ", round(modularity(kc_walktrap), 4)))
plot(karate, vertex.label = NA, layout = l, vertex.color = kc_louvain$membership, main = paste0("louvain: ", "Mod = ", round(modularity(kc_louvain), 4)))
plot(karate, vertex.label = NA, layout = l, vertex.color = kc_label_prop$membership, main = paste0("label prop: ", "Mod = ", round(modularity(kc_label_prop), 4)))
plot(karate, vertex.label = NA, layout = l, vertex.color = kc_spinglass$membership, main = paste0("spinglass: ", "Mod = ", round(modularity(kc_spinglass), 4)))
plot(karate, vertex.label = NA, layout = l, vertex.color = kc_optimal$membership, main = paste0("optimal: ", "Mod = ", round(modularity(kc_optimal), 4)))
plot(karate, vertex.label = NA, layout = l, vertex.color = kc_infomap$membership, main = paste0("infomap: ", "Mod = ", round(modularity(kc_infomap), 4)))

Ejemplo: Interacciones sociales (cont.)

fc <- as.numeric(as.factor(vertex_attr(graph = karate, name = "Faction")))
table(fc)
## fc
##  1  2 
## 16 18
# aplicación de algoritmos
kc_fast_greedy <- cluster_fast_greedy(karate)
kc_leading_eigen <- cluster_leading_eigen(karate)
# comparación
compare(comm1 = fc, comm2 = kc_fast_greedy$membership, method = "rand")
## [1] 0.9019608
compare(comm1 = fc, comm2 = kc_leading_eigen$membership, method = "rand")
## [1] 0.7557932
# tablas cruzadas
table(fc, kc_fast_greedy$membership)
##    
## fc   1  2  3
##   1  0 11  5
##   2 18  0  0
table(fc, kc_leading_eigen$membership)
##    
## fc   1  2  3  4  5
##   1 10  0  5  1  0
##   2  0 12  0  0  6

Ejemplo: Simulación

Considere el modelo de bloques aleatorios dado por \[ y_{i,j}\mid\xi_i,\xi_j,\Theta \sim\textsf{Ber}(\theta_{\phi(\xi_i,\xi_j)})\,,\qquad \text{para} 1\leq i<j\leq n\,, \] donde:

  • \(n\) : número de nodos de la red.
  • \(\xi_i\) : indicadora de cluster del nodo \(i\), i.e., \(\xi_i=k\) indica que el nodo \(i\) hace parte del cluster \(k\), para \(k=1,\ldots,K\).
  • \(\Theta\) : matriz simétrica de \(K\times K\) constituida por las probabilidades de interacción intra e inter cluster.
  • \(\phi(x,y)\) : función de valor vectorial definida por \(\phi(x,y)= (\min(x,y),\max(x,y))\).

A continuación se genera una red de acuerdo con este modelo con grupos de aproximadamente el mismo tamaño usando \(n = 100\), \(K = 4\) y \[ \Theta = \begin{bmatrix} 0.05 & 0.05 & 0.05 & 0.05 \\ 0.05 & 0.10 & 0.05 & 0.05 \\ 0.05 & 0.05 & 0.20 & 0.05 \\ 0.05 & 0.05 & 0.05 & 0.40 \\ \end{bmatrix}\,. \]

# funciones
phi <- function(x, y) c(min(x,y), max(x, y))

sim_blocks <- function (n, p, Theta) {
     # simulación redes no dirigidas (modelo de bloques estocásticos)
     K  <- length(p) 
     xi <- sample(x = 1:K, size = n, replace = T, prob = p)
     Y  <- matrix(data = 0, nrow = n, ncol = n)
     for (i in 1:(n-1)) {
          for (j in (i+1):n) {
               index <- phi(xi[i], xi[j])
               Y[i,j] <- Y[j,i] <- rbinom(n = 1, size = 1, prob = Theta[index[1], index[2]])
          }
     }
     list(Y = Y, xi = xi)
}
# argumentos
n <- 100 
p <- rep(0.25, 4)
Theta <- matrix(0.05, 4, 4)
diag(Theta) <- c(0.05, 0.1, 0.2, 0.4)
# simulación
set.seed(1)
tmp <- sim_blocks(n, p, Theta)
Y  <- tmp$Y
xi <- tmp$xi
# asignaciones ordenadas 
xi2 <- xi[order(xi)]
# matriz de adyacencia ordenada y lineas divisorias de acuerdo con las comunidades
tmp <- get_adjacency_ordered(xi = xi, A = Y)
Y <- tmp$A
d <- tmp$d
g <- graph_from_adjacency_matrix(adjmatrix = Y, mode = "undirected")
# visualización
par(mfrow = c(1,2))
set.seed(123)
plot(g, layout = layout_with_fr, vertex.label = NA, vertex.size = 6, vertex.frame.color = xi2, vertex.color = xi2, main = "")
heat.plot0(mat = Y, labs = F, tick = F)
abline(v = d+.5, h = d+.5)

Ahora, con el fin de establecer la idoneidad de los métodos de agrupamiento, se agrupa la red simulada utilizando algunos métodos de agrupamiento y luego se compara la partición real con la partición recuperado por medio de los métodos:

# agrupamiento
c1 <- cluster_fast_greedy  (g)
c2 <- cluster_leading_eigen(g)
c3 <- cluster_walktrap     (g)
c4 <- cluster_louvain      (g)
# comparación usando el índice de rand
out <- c(compare(comm1 = xi, comm2 = c1$membership, method = "rand"),
         compare(comm1 = xi, comm2 = c2$membership, method = "rand"),
         compare(comm1 = xi, comm2 = c3$membership, method = "rand"),
         compare(comm1 = xi, comm2 = c4$membership, method = "rand"))
names(out) <- c("fast_greedy","leading_eigen","walktrap","louvain")
round(out, 3)
##   fast_greedy leading_eigen      walktrap       louvain 
##         0.632         0.579         0.630         0.642

Ejemplo: Simulación (cont.)

A continuación se remite el experimento del ejemplo anterior \(N = 100\) veces:

# experimento
N <- 100
out <- matrix(NA, nrow = N, ncol = 4)
colnames(out) <- c("fast_greedy","leading_eigen","walktrap","louvain")
set.seed(123)
for (i in 1:N) {
     # simular datos
     tmp <- sim_blocks(n, p, Theta)
     g   <- graph_from_adjacency_matrix(adjmatrix = tmp$Y, mode = "undirected")
     xi  <- tmp$xi
     # agrupamiento
     c1 <- cluster_fast_greedy  (g)
     c2 <- cluster_leading_eigen(g)
     c3 <- cluster_walktrap     (g)
     c4 <- cluster_louvain      (g)
     # comparación
     out[i,] <- c(compare(comm1 = xi, comm2 = c1$membership, method = "rand"),
                  compare(comm1 = xi, comm2 = c2$membership, method = "rand"),
                  compare(comm1 = xi, comm2 = c3$membership, method = "rand"),
                  compare(comm1 = xi, comm2 = c4$membership, method = "rand"))
}
# promedio de índice de rand
round(colMeans(out), 3)
##   fast_greedy leading_eigen      walktrap       louvain 
##         0.738         0.715         0.774         0.749
# coeficiente de variación
round(apply(X = out, MARGIN = 2, FUN = sd)/colMeans(out), 3)
##   fast_greedy leading_eigen      walktrap       louvain 
##         0.044         0.049         0.042         0.037