1 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.

2 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.

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:

\(\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.

3 Particionamiento espectral

El método general de agrupación espectral consiste en utilizar un método de agrupación estándar usando los vectores propios más relevantes de la matriz Laplaciana de \(\mathbf{A}\) dada por \(\mathbf{L} = \mathbf{D} - \mathbf{A}\), donde \(\mathbf{A} = [a_{i,j}]\) es una matriz de modularidad y \(\mathbf{D} = \textsf{diag}(d_1,\ldots,d_n)\) es una matriz diagonal con entradas \[ d_i = \sum_{j} a_{i,j}\,. \]

La matriz de modularidad (modularity matrix) de un grafo \(G=(V,E)\) con matriz de adyacencia \(\mathbf{Y}\) corresponde a la matriz \(\mathbf{A} = [a_{i,j}]\), donde \[ a_{i,j} = y_{i,j} - \tfrac{1}{2m}d_id_j\,. \]

Alternativamente, la matriz de modularidad también se puede definir como \[ \mathbf{A} = \mathbf{Y} - \mathbf{P}\,, \] donde \(\mathbf{P}\) contiene las probabilidades de interacción de los vértices bajo el modelo de configuración (configuration model, modelo que genera redes aleatorias de acuerdo con una secuencia grados específica).

cluster_leading_eigen en igraph.

4 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

5 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:

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:

6 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(get.adjacency(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    : int 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 <- get.adjacency(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")

7 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)))

8 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

9 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:

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.645

10 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.736         0.713         0.759         0.748
# coeficiente de variación
round(apply(X = out, MARGIN = 2, FUN = sd)/colMeans(out), 3)
##   fast_greedy leading_eigen      walktrap       louvain 
##         0.048         0.047         0.045         0.040

11 Referencias