1 Introducción

Se asume que los vértices de una red están asociados a diferentes clases, y la propensión a formar vínculos entre pares de vértices está influenciada por la clase a la que pertenecen ambos vértices.

Los vínculos surgen como resultado de la equivalencia estructural (structural equivalence), que refleja la similitud en los roles sociales desempeñados por los vértices.

El vértice \(i\in V\) del grafo \(G=(V,E)\) pertenece a una sola clase (comunidad) de una partición \(\mathcal{P} = \{C_1,\ldots,C_Q\}\) de \(V\) con \(Q\) comunidades.

El modelo de bloques estocásticos se define como:
\[ p(\mathbf{y} \mid \boldsymbol{\theta}) = \frac{1}{\kappa} \exp{\left\{ \sum_{q,r} \theta_{q,r} \, L_{q,r}(\mathbf{y}) \right\}}, \]
donde:

2 Modelo de bloques estocásticos

En la práctica, no se conoce ni la pertenencia de los vértices a las clases ni el número de clases.

La asignación de clase (class membership) de cada vértice \(i\) se determina de manera independiente mediante una distribución de probabilidad común sobre el conjunto \(\{1, \ldots, Q\}\), representada como:
\[ z_{i,q} = \begin{cases} 1 & \text{si el vértice } i \text{ pertenece a la clase } q, \\ 0 & \text{en caso contrario}. \end{cases} \]

Por lo tanto,
\[ \boldsymbol{z}_i \mid \boldsymbol{\alpha} \overset{\text{iid}}{\sim} \textsf{Multinomial}(1, \boldsymbol{\alpha}) \quad \Longleftrightarrow \quad \textsf{Pr}(z_{i,q} = 1 \mid \alpha_q) = \alpha_q, \]
donde \(\boldsymbol{z}_i = (z_{i,1}, \ldots, z_{i,Q})\) para \(i = 1, \ldots, n\), y \(\boldsymbol{\alpha} = (\alpha_1, \ldots, \alpha_Q)\) es el vector de probabilidades de clase, con la restricción \(\sum_{q=1}^Q \alpha_q = 1\).

Condicional a los valores de \(\boldsymbol{z}_1, \ldots, \boldsymbol{z}_n\), las díadas \((i,j)\) son modeladas como condicionalmente independientes, siguiendo una distribución Bernoulli:
\[ y_{i,j} \mid \mathbf{\Pi}, \boldsymbol{z}_i, \boldsymbol{z}_j \overset{\text{ind}}{\sim} \textsf{Bernoulli}(\pi_{\xi_i, \xi_j}), \]
donde \(\mathbf{\Pi} = [\pi_{q,r}]\) es una matriz \(Q \times Q\) que contiene las probabilidades de interacción entre clases, y \(\xi_i = \xi(\boldsymbol{z}_i)\) denota la posición \(q\) en \(\boldsymbol{z}_i\) tal que \(z_{i,q} = 1\) (es decir, \(\xi_i = q\) implica que el vértice \(i\) pertenece a la clase \(q\)).

En la literatura, se han propuesto diversos métodos para aproximar o mejorar la máxima verosimilitud en este contexto. El paquete blockmodels implementa uno de estos métodos, conocido como el algoritmo EM variacional.

El algoritmo EM variacional aproxima la inferencia en modelos con variables ocultas, donde calcular la verosimilitud exacta es intratable. En el Paso E, se actualiza una distribución aproximada (variacional) para las variables latentes maximizando un límite inferior de la verosimilitud (ELBO). En el Paso M, se ajustan los parámetros del modelo maximizando el mismo límite. A diferencia del EM clásico, utiliza una aproximación eficiente de la posterior en lugar de calcularla exactamente, lo que lo hace adecuado para modelos complejos o de gran escala.

3 Ejemplo: Blogs de política

Red de referencias on-line entre blogs políticos franceses clasificados por el proyecto Observatoire Presidentielle en relación con su afiliación política.

Un enlace indica que al menos uno de los dos blogs hace referencia al otro en su página web.

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

Disponible en el paquete sand de R.

suppressMessages(suppressWarnings(library(igraph)))
suppressMessages(suppressWarnings(library(sand)))
# datos
data(fblog)
# Actualizar el grafo
fblog <- igraph::upgrade_graph(fblog)
# color de vertices
cols <- RColorBrewer::brewer.pal(n = 9, name = "Set1")
party.nums <- as.numeric(as.factor(V(fblog)$PolParty))
# grafico
par(mfrow = c(1,1), mar = 0.2*c(1,1,1,1))
set.seed(42)
plot(fblog, layout = layout_with_fr, vertex.label = NA, vertex.size = 4, vertex.color = cols[party.nums], vertex.frame.color = cols[party.nums], edge.color = adjustcolor("black",0.1))

# orden
vcount(fblog)
## [1] 192
# tamaño
ecount(fblog)
## [1] 1431
# dirigida?
is_directed(fblog)
## [1] FALSE
# paquete para ajustar SBMs
suppressMessages(suppressWarnings(library(blockmodels)))
# matriz de adyacencia
A <- as.matrix(igraph::as_adjacency_matrix(fblog))
dim(A)
## [1] 192 192
# formulación del modelo
# membership_type -> "SBM_sym" para redes no dirigidas 
# membership_type -> "SBM" para redes dirigidas
set.seed(42)
fblog.sbm <- blockmodels::BM_bernoulli(membership_type = "SBM_sym", adj = A, verbosity = 0, plotting = "")
# estimación
fblog.sbm$estimate()

La verosimilitud de clasificación de integración (ICL, Integrated Classification Likelihood) es un criterio diseñado específicamente para problemas de agrupamiento, que sigue un enfoque similar al de los criterios de información empleados en modelos estándar, como el criterio de información de Akaike (AIC) y el criterio de información bayesiano (BIC).

La ICL considera explícitamente la asignación de las clases en modelos de agrupamiento, lo que lo hace especialmente adecuado para problemas de agrupamiento o modelos de mezcla: \[ \text{ICL} = \log p(\mathbf{y}, \hat{\mathbf{z}} \mid \hat{\boldsymbol{\theta}}) - \frac{\nu}{2} \log n, \] donde:

El término \(\log p(\mathbf{y}, \hat{\mathbf{z}} \mid \hat{\boldsymbol{\theta}})\) mide la calidad del ajuste del modelo, considerando tanto los datos como las asignaciones de clase. El término de penalización \(-\frac{\nu}{2} \log n\) controla el sobreajuste, penalizando modelos más complejos en función del número de parámetros y el tamaño de los datos.

# ICL
(ICLs <- fblog.sbm$ICL)
##  [1] -5028.312 -4580.705 -4318.101 -4118.934 -3986.373 -3865.908 -3810.839
##  [8] -3768.287 -3742.687 -3720.026 -3725.161 -3752.341 -3783.983 -3833.919
## [15] -3886.516
# número de grupos optimo
(Q <- which.max(ICLs))
## [1] 10
# gráfico del ICL
par(mfrow = c(1,1), mar = c(2.75,2.75,1.5,0.5), mgp = c(1.7,0.7,0))
plot(fblog.sbm$ICL, xlab = "Q", ylab = "ICL", type = "b", pch = 16)
lines(x = c(Q,Q), y = c(min(ICLs),max(ICLs)), col = "red", lty = 2)

Es posible estimar las probabilidades de pertenencia a las comunidades, es decir, los valores esperados de \(\boldsymbol{z}_i\) condicionados a \(\mathbf{Y} = \mathbf{y}\).

Estas estimaciones permiten determinar las etiquetas de asignación a las comunidades (cluster assignments o class memberships), utilizando un criterio basado en la probabilidad máxima. En este caso, se asigna un vértice a la comunidad correspondiente si su probabilidad de pertenencia supera un umbral predefinido (por ejemplo, 0.8586).

# probabilidades estimadas de pertenencia a las comunidades
Z <- fblog.sbm$memberships[[Q]]$Z
head(x = round(Z,3), n = 10)
##        [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
##  [1,] 0.995 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
##  [2,] 0.995 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
##  [3,] 0.995 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
##  [4,] 0.995 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
##  [5,] 0.995 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
##  [6,] 0.995 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
##  [7,] 0.995 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
##  [8,] 0.001 0.001 0.995 0.001 0.001 0.001 0.001 0.001 0.001 0.001
##  [9,] 0.001 0.001 0.995 0.001 0.001 0.001 0.001 0.001 0.001 0.001
## [10,] 0.001 0.001 0.001 0.001 0.001 0.995 0.001 0.001 0.001 0.001
tail(x = round(Z,3), n = 10)
##         [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
## [183,] 0.001 0.001 0.001 0.001 0.995 0.001 0.001 0.001 0.001 0.001
## [184,] 0.001 0.001 0.001 0.001 0.995 0.001 0.001 0.001 0.001 0.001
## [185,] 0.001 0.001 0.001 0.001 0.995 0.001 0.001 0.001 0.001 0.001
## [186,] 0.001 0.001 0.001 0.001 0.995 0.001 0.001 0.001 0.001 0.001
## [187,] 0.001 0.001 0.001 0.001 0.995 0.001 0.001 0.001 0.001 0.001
## [188,] 0.001 0.001 0.001 0.001 0.995 0.001 0.001 0.001 0.001 0.001
## [189,] 0.001 0.001 0.001 0.001 0.995 0.001 0.001 0.001 0.001 0.001
## [190,] 0.001 0.001 0.001 0.001 0.995 0.001 0.001 0.001 0.001 0.001
## [191,] 0.001 0.001 0.001 0.001 0.995 0.001 0.001 0.001 0.001 0.001
## [192,] 0.001 0.001 0.001 0.001 0.995 0.001 0.001 0.001 0.001 0.001
# dimensión
dim(Z)
## [1] 192  10
# asignaciones
labs <- apply(X = Z, MARGIN = 1, FUN = which.max)
head(x = labs, n = 10)
##  [1] 1 1 1 1 1 1 1 3 3 6
tail(x = labs, n = 10)
##  [1] 5 5 5 5 5 5 5 5 5 5
length(labs)
## [1] 192
# resumen de las probabilidades maximales
summary(Z[cbind(1:vcount(fblog), labs)])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.8586  0.9953  0.9953  0.9938  0.9953  0.9953
# tamaño de las comunidades
table(labs)
## labs
##  1  2  3  4  5  6  7  8  9 10 
## 35 27 11 21 24 25  6  7 13 23
# probabilidades de los grupos
alpha <- table(labs)/vcount(fblog)
round(alpha, 3)
## labs
##     1     2     3     4     5     6     7     8     9    10 
## 0.182 0.141 0.057 0.109 0.125 0.130 0.031 0.036 0.068 0.120
# probabilidades de los grupos (ordenadas)
round(alpha[order(alpha, decreasing = T)], 3)
## labs
##     1     2     6     5    10     4     9     3     8     7 
## 0.182 0.141 0.130 0.125 0.120 0.109 0.068 0.057 0.036 0.031

Probabilidades de interacción (intection probabilities) \(\mathbf{\Pi}=[\pi_{q,r}]\):

# matriz de probabilidades de interacción
Pi <- fblog.sbm$model_parameters[[Q]]$pi
# gráfico
corrplot::corrplot(corr = Pi, type = "full", col.lim = c(0,1),  method = "shade", addgrid.col = "gray90", tl.col = "black")

Grafo y matriz de adyacencia:

# funciones
# para ordenar la matriz de adyacencia respecto a una partición
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)
}
# para graficar la matriz de adyacencia
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 de grupos
xi <- apply(X = Z, MARGIN = 1, FUN = which.max)
# matriz de adyacencia original
Y <- A
# matriz de adyacencia ordenada y lineas divisorias de acuerdo con las comunidades
tmp <- get_adjacency_ordered(xi = xi, A = Y)
# viz
par(mfrow = c(1,2), mar = 0.2*c(1,1,1,1))
# G
cols <- RColorBrewer::brewer.pal(n = 12, name = "Paired")
set.seed(42)
plot(fblog, layout = layout_with_fr, vertex.label = NA, vertex.size = 5, vertex.color = cols[labs], vertex.frame.color = cols[labs], edge.color = adjustcolor("black",0.1))
# A
heat.plot0(mat = tmp$A, tick = F, labs = NA)
abline(v = tmp$d+.5, h = tmp$d+.5)

Comparación con el agrupamiento natural:

# agrupamiento natural
party.nums <- as.numeric(as.factor(V(fblog)$PolParty))
# agrupamiento estimado
labs <- apply(X = Z, MARGIN = 1, FUN = which.max)
# comparación
round(igraph::compare(comm1 = party.nums, comm2 = labs, method = "rand"), 4)
## [1] 0.8619
round(igraph::compare(comm1 = party.nums, comm2 = labs, method = "adjusted.rand"), 4)
## [1] 0.4608
round(igraph::compare(comm1 = party.nums, comm2 = labs, method = "nmi"), 4)
## [1] 0.6565
table(party.nums, labs)
##           labs
## party.nums  1  2  3  4  5  6  7  8  9 10
##          1  2  0  0  0  0  0  0  0  0  0
##          2  2  1  0  0  0  0  0  1  7  0
##          3  7  0  0  0  0  0  0  0  0  0
##          4  1  0  0  0 24  0  0  0  0  0
##          5  9  0  0  0  0  0  0  0  1  1
##          6  7  0  0  0  0  0  0  0  0  0
##          7  6  0  0 21  0  0  6  0  2 22
##          8  0 24  0  0  0  1  0  6  1  0
##          9  1  2 11  0  0 24  0  0  2  0

Para realizar la bondad de ajuste del modelo se pueden usar estadísticos de prueba por medio de métodos de simulación.

# Cargar librerías necesarias
suppressMessages(suppressWarnings(library(knitr)))
suppressMessages(suppressWarnings(library(kableExtra)))

# Calcular las estadísticas observadas
observed_density <- edge_density(fblog)
observed_transitivity <- transitivity(fblog, type = "global")
observed_assortativity <- assortativity_degree(fblog)

# Probabilidades de interacción
labs <- apply(X = Z, MARGIN = 1, FUN = which.max)
alpha <- table(labs) / vcount(fblog)  # Distribución de comunidades
Pi <- fblog.sbm$model_parameters[[Q]]$pi
Pi <- 0.5*(t(Pi) + Pi)  # asegurarse que Pi sea simétrica

# Inicializar vectores para estadísticas simuladas
n_sim <- 1000
simulated_densities <- numeric(n_sim)
simulated_transitivities <- numeric(n_sim)
simulated_assortativities <- numeric(n_sim)

# Simular redes y calcular estadísticas
set.seed(42)
for (i in 1:n_sim) {
  # Generar red simulada
  bs <- stats::rmultinom(n = 1, size = vcount(fblog), prob = alpha)
  g_sim <- igraph::sample_sbm(n = vcount(fblog), pref.matrix = Pi, block.sizes = bs, directed = FALSE)
  
  # Calcular estadísticas de la red simulada
  simulated_densities[i] <- edge_density(g_sim)
  simulated_transitivities[i] <- transitivity(g_sim, type = "global")
  simulated_assortativities[i] <- assortativity_degree(g_sim)
}

# Calcular intervalos de confianza al 95% basados en percentiles
ci_lower <- c(
  quantile(simulated_densities, 0.025),
  quantile(simulated_transitivities, 0.025),
  quantile(simulated_assortativities, 0.025)
)

ci_upper <- c(
  quantile(simulated_densities, 0.975),
  quantile(simulated_transitivities, 0.975),
  quantile(simulated_assortativities, 0.975)
)

# Crear la tabla de resultados
resultados <- data.frame(
  Estadística = c("Densidad", "Transitividad", "Asortatividad"),
  Observado = c(observed_density, observed_transitivity, observed_assortativity),
  Media_Simulada = c(mean(simulated_densities), mean(simulated_transitivities), mean(simulated_assortativities)),
  Desviación_Simulada = c(sd(simulated_densities), sd(simulated_transitivities), sd(simulated_assortativities)),
  IC_Inferior = ci_lower,
  IC_Superior = ci_upper
)

# Presentar la tabla estilizada
resultados %>%
  kbl(
    caption = "Resumen de las estadísticas observadas y simuladas con intervalos de confianza al 95%",
    digits = 4,
    col.names = c("Estadística", "Observado", "Media", "Desviación", "IC Inferior", "IC Superior")
  ) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) %>%
  column_spec(1, bold = TRUE) %>%
  add_header_above(c(" " = 1, "Estadísticas Observadas y Simuladas" = 5))
Resumen de las estadísticas observadas y simuladas con intervalos de confianza al 95%
Estadísticas Observadas y Simuladas
Estadística Observado Media Desviación IC Inferior IC Superior
Densidad 0.0780 0.0795 0.0068 0.0671 0.0943
Transitividad 0.3858 0.3162 0.0299 0.2595 0.3777
Asortatividad 0.0124 0.0870 0.0584 -0.0080 0.2161

Referencias